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 57 typedef struct { 58 59 Field x_inflow; 60 Field x_0; 61 PetscReal stoichiometry[N_SPECIES*N_REACTIONS]; 62 PetscReal porosity; 63 PetscReal dispersivity; 64 PetscReal saturation; 65 PetscReal rate_constant[N_REACTIONS]; 66 PetscReal gradq_inflow; 67 PetscReal length[DIM]; 68 } AppCtx; 69 70 extern PetscErrorCode FormInitialGuess(DM da,AppCtx *ctx,Vec X); 71 extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo*,PetscReal,Field**,Field**,Field**,AppCtx*); 72 extern PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*); 73 extern PetscErrorCode ReactingFlowPostCheck(SNESLineSearch,Vec,Vec,Vec,PetscBool*,PetscBool*,void*); 74 75 76 PetscErrorCode SetFromOptions(AppCtx * ctx) 77 { 78 PetscErrorCode ierr; 79 PetscInt i,j; 80 81 PetscFunctionBeginUser; 82 ctx->dispersivity = 0.5; 83 ctx->porosity = 0.25; 84 ctx->saturation = 1.0; 85 ctx->gradq_inflow = 1.0; 86 ctx->rate_constant[0] = 0.5; 87 88 ctx->length[0] = 100.0; 89 ctx->length[1] = 100.0; 90 91 ctx->x_0.sp[0] = 0.0; 92 ctx->x_0.sp[1] = 0.0; 93 ctx->x_0.sp[2] = 0.0; 94 95 ctx->x_inflow.sp[0] = 0.05; 96 ctx->x_inflow.sp[1] = 0.05; 97 ctx->x_inflow.sp[2] = 0.0; 98 99 for (i = 0; i < N_REACTIONS; i++) { 100 for (j = 0; j < N_SPECIES; j++) stoich(i, j) = 0.; 101 } 102 103 /* set up a pretty easy example */ 104 stoich(0, 0) = -1.; 105 stoich(0, 1) = -2.; 106 stoich(0, 2) = 1.; 107 108 PetscInt as = N_SPECIES; 109 ierr = PetscOptionsGetReal(NULL,NULL,"-length_x",&ctx->length[0],NULL);CHKERRQ(ierr); 110 ierr = PetscOptionsGetReal(NULL,NULL,"-length_y",&ctx->length[1],NULL);CHKERRQ(ierr); 111 ierr = PetscOptionsGetReal(NULL,NULL,"-porosity",&ctx->porosity,NULL);CHKERRQ(ierr); 112 ierr = PetscOptionsGetReal(NULL,NULL,"-saturation",&ctx->saturation,NULL);CHKERRQ(ierr); 113 ierr = PetscOptionsGetReal(NULL,NULL,"-dispersivity",&ctx->dispersivity,NULL);CHKERRQ(ierr); 114 ierr = PetscOptionsGetReal(NULL,NULL,"-gradq_inflow",&ctx->gradq_inflow,NULL);CHKERRQ(ierr); 115 ierr = PetscOptionsGetReal(NULL,NULL,"-rate_constant",&ctx->rate_constant[0],NULL);CHKERRQ(ierr); 116 ierr = PetscOptionsGetRealArray(NULL,NULL,"-sp_inflow",ctx->x_inflow.sp,&as,NULL);CHKERRQ(ierr); 117 ierr = PetscOptionsGetRealArray(NULL,NULL,"-sp_0",ctx->x_0.sp,&as,NULL);CHKERRQ(ierr); 118 as = N_SPECIES; 119 ierr = PetscOptionsGetRealArray(NULL,NULL,"-stoich",ctx->stoichiometry,&as,NULL);CHKERRQ(ierr); 120 PetscFunctionReturn(0); 121 } 122 123 int main(int argc,char **argv) 124 { 125 TS ts; 126 SNES snes; 127 SNESLineSearch linesearch; 128 Vec x; 129 AppCtx ctx; 130 PetscErrorCode ierr; 131 DM da; 132 133 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 134 ierr = SetFromOptions(&ctx);CHKERRQ(ierr); 135 ierr = TSCreate(PETSC_COMM_WORLD, &ts);CHKERRQ(ierr); 136 ierr = TSSetType(ts,TSCN);CHKERRQ(ierr); 137 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 138 ierr = TSSetIFunction(ts, NULL, FormIFunction, &ctx);CHKERRQ(ierr); 139 ierr = 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);CHKERRQ(ierr); 140 ierr = DMSetFromOptions(da);CHKERRQ(ierr); 141 ierr = DMSetUp(da);CHKERRQ(ierr); 142 ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr); 143 ierr = DMDASetFieldName(da,0,"species A");CHKERRQ(ierr); 144 ierr = DMDASetFieldName(da,1,"species B");CHKERRQ(ierr); 145 ierr = DMDASetFieldName(da,2,"species C");CHKERRQ(ierr); 146 ierr = DMSetApplicationContext(da,&ctx);CHKERRQ(ierr); 147 ierr = DMCreateGlobalVector(da,&x);CHKERRQ(ierr); 148 ierr = FormInitialGuess(da, &ctx, x);CHKERRQ(ierr); 149 150 ierr = TSSetDM(ts, da);CHKERRQ(ierr); 151 ierr = TSSetMaxTime(ts,1000.0);CHKERRQ(ierr); 152 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 153 ierr = TSSetTimeStep(ts,1.0);CHKERRQ(ierr); 154 ierr = TSSetSolution(ts,x);CHKERRQ(ierr); 155 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 156 157 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr); 158 ierr = SNESGetLineSearch(snes,&linesearch);CHKERRQ(ierr); 159 ierr = SNESLineSearchSetPostCheck(linesearch, ReactingFlowPostCheck, (void*)&ctx);CHKERRQ(ierr); 160 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 161 ierr = TSSolve(ts,x);CHKERRQ(ierr); 162 163 ierr = VecDestroy(&x);CHKERRQ(ierr); 164 ierr = TSDestroy(&ts);CHKERRQ(ierr); 165 ierr = DMDestroy(&da);CHKERRQ(ierr); 166 ierr = PetscFinalize(); 167 PetscFunctionReturn(0); 168 } 169 170 /* ------------------------------------------------------------------- */ 171 172 PetscErrorCode FormInitialGuess(DM da,AppCtx *ctx,Vec X) 173 { 174 PetscInt i,j,l,Mx,My,xs,ys,xm,ym; 175 PetscErrorCode ierr; 176 Field **x; 177 178 PetscFunctionBeginUser; 179 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 180 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE); 181 182 ierr = DMDAVecGetArray(da,X,&x);CHKERRQ(ierr); 183 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 184 185 for (j=ys; j<ys+ym; j++) { 186 for (i=xs; i<xs+xm; i++) { 187 for (l = 0; l < N_SPECIES; l++) { 188 if (i == 0) { 189 if (l == 0) x[j][i].sp[l] = (ctx->x_inflow.sp[l]*((PetscScalar)j) / (My - 1)); 190 else if (l == 1) x[j][i].sp[l] = (ctx->x_inflow.sp[l]*(1. - ((PetscScalar)j) / (My - 1))); 191 else x[j][i].sp[l] = ctx->x_0.sp[l]; 192 } 193 } 194 } 195 } 196 ierr = DMDAVecRestoreArray(da,X,&x);CHKERRQ(ierr); 197 PetscFunctionReturn(0); 198 199 } 200 201 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info,PetscScalar ptime,Field **x,Field **xt,Field **f,AppCtx *ctx) 202 { 203 PetscErrorCode ierr; 204 PetscInt i,j,l,m; 205 PetscReal hx,hy,dhx,dhy,hxdhy,hydhx,scale; 206 PetscScalar u,uxx,uyy; 207 PetscScalar vx, vy,sxp,syp,sxm,sym,avx,vxp,vxm,avy,vyp,vym,f_advect; 208 PetscScalar rate; 209 210 PetscFunctionBeginUser; 211 hx = ctx->length[0]/((PetscReal)(info->mx-1)); 212 hy = ctx->length[1]/((PetscReal)(info->my-1)); 213 214 dhx = 1. / hx; 215 dhy = 1. / hy; 216 hxdhy = hx/hy; 217 hydhx = hy/hx; 218 scale = hx*hy; 219 220 for (j=info->ys; j<info->ys+info->ym; j++) { 221 for (i=info->xs; i<info->xs+info->xm; i++) { 222 vx = ctx->gradq_inflow*ctx->porosity*ctx->saturation; 223 vy = 0.0*dhy; 224 avx = PetscAbsScalar(vx); 225 vxp = .5*(vx+avx); vxm = .5*(vx-avx); 226 avy = PetscAbsScalar(vy); 227 vyp = .5*(vy+avy); vym = .5*(vy-avy); 228 /* chemical reactions */ 229 for (l = 0; l < N_SPECIES; l++) { 230 /* determine the velocites as the gradients of the pressure */ 231 if (i == 0) { 232 sxp = (x[j][i+1].sp[l] - x[j][i].sp[l])*dhx; 233 sxm = sxp; 234 } else if (i == info->mx - 1) { 235 sxp = (x[j][i].sp[l] - x[j][i-1].sp[l])*dhx; 236 sxm = sxp; 237 } else { 238 sxm = (x[j][i+1].sp[l] - x[j][i].sp[l])*dhx; 239 sxp = (x[j][i].sp[l] - x[j][i-1].sp[l])*dhx; 240 } 241 if (j == 0) { 242 syp = (x[j+1][i].sp[l] - x[j][i].sp[l])*dhy; 243 sym = syp; 244 } else if (j == info->my - 1) { 245 syp = (x[j][i].sp[l] - x[j-1][i].sp[l])*dhy; 246 sym = syp; 247 } else { 248 sym = (x[j+1][i].sp[l] - x[j][i].sp[l])*dhy; 249 syp = (x[j][i].sp[l] - x[j-1][i].sp[l])*dhy; 250 } /* 4 flops per species*point */ 251 252 if (i == 0) { 253 if (l == 0) f[j][i].sp[l] = (x[j][i].sp[l] - ctx->x_inflow.sp[l]*((PetscScalar)j) / (info->my - 1)); 254 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))); 255 else f[j][i].sp[l] = x[j][i].sp[l]; 256 257 } else { 258 f[j][i].sp[l] = xt[j][i].sp[l]*scale; 259 u = x[j][i].sp[l]; 260 if (j == 0) uyy = u - x[j+1][i].sp[l]; 261 else if (j == info->my - 1) uyy = u - x[j-1][i].sp[l]; 262 else uyy = (2.0*u - x[j-1][i].sp[l] - x[j+1][i].sp[l])*hxdhy; 263 264 if (i != info->mx - 1) uxx = (2.0*u - x[j][i-1].sp[l] - x[j][i+1].sp[l])*hydhx; 265 else uxx = u - x[j][i-1].sp[l]; 266 /* 10 flops per species*point */ 267 268 f_advect = 0.; 269 f_advect += scale*(vxp*sxp + vxm*sxm); 270 f_advect += scale*(vyp*syp + vym*sym); 271 f[j][i].sp[l] += f_advect + ctx->dispersivity*(uxx + uyy); 272 /* 14 flops per species*point */ 273 } 274 } 275 /* reaction */ 276 if (i != 0) { 277 for (m = 0; m < N_REACTIONS; m++) { 278 rate = ctx->rate_constant[m]; 279 for (l = 0; l < N_SPECIES; l++) { 280 if (stoich(m, l) < 0) { 281 /* assume an elementary reaction */ 282 rate *= PetscPowScalar(x[j][i].sp[l], PetscAbsScalar(stoich(m, l))); 283 /* ~10 flops per reaction per species per point */ 284 } 285 } 286 for (l = 0; l < N_SPECIES; l++) { 287 f[j][i].sp[l] += -scale*stoich(m, l)*rate; /* Reaction term */ 288 /* 3 flops per reaction per species per point */ 289 } 290 } 291 } 292 } 293 } 294 ierr = PetscLogFlops((N_SPECIES*(28 + 13*N_REACTIONS))*info->ym*info->xm);CHKERRQ(ierr); 295 PetscFunctionReturn(0); 296 } 297 298 299 PetscErrorCode ReactingFlowPostCheck(SNESLineSearch linesearch, Vec X, Vec Y, Vec W, PetscBool *changed_y, PetscBool *changed_w, void *vctx) 300 { 301 PetscInt i,j,l,Mx,My,xs,ys,xm,ym; 302 PetscErrorCode ierr; 303 Field **x; 304 SNES snes; 305 DM da; 306 PetscScalar min; 307 308 PetscFunctionBeginUser; 309 *changed_w = PETSC_FALSE; 310 ierr = VecMin(X,NULL,&min);CHKERRQ(ierr); 311 if (min >= 0.) PetscFunctionReturn(0); 312 313 *changed_w = PETSC_TRUE; 314 ierr = SNESLineSearchGetSNES(linesearch, &snes);CHKERRQ(ierr); 315 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 316 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE, 317 PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE); 318 ierr = DMDAVecGetArray(da,W,&x);CHKERRQ(ierr); 319 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr); 320 for (j=ys; j<ys+ym; j++) { 321 for (i=xs; i<xs+xm; i++) { 322 for (l = 0; l < N_SPECIES; l++) { 323 if (x[j][i].sp[l] < 0.) x[j][i].sp[l] = 0.; 324 } 325 } 326 } 327 ierr = DMDAVecRestoreArray(da,W,&x);CHKERRQ(ierr); 328 PetscFunctionReturn(0); 329 } 330 331 332 PetscErrorCode FormIFunction(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void *user) 333 { 334 DMDALocalInfo info; 335 Field **u,**udot,**fu; 336 PetscErrorCode ierr; 337 Vec localX; 338 DM da; 339 340 PetscFunctionBeginUser; 341 ierr = TSGetDM(ts,(DM*)&da);CHKERRQ(ierr); 342 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 343 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 344 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 345 ierr = DMDAGetLocalInfo(da,&info);CHKERRQ(ierr); 346 ierr = DMDAVecGetArrayRead(da,localX,&u);CHKERRQ(ierr); 347 ierr = DMDAVecGetArrayRead(da,Xdot,&udot);CHKERRQ(ierr); 348 ierr = DMDAVecGetArray(da,F,&fu);CHKERRQ(ierr); 349 ierr = FormIFunctionLocal(&info,ptime,u,udot,fu,(AppCtx*)user);CHKERRQ(ierr); 350 ierr = DMDAVecRestoreArrayRead(da,localX,&u);CHKERRQ(ierr); 351 ierr = DMDAVecRestoreArrayRead(da,Xdot,&udot);CHKERRQ(ierr); 352 ierr = DMDAVecRestoreArray(da,F,&fu);CHKERRQ(ierr); 353 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 354 PetscFunctionReturn(0); 355 } 356