1 |
gross |
1361 |
/* $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 |
gross |
1362 |
/* Paso: FCTransportProblem */ |
17 |
gross |
1361 |
|
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 |
gross |
1362 |
/**************************************************************/ |
31 |
gross |
1361 |
|
32 |
gross |
1362 |
/* 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 |
gross |
1361 |
|
38 |
gross |
1362 |
Paso_SystemMatrix_free(in->transport_matrix); |
39 |
gross |
1407 |
Paso_SystemMatrix_free(in->mass_matrix); |
40 |
|
|
Paso_SystemMatrix_free(in->iteration_matrix); |
41 |
gross |
1362 |
Paso_MPIInfo_free(in->mpi_info); |
42 |
gross |
1364 |
MEMFREE(in->u); |
43 |
gross |
1407 |
MEMFREE(in->main_iptr); |
44 |
gross |
1362 |
MEMFREE(in->lumped_mass_matrix); |
45 |
gross |
1407 |
MEMFREE(in->main_diagonal_low_order_transport_matrix); |
46 |
gross |
1362 |
MEMFREE(in); |
47 |
|
|
} |
48 |
|
|
} |
49 |
|
|
} |
50 |
gross |
1361 |
|
51 |
gross |
1362 |
Paso_FCTransportProblem* Paso_FCTransportProblem_getReference(Paso_FCTransportProblem* in) { |
52 |
gross |
1361 |
if (in!=NULL) { |
53 |
gross |
1362 |
++(in->reference_counter); |
54 |
gross |
1361 |
} |
55 |
gross |
1367 |
return in; |
56 |
gross |
1362 |
} |
57 |
|
|
|
58 |
|
|
Paso_SystemMatrix* Paso_FCTransportProblem_borrowTransportMatrix(Paso_FCTransportProblem* in) { |
59 |
|
|
return in->transport_matrix; |
60 |
gross |
1361 |
} |
61 |
gross |
1407 |
Paso_SystemMatrix* Paso_FCTransportProblem_borrowMassMatrix(Paso_FCTransportProblem* in) { |
62 |
|
|
return in->mass_matrix; |
63 |
gross |
1362 |
} |
64 |
gross |
1361 |
|
65 |
gross |
1362 |
double* Paso_FCTransportProblem_borrowLumpedMassMatrix(Paso_FCTransportProblem* in) { |
66 |
|
|
return in->lumped_mass_matrix; |
67 |
|
|
} |
68 |
gross |
1361 |
|
69 |
gross |
1362 |
dim_t Paso_FCTransportProblem_getTotalNumRows(Paso_FCTransportProblem* in) { |
70 |
|
|
return Paso_SystemMatrix_getTotalNumRows(in->transport_matrix); |
71 |
|
|
} |
72 |
gross |
1361 |
|
73 |
gross |
1407 |
Paso_FCTransportProblem* Paso_FCTransportProblem_alloc(double theta, Paso_SystemMatrixPattern *pattern, int block_size |
74 |
gross |
1361 |
|
75 |
gross |
1362 |
|
76 |
|
|
) { |
77 |
|
|
Paso_SystemMatrixType matrix_type=MATRIX_FORMAT_DEFAULT+MATRIX_FORMAT_BLK1; /* at the moment only block size 1 is supported */ |
78 |
|
|
Paso_FCTransportProblem* out=NULL; |
79 |
|
|
dim_t n,i; |
80 |
trankine |
1462 |
index_t iptr,iptr_main; |
81 |
gross |
1362 |
|
82 |
|
|
if ((theta<0.) || (theta >1.)) { |
83 |
|
|
Paso_setError(TYPE_ERROR,"Paso_FCTransportProblem_alloc: theta needs to be between 0. and. 1."); |
84 |
gross |
1361 |
return NULL; |
85 |
gross |
1362 |
} |
86 |
gross |
1361 |
|
87 |
gross |
1362 |
out=MEMALLOC(1,Paso_FCTransportProblem); |
88 |
|
|
if (Paso_checkPtr(out)) return NULL; |
89 |
gross |
1361 |
|
90 |
gross |
1407 |
|
91 |
gross |
1362 |
out->theta=theta; |
92 |
gross |
1407 |
out->dt_max=LARGE_POSITIVE_FLOAT; |
93 |
gross |
1363 |
out->valid_matrices=FALSE; |
94 |
gross |
1362 |
out->transport_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size); |
95 |
|
|
Paso_SystemMatrix_allocBuffer(out->transport_matrix); |
96 |
gross |
1407 |
out->mass_matrix=Paso_SystemMatrix_alloc(matrix_type,pattern,block_size,block_size); |
97 |
|
|
out->iteration_matrix=NULL; |
98 |
|
|
|
99 |
gross |
1362 |
out->mpi_info=Paso_MPIInfo_getReference(pattern->mpi_info); |
100 |
gross |
1407 |
out->u=NULL; |
101 |
|
|
out->u_min=0.; |
102 |
gross |
1362 |
out->main_iptr=NULL; |
103 |
|
|
out->lumped_mass_matrix=NULL; |
104 |
gross |
1407 |
out->main_diagonal_low_order_transport_matrix=NULL; |
105 |
gross |
1361 |
|
106 |
gross |
1362 |
if (Paso_noError()) { |
107 |
|
|
n=Paso_SystemMatrix_getTotalNumRows(out->transport_matrix); |
108 |
|
|
|
109 |
gross |
1407 |
out->u=MEMALLOC(n,double); |
110 |
gross |
1362 |
out->main_iptr=MEMALLOC(n,index_t); |
111 |
|
|
out->lumped_mass_matrix=MEMALLOC(n,double); |
112 |
gross |
1407 |
out->main_diagonal_low_order_transport_matrix=MEMALLOC(n,double); |
113 |
gross |
1362 |
|
114 |
gross |
1407 |
if ( ! (Paso_checkPtr(out->u) || Paso_checkPtr(out->main_iptr) || |
115 |
|
|
Paso_checkPtr(out->lumped_mass_matrix) || Paso_checkPtr(out->main_diagonal_low_order_transport_matrix)) ) { |
116 |
gross |
1362 |
|
117 |
gross |
1367 |
#pragma omp parallel for schedule(static) private(i) |
118 |
|
|
for (i = 0; i < n; ++i) { |
119 |
|
|
out->lumped_mass_matrix[i]=0.; |
120 |
gross |
1407 |
out->main_diagonal_low_order_transport_matrix[i]=0.; |
121 |
gross |
1367 |
out->u[i]=0.; |
122 |
|
|
} |
123 |
gross |
1362 |
/* identify the main diagonals */ |
124 |
ksteube |
1558 |
#pragma omp parallel for schedule(static) private(i,iptr,iptr_main) |
125 |
gross |
1362 |
for (i = 0; i < n; ++i) { |
126 |
|
|
iptr_main=pattern->mainPattern->ptr[0]-1; |
127 |
|
|
for (iptr=pattern->mainPattern->ptr[i];iptr<pattern->mainPattern->ptr[i+1]; iptr++) { |
128 |
|
|
if (pattern->mainPattern->index[iptr]==i) { |
129 |
|
|
iptr_main=iptr; |
130 |
|
|
break; |
131 |
|
|
} |
132 |
|
|
} |
133 |
|
|
out->main_iptr[i]=iptr_main; |
134 |
|
|
if (iptr_main==pattern->mainPattern->ptr[0]-1) |
135 |
|
|
Paso_setError(VALUE_ERROR, "Paso_FCTransportProblem_alloc: no main diagonal"); |
136 |
|
|
} |
137 |
|
|
|
138 |
|
|
} |
139 |
|
|
|
140 |
gross |
1361 |
} |
141 |
|
|
if (Paso_noError()) { |
142 |
|
|
return out; |
143 |
|
|
} else { |
144 |
gross |
1362 |
Paso_FCTransportProblem_free(out); |
145 |
gross |
1361 |
return NULL; |
146 |
|
|
} |
147 |
|
|
} |
148 |
gross |
1364 |
|
149 |
|
|
void Paso_FCTransportProblem_checkinSolution(Paso_FCTransportProblem* in, double* u) { |
150 |
|
|
dim_t i, n; |
151 |
gross |
1407 |
double local_u_min,u_min; |
152 |
gross |
1364 |
|
153 |
|
|
n=Paso_FCTransportProblem_getTotalNumRows(in); |
154 |
gross |
1407 |
u_min=LARGE_POSITIVE_FLOAT; |
155 |
|
|
#pragma omp parallel private(local_u_min) |
156 |
|
|
{ |
157 |
|
|
local_u_min=0.; |
158 |
|
|
#pragma omp for schedule(static) private(i) |
159 |
gross |
1410 |
for (i=0;i<n;++i) local_u_min=MIN(local_u_min,u[i]); |
160 |
gross |
1407 |
#pragma omp critical |
161 |
|
|
{ |
162 |
|
|
u_min=MIN(u_min,local_u_min); |
163 |
|
|
} |
164 |
gross |
1364 |
} |
165 |
gross |
1407 |
#ifdef PASO_MPI |
166 |
|
|
local_u_min=u_min; |
167 |
ksteube |
1415 |
MPI_Allreduce(&local_u_min,&u_min, 1, MPI_DOUBLE, MPI_MIN, in->mpi_info->comm); |
168 |
gross |
1407 |
#endif |
169 |
|
|
in->u_min=u_min; |
170 |
gross |
1410 |
#pragma omp parallel for schedule(static) private(i) |
171 |
|
|
for (i=0;i<n;++i) { |
172 |
|
|
in->u[i]=u[i]-u_min; |
173 |
|
|
} |
174 |
gross |
1364 |
} |