1 static char help[] = "Time-Dependent Reactive Flow example in 2D with Darcy Flow"; 2 3 /*T 4 Concepts: Solving a multicomponent time-dependent reactive flow system 5 Concepts: DMDA with timestepping 6 Processors: n 7 T*/ 8 9 /* 10 11 This example solves the elementary chemical reaction: 12 13 SP_A + 2SP_B = SP_C 14 15 Subject to predetermined flow modeled as though it were in a porous media. 16 This flow is modeled as advection and diffusion of the three species as 17 18 v = porosity*saturation*grad(q) 19 20 and the time-dependent equation solved for each species as 21 advection + diffusion + reaction: 22 23 v dot grad SP_i + dSP_i / dt + dispersivity*div*grad(SP_i) + R(SP_i) = 0 24 25 The following options are available: 26 27 -length_x - The length of the domain in the direction of the flow 28 -length_y - The length of the domain in the direction orthogonal to the flow 29 30 -gradq_inflow - The inflow speed as if the saturation and porosity were 1. 31 -saturation - The saturation of the porous medium. 32 -porosity - The porosity of the medium. 33 -dispersivity - The dispersivity of the flow. 34 -rate_constant - The rate constant for the chemical reaction 35 -stoich - The stoichiometry matrix for the reaction 36 -sp_inflow - The species concentrations at the inflow 37 -sp_0 - The species concentrations at time 0 38 39 */ 40 41 #include <petscdm.h> 42 #include <petscdmda.h> 43 #include <petscsnes.h> 44 #include <petscts.h> 45 46 #define N_SPECIES 3 47 #define N_REACTIONS 1 48 #define DIM 2 49 50 #define stoich(i, j) ctx->stoichiometry[N_SPECIES*i + j] 51 52 typedef struct { 53 PetscScalar sp[N_SPECIES]; 54 } Field; 55 56 typedef struct { 57 58 Field x_inflow; 59 Field x_0; 60 PetscReal stoichiometry[N_SPECIES*N_REACTIONS]; 61 PetscReal porosity; 62 PetscReal dispersivity; 63 PetscReal saturation; 64 PetscReal rate_constant[N_REACTIONS]; 65 PetscReal gradq_inflow; 66 PetscReal length[DIM]; 67 } AppCtx; 68 69 extern PetscErrorCode FormInitialGuess(DM da,AppCtx *ctx,Vec X); 70 extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo*,PetscReal,Field**,Field**,Field**,AppCtx*); 71 extern PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*); 72 extern PetscErrorCode ReactingFlowPostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 73 74 PetscErrorCode SetFromOptions(AppCtx * ctx) 75 { 76 PetscInt i,j; 77 78 PetscFunctionBeginUser; 79 ctx->dispersivity = 0.5; 80 ctx->porosity = 0.25; 81 ctx->saturation = 1.0; 82 ctx->gradq_inflow = 1.0; 83 ctx->rate_constant[0] = 0.5; 84 85 ctx->length[0] = 100.0; 86 ctx->length[1] = 100.0; 87 88 ctx->x_0.sp[0] = 0.0; 89 ctx->x_0.sp[1] = 0.0; 90 ctx->x_0.sp[2] = 0.0; 91 92 ctx->x_inflow.sp[0] = 0.05; 93 ctx->x_inflow.sp[1] = 0.05; 94 ctx->x_inflow.sp[2] = 0.0; 95 96 for (i = 0; i < N_REACTIONS; i++) { 97 for (j = 0; j < N_SPECIES; j++) stoich(i, j) = 0.; 98 } 99 100 /* set up a pretty easy example */ 101 stoich(0, 0) = -1.; 102 stoich(0, 1) = -2.; 103 stoich(0, 2) = 1.; 104 105 PetscInt as = N_SPECIES; 106 PetscCall(PetscOptionsGetReal(NULL,NULL,"-length_x",&ctx->length[0],NULL)); 107 PetscCall(PetscOptionsGetReal(NULL,NULL,"-length_y",&ctx->length[1],NULL)); 108 PetscCall(PetscOptionsGetReal(NULL,NULL,"-porosity",&ctx->porosity,NULL)); 109 PetscCall(PetscOptionsGetReal(NULL,NULL,"-saturation",&ctx->saturation,NULL)); 110 PetscCall(PetscOptionsGetReal(NULL,NULL,"-dispersivity",&ctx->dispersivity,NULL)); 111 PetscCall(PetscOptionsGetReal(NULL,NULL,"-gradq_inflow",&ctx->gradq_inflow,NULL)); 112 PetscCall(PetscOptionsGetReal(NULL,NULL,"-rate_constant",&ctx->rate_constant[0],NULL)); 113 PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-sp_inflow",ctx->x_inflow.sp,&as,NULL)); 114 PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-sp_0",ctx->x_0.sp,&as,NULL)); 115 as = N_SPECIES; 116 PetscCall(PetscOptionsGetRealArray(NULL,NULL,"-stoich",ctx->stoichiometry,&as,NULL)); 117 PetscFunctionReturn(0); 118 } 119 120 int main(int argc,char **argv) 121 { 122 TS ts; 123 SNES snes; 124 SNESLineSearch linesearch; 125 Vec x; 126 AppCtx ctx; 127 DM da; 128 129 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 130 PetscCall(SetFromOptions(&ctx)); 131 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 132 PetscCall(TSSetType(ts,TSCN)); 133 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 134 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &ctx)); 135 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-4,-4,PETSC_DECIDE,PETSC_DECIDE,N_SPECIES,1,NULL,NULL,&da)); 136 PetscCall(DMSetFromOptions(da)); 137 PetscCall(DMSetUp(da)); 138 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 139 PetscCall(DMDASetFieldName(da,0,"species A")); 140 PetscCall(DMDASetFieldName(da,1,"species B")); 141 PetscCall(DMDASetFieldName(da,2,"species C")); 142 PetscCall(DMSetApplicationContext(da,&ctx)); 143 PetscCall(DMCreateGlobalVector(da,&x)); 144 PetscCall(FormInitialGuess(da, &ctx, x)); 145 146 PetscCall(TSSetDM(ts, da)); 147 PetscCall(TSSetMaxTime(ts,1000.0)); 148 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 149 PetscCall(TSSetTimeStep(ts,1.0)); 150 PetscCall(TSSetSolution(ts,x)); 151 PetscCall(TSSetFromOptions(ts)); 152 153 PetscCall(TSGetSNES(ts,&snes)); 154 PetscCall(SNESGetLineSearch(snes,&linesearch)); 155 PetscCall(SNESLineSearchSetPostCheck(linesearch, ReactingFlowPostCheck, (void*)&ctx)); 156 PetscCall(SNESSetFromOptions(snes)); 157 PetscCall(TSSolve(ts,x)); 158 159 PetscCall(VecDestroy(&x)); 160 PetscCall(TSDestroy(&ts)); 161 PetscCall(DMDestroy(&da)); 162 PetscCall(PetscFinalize()); 163 PetscFunctionReturn(0); 164 } 165 166 /* ------------------------------------------------------------------- */ 167 168 PetscErrorCode FormInitialGuess(DM da,AppCtx *ctx,Vec X) 169 { 170 PetscInt i,j,l,Mx,My,xs,ys,xm,ym; 171 PetscErrorCode ierr; 172 Field **x; 173 174 PetscFunctionBeginUser; 175 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 176 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE); 177 178 PetscCall(DMDAVecGetArray(da,X,&x)); 179 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 180 181 for (j=ys; j<ys+ym; j++) { 182 for (i=xs; i<xs+xm; i++) { 183 for (l = 0; l < N_SPECIES; l++) { 184 if (i == 0) { 185 if (l == 0) x[j][i].sp[l] = (ctx->x_inflow.sp[l]*((PetscScalar)j) / (My - 1)); 186 else if (l == 1) x[j][i].sp[l] = (ctx->x_inflow.sp[l]*(1. - ((PetscScalar)j) / (My - 1))); 187 else x[j][i].sp[l] = ctx->x_0.sp[l]; 188 } 189 } 190 } 191 } 192 PetscCall(DMDAVecRestoreArray(da,X,&x)); 193 PetscFunctionReturn(0); 194 195 } 196 197 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscScalar ptime,Field **x,Field **xt,Field **f,AppCtx *ctx) 198 { 199 PetscInt i,j,l,m; 200 PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,scale; 201 PetscScalar u,uxx,uyy; 202 PetscScalar vx, vy,sxp,syp,sxm,sym,avx,vxp,vxm,avy,vyp,vym,f_advect; 203 PetscScalar rate; 204 205 PetscFunctionBeginUser; 206 hx = ctx->length[0]/((PetscReal)(info->mx-1)); 207 hy = ctx->length[1]/((PetscReal)(info->my-1)); 208 209 dhx = 1. / hx; 210 dhy = 1. / hy; 211 hxdhy = hx/hy; 212 hydhx = hy/hx; 213 scale = hx*hy; 214 215 for (j=info->ys; j<info->ys+info->ym; j++) { 216 for (i=info->xs; i<info->xs+info->xm; i++) { 217 vx = ctx->gradq_inflow*ctx->porosity*ctx->saturation; 218 vy = 0.0*dhy; 219 avx = PetscAbsScalar(vx); 220 vxp = .5*(vx+avx); vxm = .5*(vx-avx); 221 avy = PetscAbsScalar(vy); 222 vyp = .5*(vy+avy); vym = .5*(vy-avy); 223 /* chemical reactions */ 224 for (l = 0; l < N_SPECIES; l++) { 225 /* determine the velocites as the gradients of the pressure */ 226 if (i == 0) { 227 sxp = (x[j][i+1].sp[l] - x[j][i].sp[l])*dhx; 228 sxm = sxp; 229 } else if (i == info->mx - 1) { 230 sxp = (x[j][i].sp[l] - x[j][i-1].sp[l])*dhx; 231 sxm = sxp; 232 } else { 233 sxm = (x[j][i+1].sp[l] - x[j][i].sp[l])*dhx; 234 sxp = (x[j][i].sp[l] - x[j][i-1].sp[l])*dhx; 235 } 236 if (j == 0) { 237 syp = (x[j+1][i].sp[l] - x[j][i].sp[l])*dhy; 238 sym = syp; 239 } else if (j == info->my - 1) { 240 syp = (x[j][i].sp[l] - x[j-1][i].sp[l])*dhy; 241 sym = syp; 242 } else { 243 sym = (x[j+1][i].sp[l] - x[j][i].sp[l])*dhy; 244 syp = (x[j][i].sp[l] - x[j-1][i].sp[l])*dhy; 245 } /* 4 flops per species*point */ 246 247 if (i == 0) { 248 if (l == 0) f[j][i].sp[l] = (x[j][i].sp[l] - ctx->x_inflow.sp[l]*((PetscScalar)j) / (info->my - 1)); 249 else if (l == 1) f[j][i].sp[l] = (x[j][i].sp[l] - ctx->x_inflow.sp[l]*(1. - ((PetscScalar)j) / (info->my - 1))); 250 else f[j][i].sp[l] = x[j][i].sp[l]; 251 252 } else { 253 f[j][i].sp[l] = xt[j][i].sp[l]*scale; 254 u = x[j][i].sp[l]; 255 if (j == 0) uyy = u - x[j+1][i].sp[l]; 256 else if (j == info->my - 1) uyy = u - x[j-1][i].sp[l]; 257 else uyy = (2.0*u - x[j-1][i].sp[l] - x[j+1][i].sp[l])*hxdhy; 258 259 if (i != info->mx - 1) uxx = (2.0*u - x[j][i-1].sp[l] - x[j][i+1].sp[l])*hydhx; 260 else uxx = u - x[j][i-1].sp[l]; 261 /* 10 flops per species*point */ 262 263 f_advect = 0.; 264 f_advect += scale*(vxp*sxp + vxm*sxm); 265 f_advect += scale*(vyp*syp + vym*sym); 266 f[j][i].sp[l] += f_advect + ctx->dispersivity*(uxx + uyy); 267 /* 14 flops per species*point */ 268 } 269 } 270 /* reaction */ 271 if (i != 0) { 272 for (m = 0; m < N_REACTIONS; m++) { 273 rate = ctx->rate_constant[m]; 274 for (l = 0; l < N_SPECIES; l++) { 275 if (stoich(m, l) < 0) { 276 /* assume an elementary reaction */ 277 rate *= PetscPowScalar(x[j][i].sp[l], PetscAbsScalar(stoich(m, l))); 278 /* ~10 flops per reaction per species per point */ 279 } 280 } 281 for (l = 0; l < N_SPECIES; l++) { 282 f[j][i].sp[l] += -scale*stoich(m, l)*rate; /* Reaction term */ 283 /* 3 flops per reaction per species per point */ 284 } 285 } 286 } 287 } 288 } 289 PetscCall(PetscLogFlops((N_SPECIES*(28.0 + 13.0*N_REACTIONS))*info->ym*info->xm)); 290 PetscFunctionReturn(0); 291 } 292 293 PetscErrorCode ReactingFlowPostCheck(SNESLineSearch linesearch, Vec X, Vec Y, Vec W, PetscBool *changed_y, PetscBool *changed_w, void *vctx) 294 { 295 PetscInt i,j,l,Mx,My,xs,ys,xm,ym; 296 PetscErrorCode ierr; 297 Field **x; 298 SNES snes; 299 DM da; 300 PetscScalar min; 301 302 PetscFunctionBeginUser; 303 *changed_w = PETSC_FALSE; 304 PetscCall(VecMin(X,NULL,&min)); 305 if (min >= 0.) PetscFunctionReturn(0); 306 307 *changed_w = PETSC_TRUE; 308 PetscCall(SNESLineSearchGetSNES(linesearch, &snes)); 309 PetscCall(SNESGetDM(snes,&da)); 310 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 311 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE); 312 PetscCall(DMDAVecGetArray(da,W,&x)); 313 PetscCall(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 314 for (j=ys; j<ys+ym; j++) { 315 for (i=xs; i<xs+xm; i++) { 316 for (l = 0; l < N_SPECIES; l++) { 317 if (x[j][i].sp[l] < 0.) x[j][i].sp[l] = 0.; 318 } 319 } 320 } 321 PetscCall(DMDAVecRestoreArray(da,W,&x)); 322 PetscFunctionReturn(0); 323 } 324 325 PetscErrorCode FormIFunction(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *user) 326 { 327 DMDALocalInfo info; 328 Field **u,**udot,**fu; 329 Vec localX; 330 DM da; 331 332 PetscFunctionBeginUser; 333 PetscCall(TSGetDM(ts,(DM*)&da)); 334 PetscCall(DMGetLocalVector(da,&localX)); 335 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 336 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 337 PetscCall(DMDAGetLocalInfo(da,&info)); 338 PetscCall(DMDAVecGetArrayRead(da,localX,&u)); 339 PetscCall(DMDAVecGetArrayRead(da,Xdot,&udot)); 340 PetscCall(DMDAVecGetArray(da,F,&fu)); 341 PetscCall(FormIFunctionLocal(&info,ptime,u,udot,fu,(AppCtx*)user)); 342 PetscCall(DMDAVecRestoreArrayRead(da,localX,&u)); 343 PetscCall(DMDAVecRestoreArrayRead(da,Xdot,&udot)); 344 PetscCall(DMDAVecRestoreArray(da,F,&fu)); 345 PetscCall(DMRestoreLocalVector(da,&localX)); 346 PetscFunctionReturn(0); 347 } 348