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 Field x_inflow; 52 Field x_0; 53 PetscReal stoichiometry[N_SPECIES * N_REACTIONS]; 54 PetscReal porosity; 55 PetscReal dispersivity; 56 PetscReal saturation; 57 PetscReal rate_constant[N_REACTIONS]; 58 PetscReal gradq_inflow; 59 PetscReal length[DIM]; 60 } AppCtx; 61 62 extern PetscErrorCode FormInitialGuess(DM da, AppCtx *ctx, Vec X); 63 extern PetscErrorCode FormIFunctionLocal(DMDALocalInfo *, PetscReal, Field **, Field **, Field **, AppCtx *); 64 extern PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *); 65 extern PetscErrorCode ReactingFlowPostCheck(SNESLineSearch, Vec, Vec, Vec, PetscBool *, PetscBool *, void *); 66 67 PetscErrorCode SetFromOptions(AppCtx *ctx) 68 { 69 PetscInt i, j; 70 71 PetscFunctionBeginUser; 72 ctx->dispersivity = 0.5; 73 ctx->porosity = 0.25; 74 ctx->saturation = 1.0; 75 ctx->gradq_inflow = 1.0; 76 ctx->rate_constant[0] = 0.5; 77 78 ctx->length[0] = 100.0; 79 ctx->length[1] = 100.0; 80 81 ctx->x_0.sp[0] = 0.0; 82 ctx->x_0.sp[1] = 0.0; 83 ctx->x_0.sp[2] = 0.0; 84 85 ctx->x_inflow.sp[0] = 0.05; 86 ctx->x_inflow.sp[1] = 0.05; 87 ctx->x_inflow.sp[2] = 0.0; 88 89 for (i = 0; i < N_REACTIONS; i++) { 90 for (j = 0; j < N_SPECIES; j++) stoich(i, j) = 0.; 91 } 92 93 /* set up a pretty easy example */ 94 stoich(0, 0) = -1.; 95 stoich(0, 1) = -2.; 96 stoich(0, 2) = 1.; 97 98 PetscInt as = N_SPECIES; 99 PetscCall(PetscOptionsGetReal(NULL, NULL, "-length_x", &ctx->length[0], NULL)); 100 PetscCall(PetscOptionsGetReal(NULL, NULL, "-length_y", &ctx->length[1], NULL)); 101 PetscCall(PetscOptionsGetReal(NULL, NULL, "-porosity", &ctx->porosity, NULL)); 102 PetscCall(PetscOptionsGetReal(NULL, NULL, "-saturation", &ctx->saturation, NULL)); 103 PetscCall(PetscOptionsGetReal(NULL, NULL, "-dispersivity", &ctx->dispersivity, NULL)); 104 PetscCall(PetscOptionsGetReal(NULL, NULL, "-gradq_inflow", &ctx->gradq_inflow, NULL)); 105 PetscCall(PetscOptionsGetReal(NULL, NULL, "-rate_constant", &ctx->rate_constant[0], NULL)); 106 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-sp_inflow", ctx->x_inflow.sp, &as, NULL)); 107 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-sp_0", ctx->x_0.sp, &as, NULL)); 108 as = N_SPECIES; 109 PetscCall(PetscOptionsGetRealArray(NULL, NULL, "-stoich", ctx->stoichiometry, &as, NULL)); 110 PetscFunctionReturn(PETSC_SUCCESS); 111 } 112 113 int main(int argc, char **argv) 114 { 115 TS ts; 116 SNES snes; 117 SNESLineSearch linesearch; 118 Vec x; 119 AppCtx ctx; 120 DM da; 121 122 PetscFunctionBeginUser; 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(PETSC_SUCCESS); 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, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 169 PetscCall(DMDAVecGetArray(da, X, &x)); 170 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 171 172 for (j = ys; j < ys + ym; j++) { 173 for (i = xs; i < xs + xm; i++) { 174 for (l = 0; l < N_SPECIES; l++) { 175 if (i == 0) { 176 if (l == 0) x[j][i].sp[l] = (ctx->x_inflow.sp[l] * ((PetscScalar)j) / (My - 1)); 177 else if (l == 1) x[j][i].sp[l] = (ctx->x_inflow.sp[l] * (1. - ((PetscScalar)j) / (My - 1))); 178 else x[j][i].sp[l] = ctx->x_0.sp[l]; 179 } 180 } 181 } 182 } 183 PetscCall(DMDAVecRestoreArray(da, X, &x)); 184 PetscFunctionReturn(PETSC_SUCCESS); 185 } 186 187 PetscErrorCode FormIFunctionLocal(DMDALocalInfo *info, PetscScalar ptime, Field **x, Field **xt, Field **f, AppCtx *ctx) 188 { 189 PetscInt i, j, l, m; 190 PetscReal hx, hy, dhx, dhy, hxdhy, hydhx, scale; 191 PetscScalar u, uxx, uyy; 192 PetscScalar vx, vy, sxp, syp, sxm, sym, avx, vxp, vxm, avy, vyp, vym, f_advect; 193 PetscScalar rate; 194 195 PetscFunctionBeginUser; 196 hx = ctx->length[0] / ((PetscReal)(info->mx - 1)); 197 hy = ctx->length[1] / ((PetscReal)(info->my - 1)); 198 199 dhx = 1. / hx; 200 dhy = 1. / hy; 201 hxdhy = hx / hy; 202 hydhx = hy / hx; 203 scale = hx * hy; 204 205 for (j = info->ys; j < info->ys + info->ym; j++) { 206 for (i = info->xs; i < info->xs + info->xm; i++) { 207 vx = ctx->gradq_inflow * ctx->porosity * ctx->saturation; 208 vy = 0.0 * dhy; 209 avx = PetscAbsScalar(vx); 210 vxp = .5 * (vx + avx); 211 vxm = .5 * (vx - avx); 212 avy = PetscAbsScalar(vy); 213 vyp = .5 * (vy + avy); 214 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 302 PetscCall(DMDAVecGetArray(da, W, &x)); 303 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 304 for (j = ys; j < ys + ym; j++) { 305 for (i = xs; i < xs + xm; i++) { 306 for (l = 0; l < N_SPECIES; l++) { 307 if (x[j][i].sp[l] < 0.) x[j][i].sp[l] = 0.; 308 } 309 } 310 } 311 PetscCall(DMDAVecRestoreArray(da, W, &x)); 312 PetscFunctionReturn(PETSC_SUCCESS); 313 } 314 315 PetscErrorCode FormIFunction(TS ts, PetscReal ptime, Vec X, Vec Xdot, Vec F, void *user) 316 { 317 DMDALocalInfo info; 318 Field **u, **udot, **fu; 319 Vec localX; 320 DM da; 321 322 PetscFunctionBeginUser; 323 PetscCall(TSGetDM(ts, (DM *)&da)); 324 PetscCall(DMGetLocalVector(da, &localX)); 325 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX)); 326 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX)); 327 PetscCall(DMDAGetLocalInfo(da, &info)); 328 PetscCall(DMDAVecGetArrayRead(da, localX, &u)); 329 PetscCall(DMDAVecGetArrayRead(da, Xdot, &udot)); 330 PetscCall(DMDAVecGetArray(da, F, &fu)); 331 PetscCall(FormIFunctionLocal(&info, ptime, u, udot, fu, (AppCtx *)user)); 332 PetscCall(DMDAVecRestoreArrayRead(da, localX, &u)); 333 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &udot)); 334 PetscCall(DMDAVecRestoreArray(da, F, &fu)); 335 PetscCall(DMRestoreLocalVector(da, &localX)); 336 PetscFunctionReturn(PETSC_SUCCESS); 337 } 338