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
SetFromOptions(AppCtx * ctx)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
main(int argc,char ** argv)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, NULL, 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
FormInitialGuess(DM da,AppCtx * ctx,Vec X)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
FormIFunctionLocal(DMDALocalInfo * info,PetscScalar ptime,Field ** x,Field ** xt,Field ** f,AppCtx * ctx)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
ReactingFlowPostCheck(SNESLineSearch linesearch,Vec X,Vec Y,Vec W,PetscBool * changed_y,PetscBool * changed_w,void * vctx)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
FormIFunction(TS ts,PetscReal ptime,Vec X,Vec Xdot,Vec F,void * user)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