| 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 |
} |
Parent Directory
| 
Revision Log