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