Actual source code: ex25.c

```  1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d formulated as a PDAE. Demonstrates solving PDEs with algebraic constraints (PDAE).\n";
2: /*
3:    u_t - alpha u_xx = A + u^2 v - (B+1) u
4:    v_t - alpha v_xx = B u - u^2 v
5:    0 < x < 1;
6:    A = 1, B = 3, alpha = 1/50

8:    Initial conditions:
9:    u(x,0) = 1 + sin(2 pi x)
10:    v(x,0) = 3

12:    Boundary conditions:
13:    u(0,t) = u(1,t) = 1
14:    v(0,t) = v(1,t) = 3
15: */

17: #include <petscdm.h>
18: #include <petscdmda.h>
19: #include <petscts.h>

21: typedef struct {
22:   PetscScalar u, v;
23: } Field;

25: typedef struct _User *User;
26: struct _User {
27:   PetscReal A, B;          /* Reaction coefficients */
28:   PetscReal alpha;         /* Diffusion coefficient */
29:   PetscReal uleft, uright; /* Dirichlet boundary conditions */
30:   PetscReal vleft, vright; /* Dirichlet boundary conditions */
31: };

33: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
34: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
35: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
36: static PetscErrorCode FormInitialSolution(TS, Vec, void *);

38: int main(int argc, char **argv)
39: {
40:   TS                ts; /* nonlinear solver */
41:   Vec               X;  /* solution, residual vectors */
42:   Mat               J;  /* Jacobian matrix */
43:   PetscInt          steps, mx;
44:   DM                da;
45:   PetscReal         ftime, hx, dt;
46:   struct _User      user; /* user-defined work context */
47:   TSConvergedReason reason;

49:   PetscFunctionBeginUser;
50:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52:      Create distributed array (DMDA) to manage parallel grid and vectors
53:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54:   PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da));
55:   PetscCall(DMSetFromOptions(da));
56:   PetscCall(DMSetUp(da));

58:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59:      Extract global vectors from DMDA;
60:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61:   PetscCall(DMCreateGlobalVector(da, &X));

63:   /* Initialize user application context */
64:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
65:   {
66:     user.A      = 1;
67:     user.B      = 3;
68:     user.alpha  = 0.02;
69:     user.uleft  = 1;
70:     user.uright = 1;
71:     user.vleft  = 3;
72:     user.vright = 3;
73:     PetscCall(PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL));
74:     PetscCall(PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL));
75:     PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL));
76:     PetscCall(PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL));
77:     PetscCall(PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL));
78:     PetscCall(PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL));
79:     PetscCall(PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL));
80:   }
81:   PetscOptionsEnd();

83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84:      Create timestepping solver context
85:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
87:   PetscCall(TSSetDM(ts, da));
88:   PetscCall(TSSetType(ts, TSARKIMEX));
89:   PetscCall(TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1));
90:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
91:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
92:   PetscCall(DMSetMatType(da, MATAIJ));
93:   PetscCall(DMCreateMatrix(da, &J));
94:   PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));

96:   ftime = 10.0;
97:   PetscCall(TSSetMaxTime(ts, ftime));
98:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:      Set initial conditions
102:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103:   PetscCall(FormInitialSolution(ts, X, &user));
104:   PetscCall(TSSetSolution(ts, X));
105:   PetscCall(VecGetSize(X, &mx));
106:   hx = 1.0 / (PetscReal)(mx / 2 - 1);
107:   dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
108:   PetscCall(TSSetTimeStep(ts, dt));

110:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111:      Set runtime options
112:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113:   PetscCall(TSSetFromOptions(ts));

115:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116:      Solve nonlinear system
117:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118:   PetscCall(TSSolve(ts, X));
119:   PetscCall(TSGetSolveTime(ts, &ftime));
120:   PetscCall(TSGetStepNumber(ts, &steps));
121:   PetscCall(TSGetConvergedReason(ts, &reason));
122:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:      Free work space.
126:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
127:   PetscCall(MatDestroy(&J));
128:   PetscCall(VecDestroy(&X));
129:   PetscCall(TSDestroy(&ts));
130:   PetscCall(DMDestroy(&da));
131:   PetscCall(PetscFinalize());
132:   return 0;
133: }

135: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
136: {
137:   User          user = (User)ptr;
138:   DM            da;
139:   DMDALocalInfo info;
140:   PetscInt      i;
141:   Field        *x, *xdot, *f;
142:   PetscReal     hx;
143:   Vec           Xloc;

145:   PetscFunctionBeginUser;
146:   PetscCall(TSGetDM(ts, &da));
147:   PetscCall(DMDAGetLocalInfo(da, &info));
148:   hx = 1.0 / (PetscReal)(info.mx - 1);

150:   /*
151:      Scatter ghost points to local vector,using the 2-step process
152:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
153:      By placing code between these two statements, computations can be
154:      done while messages are in transition.
155:   */
156:   PetscCall(DMGetLocalVector(da, &Xloc));
157:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
158:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));

160:   /* Get pointers to vector data */
163:   PetscCall(DMDAVecGetArray(da, F, &f));

165:   /* Compute function over the locally owned part of the grid */
166:   for (i = info.xs; i < info.xs + info.xm; i++) {
167:     if (i == 0) {
168:       f[i].u = hx * (x[i].u - user->uleft);
169:       f[i].v = hx * (x[i].v - user->vleft);
170:     } else if (i == info.mx - 1) {
171:       f[i].u = hx * (x[i].u - user->uright);
172:       f[i].v = hx * (x[i].v - user->vright);
173:     } else {
174:       f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
175:       f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
176:     }
177:   }

179:   /* Restore vectors */
182:   PetscCall(DMDAVecRestoreArray(da, F, &f));
183:   PetscCall(DMRestoreLocalVector(da, &Xloc));
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
188: {
189:   User          user = (User)ptr;
190:   DM            da;
191:   DMDALocalInfo info;
192:   PetscInt      i;
193:   PetscReal     hx;
194:   Field        *x, *f;

196:   PetscFunctionBeginUser;
197:   PetscCall(TSGetDM(ts, &da));
198:   PetscCall(DMDAGetLocalInfo(da, &info));
199:   hx = 1.0 / (PetscReal)(info.mx - 1);

201:   /* Get pointers to vector data */
203:   PetscCall(DMDAVecGetArray(da, F, &f));

205:   /* Compute function over the locally owned part of the grid */
206:   for (i = info.xs; i < info.xs + info.xm; i++) {
207:     PetscScalar u = x[i].u, v = x[i].v;
208:     f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u);
209:     f[i].v = hx * (user->B * u - u * u * v);
210:   }

212:   /* Restore vectors */
214:   PetscCall(DMDAVecRestoreArray(da, F, &f));
215:   PetscFunctionReturn(PETSC_SUCCESS);
216: }

218: /* --------------------------------------------------------------------- */
219: /*
220:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
221: */
222: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
223: {
224:   User          user = (User)ptr;
225:   DMDALocalInfo info;
226:   PetscInt      i;
227:   PetscReal     hx;
228:   DM            da;
229:   Field        *x, *xdot;

231:   PetscFunctionBeginUser;
232:   PetscCall(TSGetDM(ts, &da));
233:   PetscCall(DMDAGetLocalInfo(da, &info));
234:   hx = 1.0 / (PetscReal)(info.mx - 1);

236:   /* Get pointers to vector data */

240:   /* Compute function over the locally owned part of the grid */
241:   for (i = info.xs; i < info.xs + info.xm; i++) {
242:     if (i == 0 || i == info.mx - 1) {
243:       const PetscInt    row = i, col = i;
244:       const PetscScalar vals[2][2] = {
245:         {hx, 0 },
246:         {0,  hx}
247:       };
248:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES));
249:     } else {
250:       const PetscInt    row = i, col[] = {i - 1, i, i + 1};
251:       const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
252:       const PetscScalar vals[2][3][2] = {
253:         {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
254:         {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
255:       };
256:       PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES));
257:     }
258:   }

260:   /* Restore vectors */

264:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
265:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
266:   if (J != Jpre) {
267:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
268:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
269:   }
270:   PetscFunctionReturn(PETSC_SUCCESS);
271: }

273: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
274: {
275:   User          user = (User)ctx;
276:   DM            da;
277:   PetscInt      i;
278:   DMDALocalInfo info;
279:   Field        *x;
280:   PetscReal     hx;

282:   PetscFunctionBeginUser;
283:   PetscCall(TSGetDM(ts, &da));
284:   PetscCall(DMDAGetLocalInfo(da, &info));
285:   hx = 1.0 / (PetscReal)(info.mx - 1);

287:   /* Get pointers to vector data */
288:   PetscCall(DMDAVecGetArray(da, X, &x));

290:   /* Compute function over the locally owned part of the grid */
291:   for (i = info.xs; i < info.xs + info.xm; i++) {
292:     PetscReal xi = i * hx;
293:     x[i].u       = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi);
294:     x[i].v       = user->vleft * (1. - xi) + user->vright * xi;
295:   }
296:   PetscCall(DMDAVecRestoreArray(da, X, &x));
297:   PetscFunctionReturn(PETSC_SUCCESS);
298: }

300: /*TEST

302:     test:
303:       args: -nox -da_grid_x 20 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none

305: TEST*/
```