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: Flux correction transport solver |
17 |
* |
18 |
* solves Mu_t=Du+Ku+q |
19 |
* |
20 |
* where is D is diffusive (not checked) |
21 |
* - D is symmetric |
22 |
* - row sums are equal to zero. |
23 |
* and K is the advective part. |
24 |
* |
25 |
* u(0) >= 0 |
26 |
* |
27 |
* intially fctp->transport_matrix defines the diffusive part |
28 |
* but the matrix is updated by the adevctive part + artificial diffusion |
29 |
* |
30 |
*/ |
31 |
/**************************************************************/ |
32 |
|
33 |
/* Author: l.gross@uq.edu.au */ |
34 |
|
35 |
/**************************************************************/ |
36 |
|
37 |
#include "Paso.h" |
38 |
#include "Solver.h" |
39 |
#include "SolverFCT.h" |
40 |
#include "escript/blocktimer.h" |
41 |
#ifdef _OPENMP |
42 |
#include <omp.h> |
43 |
#endif |
44 |
#ifdef PASO_MPI |
45 |
#include <mpi.h> |
46 |
#endif |
47 |
|
48 |
/*********************************************************************************** |
49 |
|
50 |
solve (until convergence) |
51 |
|
52 |
(*) [M-theta*dt*L] du = M*u_last + (1-theta)*dt*F(u_last) - M*u + theta*dt F(u) |
53 |
u <- u+du |
54 |
|
55 |
with F(u) = L u + f_a(u) (f_a anti-diffusion) |
56 |
and L = D + K + D_K stored in transport_matrix |
57 |
D = Diffusive part (on input stored in transport_matrix) |
58 |
K = flux_matrix |
59 |
D_K = artificial diffusion introduced by K |
60 |
f_a = anti-diffusion introduced by K and u |
61 |
u_last=u of last time step |
62 |
|
63 |
For the case theta==0 (explicit case) no iteration is required. One sets |
64 |
|
65 |
M*u= M*u_last + dt*F(u_last) |
66 |
|
67 |
with F(u) = L u + f_a(u). |
68 |
|
69 |
|
70 |
For the case theta>0 we set |
71 |
|
72 |
A=L-M/(theta*dt) (stored into transport_matrix) |
73 |
F(u)=L u + f_a(u) = M/(theta*dt) u + A u + f_a(u) |
74 |
b(u)= - M/(theta*dt) * u + F(u) = A u + f_a(u) |
75 |
b_last=M/(theta*dt)*u_last + (1-theta)*dt * F(u_last) |
76 |
=M/(theta*dt)*u_last + (1-theta)/theta * [ M/(theta*dt) u_last + A u_last + f_a(u_last) ] |
77 |
=(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last) |
78 |
so (*) takes the form |
79 |
|
80 |
b_last=(1/theta**2*dt)*M u_last + (1-theta)/theta*b(u_last) |
81 |
while (|du| > tol * |u|) : |
82 |
A*du=b_last + b(u) |
83 |
u <- u-du |
84 |
|
85 |
*/ |
86 |
|
87 |
void Paso_SolverFCT_solve(Paso_FCTransportProblem* fctp, double* u, double dt, double* source, Paso_Options* options) { |
88 |
|
89 |
index_t i; |
90 |
int n_substeps,n; |
91 |
double dt2=fctp->dt_max, dt2_loc, rtmp,rtmp2,t; |
92 |
if (dt<=0.) { |
93 |
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive."); |
94 |
} |
95 |
if (! fctp->valid_matrices) { |
96 |
|
97 |
/* extract the row sum of the advective part */ |
98 |
Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix); |
99 |
/* add the advective part + artificial diffusion to the diffusive part = transport-matrix */ |
100 |
Paso_FCTransportProblem_addAdvectivePart(fctp,1.); |
101 |
/* create a copy of the main diagonal entires of the transport-matrix */ |
102 |
#pragma omp parallel for schedule(static) private(i) |
103 |
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
104 |
fctp->transport_matrix_diagonal[i]= |
105 |
fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]; |
106 |
} |
107 |
|
108 |
Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm"); |
109 |
Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm"); |
110 |
/*=================================================================== * |
111 |
* |
112 |
* calculate time step size: |
113 |
*/ |
114 |
dt2=fctp->dt_max; |
115 |
if (fctp->theta < 1.) { |
116 |
dt2=LARGE_POSITIVE_FLOAT; |
117 |
#pragma omp parallel private(dt2_loc) |
118 |
{ |
119 |
dt2_loc=LARGE_POSITIVE_FLOAT; |
120 |
#pragma omp for schedule(static) private(i,rtmp,rtmp2) |
121 |
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
122 |
rtmp=fctp->transport_matrix_diagonal[i]; |
123 |
rtmp2=fctp->lumped_mass_matrix[i]; |
124 |
if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) { |
125 |
dt2_loc=MIN(dt2_loc,-rtmp2/rtmp); |
126 |
} |
127 |
} |
128 |
#pragma omp critcal |
129 |
{ |
130 |
dt2=MIN(dt2,dt2_loc); |
131 |
} |
132 |
} |
133 |
#ifdef PASO_MPI |
134 |
dt2_loc = dt2; |
135 |
MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm); |
136 |
#endif |
137 |
dt2*=1./(1.-fctp->theta); |
138 |
if (fctp->dt_max>0.) dt2=MIN(dt2,fctp->dt_max); |
139 |
} |
140 |
if (dt2 > 0.) { |
141 |
fctp->dt=dt2; |
142 |
} else { |
143 |
fctp->dt=fctp->dt_max; |
144 |
} |
145 |
if (fctp->dt < 0.) { |
146 |
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt_max must be positive."); |
147 |
} else { |
148 |
printf("minimum time step size is %e (theta = %e).\n",fctp->dt,fctp->theta); |
149 |
} |
150 |
/* =========================================================== |
151 |
* |
152 |
* |
153 |
*/ |
154 |
if (Paso_noError()) { |
155 |
|
156 |
} |
157 |
fctp->valid_matrices=Paso_noError(); |
158 |
} |
159 |
|
160 |
/* b_last=M*u_last + (1-theta) * F(u_last) */ |
161 |
|
162 |
/* decide on substepping */ |
163 |
n_substeps=ceil(dt/fctp->dt); |
164 |
dt2=dt/n_substeps; |
165 |
printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt); |
166 |
|
167 |
/* */ |
168 |
Paso_SystemMatrix_allocBuffer(fctp->transport_matrix); |
169 |
|
170 |
if (fctp->theta>0) { |
171 |
#pragma omp parallel for schedule(static) private(i) |
172 |
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
173 |
fctp->transport_matrix_diagonal[i]= |
174 |
fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]; |
175 |
} |
176 |
} else { |
177 |
/* n_substeps=1; */ |
178 |
/* u= u_last + M^-1*dt*F(u_last) */ |
179 |
t=dt2; |
180 |
n=0; |
181 |
while(n<n_substeps) { |
182 |
printf("substep step %d at t=%e\n",n+1,t); |
183 |
Paso_FCTransportProblem_setFlux(fctp,fctp->u,u); /* u stores F(u_last)*/ |
184 |
#pragma omp parallel for schedule(static) private(i) |
185 |
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
186 |
rtmp=fctp->u[i]; |
187 |
fctp->u[i]=fctp->u[i]+dt2*u[i]/fctp->lumped_mass_matrix[i]; |
188 |
printf("%d : %e %e %e\n",i,u[i],fctp->u[i],rtmp); |
189 |
} |
190 |
t+=dt2; |
191 |
n++; |
192 |
} |
193 |
} |
194 |
if (Paso_noError()) { |
195 |
#pragma omp parallel for schedule(static) private(i) |
196 |
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
197 |
u[i]=fctp->u[i]; |
198 |
} |
199 |
} |
200 |
} |