38 |
#include "Solver.h" |
#include "Solver.h" |
39 |
#include "SolverFCT.h" |
#include "SolverFCT.h" |
40 |
#include "escript/blocktimer.h" |
#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) { |
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 |
|
printf("FFFF %e",fctp->dt_max); |
93 |
if (dt<=0.) { |
if (dt<=0.) { |
94 |
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive."); |
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt must be positive."); |
95 |
} |
} |
97 |
|
|
98 |
/* extract the row sum of the advective part */ |
/* extract the row sum of the advective part */ |
99 |
Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix); |
Paso_SystemMatrix_rowSum(fctp->flux_matrix,fctp->row_sum_flux_matrix); |
100 |
|
/* add the advective part + artificial diffusion to the diffusive part = transport-matrix */ |
|
/* add the advective part + artificial diffusion to the diffusive part */ |
|
101 |
Paso_FCTransportProblem_addAdvectivePart(fctp,1.); |
Paso_FCTransportProblem_addAdvectivePart(fctp,1.); |
102 |
|
/* create a copy of the main diagonal entires of the transport-matrix */ |
103 |
|
#pragma omp parallel for schedule(static) private(i) |
104 |
|
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
105 |
|
fctp->transport_matrix_diagonal[i]= |
106 |
|
fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]; |
107 |
|
} |
108 |
|
|
109 |
|
Paso_SystemMatrix_saveMM(fctp->flux_matrix,"flux.mm"); |
110 |
|
Paso_SystemMatrix_saveMM(fctp->transport_matrix,"trans.mm"); |
111 |
|
/*=================================================================== * |
112 |
|
* |
113 |
|
* calculate time step size: |
114 |
|
*/ |
115 |
|
dt2=fctp->dt_max; |
116 |
|
if (fctp->theta < 1.) { |
117 |
|
dt2=LARGE_POSITIVE_FLOAT; |
118 |
|
#pragma omp parallel private(dt2_loc) |
119 |
|
{ |
120 |
|
dt2_loc=LARGE_POSITIVE_FLOAT; |
121 |
|
#pragma omp for schedule(static) private(i,rtmp,rtmp2) |
122 |
|
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
123 |
|
rtmp=fctp->transport_matrix_diagonal[i]; |
124 |
|
rtmp2=fctp->lumped_mass_matrix[i]; |
125 |
|
if ( (rtmp<0 && rtmp2>=0.) || (rtmp>0 && rtmp2<=0.) ) { |
126 |
|
dt2_loc=MIN(dt2_loc,-rtmp2/rtmp); |
127 |
|
} |
128 |
|
} |
129 |
|
#pragma omp critcal |
130 |
|
{ |
131 |
|
dt2=MIN(dt2,dt2_loc); |
132 |
|
} |
133 |
|
} |
134 |
|
#ifdef PASO_MPI |
135 |
|
dt2_loc = dt2; |
136 |
|
MPI_Allreduce(&dt2_loc, &dt2, 1, MPI_DOUBLE, MPI_MIN, fctp->mpi_info->comm); |
137 |
|
#endif |
138 |
|
dt2*=1./(1.-fctp->theta); |
139 |
|
if (fctp->dt_max>0.) dt2=MIN(dt2,fctp->dt_max); |
140 |
|
} |
141 |
|
if (dt2 > 0.) { |
142 |
|
fctp->dt=dt2; |
143 |
|
} else { |
144 |
|
fctp->dt=fctp->dt_max; |
145 |
|
} |
146 |
|
if (fctp->dt < 0.) { |
147 |
|
Paso_setError(TYPE_ERROR,"Paso_SolverFCT_solve: dt_max must be positive."); |
148 |
|
} else { |
149 |
|
printf("minimum time step size is %e (theta = %e).\n",fctp->dt,fctp->theta); |
150 |
|
} |
151 |
|
/* =========================================================== |
152 |
|
* |
153 |
|
* |
154 |
|
*/ |
155 |
|
if (Paso_noError()) { |
156 |
|
|
157 |
if (Paso_noError()) fctp->valid_matrices=TRUE; |
} |
158 |
|
fctp->valid_matrices=Paso_noError(); |
159 |
} |
} |
160 |
|
|
161 |
} |
/* b_last=M*u_last + (1-theta) * F(u_last) */ |
162 |
|
|
163 |
|
/* decide on substepping */ |
164 |
|
n_substeps=ceil(dt/fctp->dt); |
165 |
|
dt2=dt/n_substeps; |
166 |
|
printf("%d time steps of size is %e (theta = %e, dt_max=%e).\n",n_substeps, dt2,fctp->theta,fctp->dt); |
167 |
|
|
168 |
|
/* */ |
169 |
|
Paso_SystemMatrix_allocBuffer(fctp->transport_matrix); |
170 |
|
|
171 |
|
if (fctp->theta>0) { |
172 |
|
#pragma omp parallel for schedule(static) private(i) |
173 |
|
for (i=0;i<Paso_SystemMatrix_getTotalNumRows(fctp->flux_matrix);++i) { |
174 |
|
fctp->transport_matrix_diagonal[i]= |
175 |
|
fctp->transport_matrix->mainBlock->val[fctp->main_iptr[i]]; |
176 |
|
} |
177 |
|
} else { |
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 |
|
} |