1 |
/* $Id:$ */ |
2 |
|
3 |
/******************************************************* |
4 |
* |
5 |
* Copyright 2007 by University of Queensland |
6 |
* |
7 |
* http://esscc.uq.edu.au |
8 |
* Primary Business: Queensland, Australia |
9 |
* Licensed under the Open Software License version 3.0 |
10 |
* http://www.opensource.org/licenses/osl-3.0.php |
11 |
* |
12 |
*******************************************************/ |
13 |
|
14 |
/**************************************************************/ |
15 |
|
16 |
/* Paso: FCTransportProblem */ |
17 |
|
18 |
/**************************************************************/ |
19 |
|
20 |
/* Author: l.gross@uq.edu.au */ |
21 |
|
22 |
/**************************************************************/ |
23 |
|
24 |
|
25 |
#include "Paso.h" |
26 |
#include "SolverFCT.h" |
27 |
#include "PasoUtil.h" |
28 |
|
29 |
|
30 |
/**************************************************************/ |
31 |
|
32 |
/* free all memory used by */ |
33 |
void Paso_FCTransportProblem_free(Paso_FCTransportProblem* in) { |
34 |
if (in!=NULL) { |
35 |
in->reference_counter--; |
36 |
if (in->reference_counter<=0) { |
37 |
|
38 |
Paso_SystemMatrix_free(in->transport_matrix); |
39 |
Paso_SystemMatrix_free(in->flux_matrix); |
40 |
Paso_MPIInfo_free(in->mpi_info); |
41 |
|
42 |
MEMFREE(in->u); |
43 |
MEMFREE(in->lumped_mass_matrix); |
44 |
MEMFREE(in->row_sum_flux_matrix); |
45 |
MEMFREE(in->colorOf); |
46 |
MEMFREE(in->main_iptr); |
47 |
MEMFREE(in); |
48 |
} |
49 |
} |
50 |
} |
51 |
|
52 |
Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) { |
53 |
if (in!=NULL) { |
54 |
++(in->reference_counter); |
55 |
} |
56 |
} |
57 |
|
58 |
Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) { |
59 |
return in->transport_matrix; |
60 |
} |
61 |
|
62 |
Paso_SystemMatrix* Paso_FCTransportProblem_borrowFluxMatrix(Paso_FCTransportProblem* in) { |
63 |
return in->flux_matrix; |
64 |
} |
65 |
|
66 |
double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) { |
67 |
return in->lumped_mass_matrix; |
68 |
} |
69 |
|
70 |
dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) { |
71 |
return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix); |
72 |
} |
73 |
|
74 |
Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, double dt_max, Paso_SystemMatrixPattern *pattern, int block_size |
75 |
|
76 |
|
77 |
) { |
78 |
Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */ |
79 |
Paso_FCTransportProblem* out=NULL; |
80 |
dim_t n,i; |
81 |
index_t iptr,iptr_main,k; |
82 |
|
83 |
if ((theta<0.) || (theta >1.)) { |
84 |
Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1."); |
85 |
return NULL; |
86 |
} |
87 |
|
88 |
out=MEMALLOC(1,Paso_FCTransportProblem); |
89 |
if (Paso_checkPtr(out)) return NULL; |
90 |
|
91 |
out->theta=theta; |
92 |
out->valid_matrices=FALSE; |
93 |
out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size); |
94 |
Paso_SystemMatrix_allocBuffer(out->transport_matrix); |
95 |
out->flux_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size); |
96 |
out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info); |
97 |
|
98 |
out->colorOf=NULL; |
99 |
out->main_iptr=NULL; |
100 |
out->lumped_mass_matrix=NULL; |
101 |
out->row_sum_flux_matrix=NULL; |
102 |
|
103 |
if (Paso_noError()) { |
104 |
n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix); |
105 |
|
106 |
out->colorOf=MEMALLOC(n,index_t); |
107 |
out->main_iptr=MEMALLOC(n,index_t); |
108 |
out->lumped_mass_matrix=MEMALLOC(n,double); |
109 |
out->row_sum_flux_matrix=MEMALLOC(n,double); |
110 |
out->u=MEMALLOC(n,double); |
111 |
|
112 |
if ( ! (Paso_checkPtr(out->colorOf) || Paso_checkPtr(out->main_iptr) || |
113 |
Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->row_sum_flux_matrix) || Paso_checkPtr(out->u)) ) { |
114 |
|
115 |
printf("Paso_SolverFCT_getFCTransportProblem: Revise coloring!!\n"); |
116 |
Paso_Pattern_color(pattern->mainPattern,&(out->num_colors),out->colorOf); |
117 |
|
118 |
|
119 |
/* identify the main diagonals */ |
120 |
#pragma omp parallel for schedule(static) private(i,iptr,iptr_main,k) |
121 |
for (i = 0; i < n; ++i) { |
122 |
for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; ++iptr) { |
123 |
iptr_main=pattern->mainPattern->ptr[0]-1; |
124 |
for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) { |
125 |
if (pattern->mainPattern->index[iptr]==i) { |
126 |
iptr_main=iptr; |
127 |
break; |
128 |
} |
129 |
} |
130 |
out->main_iptr[i]=iptr_main; |
131 |
if (iptr_main==pattern->mainPattern->ptr[0]-1) |
132 |
Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal"); |
133 |
} |
134 |
} |
135 |
|
136 |
} |
137 |
|
138 |
} |
139 |
if (Paso_noError()) { |
140 |
return out; |
141 |
} else { |
142 |
Paso_FCTransportProblem_free(out); |
143 |
return NULL; |
144 |
} |
145 |
} |
146 |
|
147 |
void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) { |
148 |
dim_t i, n; |
149 |
|
150 |
n=Paso_FCTransportProblem_getTotalNumRows(in); |
151 |
#pragma omp parallel for schedule(static) private(i) |
152 |
for (i = 0; i < n; ++i) { |
153 |
in->u[i]=u[i]; |
154 |
} |
155 |
} |