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 */
161: PetscCall(DMDAVecGetArrayRead(da, Xloc, &x));
162: PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
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 */
180: PetscCall(DMDAVecRestoreArrayRead(da, Xloc, &x));
181: PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
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 */
202: PetscCall(DMDAVecGetArrayRead(da, X, &x));
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 */
213: PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
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 */
237: PetscCall(DMDAVecGetArrayRead(da, X, &x));
238: PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
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 */
261: PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
262: PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
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*/