Actual source code: ex27.c
1: static char help[] = "Time-Dependent Reactive Flow example in 2D with Darcy Flow";
3: /*
5: This example solves the elementary chemical reaction:
7: SP_A + 2SP_B = SP_C
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
12: v = porosity*saturation*grad(q)
14: and the time-dependent equation solved for each species as
15: advection + diffusion + reaction:
17: v dot grad SP_i + dSP_i / dt + dispersivity*div*grad(SP_i) + R(SP_i) = 0
19: The following options are available:
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
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
33: */
35: #include <petscdm.h>
36: #include <petscdmda.h>
37: #include <petscsnes.h>
38: #include <petscts.h>
40: #define N_SPECIES 3
41: #define N_REACTIONS 1
42: #define DIM 2
44: #define stoich(i, j) ctx->stoichiometry[N_SPECIES * i + j]
46: typedef struct {
47: PetscScalar sp[N_SPECIES];
48: } Field;
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;
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 *);
67: PetscErrorCode SetFromOptions(AppCtx *ctx)
68: {
69: PetscInt i, j;
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;
78: ctx->length[0] = 100.0;
79: ctx->length[1] = 100.0;
81: ctx->x_0.sp[0] = 0.0;
82: ctx->x_0.sp[1] = 0.0;
83: ctx->x_0.sp[2] = 0.0;
85: ctx->x_inflow.sp[0] = 0.05;
86: ctx->x_inflow.sp[1] = 0.05;
87: ctx->x_inflow.sp[2] = 0.0;
89: for (i = 0; i < N_REACTIONS; i++) {
90: for (j = 0; j < N_SPECIES; j++) stoich(i, j) = 0.;
91: }
93: /* set up a pretty easy example */
94: stoich(0, 0) = -1.;
95: stoich(0, 1) = -2.;
96: stoich(0, 2) = 1.;
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: }
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;
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));
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));
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));
153: PetscCall(VecDestroy(&x));
154: PetscCall(TSDestroy(&ts));
155: PetscCall(DMDestroy(&da));
156: PetscCall(PetscFinalize());
157: PetscFunctionReturn(PETSC_SUCCESS);
158: }
160: /* ------------------------------------------------------------------- */
162: PetscErrorCode FormInitialGuess(DM da, AppCtx *ctx, Vec X)
163: {
164: PetscInt i, j, l, Mx, My, xs, ys, xm, ym;
165: Field **x;
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));
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: }
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;
195: PetscFunctionBeginUser;
196: hx = ctx->length[0] / ((PetscReal)(info->mx - 1));
197: hy = ctx->length[1] / ((PetscReal)(info->my - 1));
199: dhx = 1. / hx;
200: dhy = 1. / hy;
201: hxdhy = hx / hy;
202: hydhx = hy / hx;
203: scale = hx * hy;
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 */
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];
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;
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 */
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: }
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;
293: PetscFunctionBeginUser;
294: *changed_w = PETSC_FALSE;
295: PetscCall(VecMin(X, NULL, &min));
296: if (min >= 0.) PetscFunctionReturn(PETSC_SUCCESS);
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: }
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;
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: }