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;
50: PetscInitialize(&argc, &argv, (char *)0, help);
51: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: Create distributed array (DMDA) to manage parallel grid and vectors
53: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
54: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da);
55: DMSetFromOptions(da);
56: DMSetUp(da);
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Extract global vectors from DMDA;
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61: 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: PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL);
74: PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL);
75: PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL);
76: PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL);
77: PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL);
78: PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL);
79: PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL);
80: }
81: PetscOptionsEnd();
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create timestepping solver context
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: TSCreate(PETSC_COMM_WORLD, &ts);
87: TSSetDM(ts, da);
88: TSSetType(ts, TSARKIMEX);
89: TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1);
90: TSSetRHSFunction(ts, NULL, FormRHSFunction, &user);
91: TSSetIFunction(ts, NULL, FormIFunction, &user);
92: DMSetMatType(da, MATAIJ);
93: DMCreateMatrix(da, &J);
94: TSSetIJacobian(ts, J, J, FormIJacobian, &user);
96: ftime = 10.0;
97: TSSetMaxTime(ts, ftime);
98: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Set initial conditions
102: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
103: FormInitialSolution(ts, X, &user);
104: TSSetSolution(ts, X);
105: VecGetSize(X, &mx);
106: hx = 1.0 / (PetscReal)(mx / 2 - 1);
107: dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
108: TSSetTimeStep(ts, dt);
110: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111: Set runtime options
112: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113: TSSetFromOptions(ts);
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Solve nonlinear system
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118: TSSolve(ts, X);
119: TSGetSolveTime(ts, &ftime);
120: TSGetStepNumber(ts, &steps);
121: TSGetConvergedReason(ts, &reason);
122: 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: MatDestroy(&J);
128: VecDestroy(&X);
129: TSDestroy(&ts);
130: DMDestroy(&da);
131: 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;
146: TSGetDM(ts, &da);
147: 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: DMGetLocalVector(da, &Xloc);
157: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
158: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
160: /* Get pointers to vector data */
161: DMDAVecGetArrayRead(da, Xloc, &x);
162: DMDAVecGetArrayRead(da, Xdot, &xdot);
163: 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: DMDAVecRestoreArrayRead(da, Xloc, &x);
181: DMDAVecRestoreArrayRead(da, Xdot, &xdot);
182: DMDAVecRestoreArray(da, F, &f);
183: DMRestoreLocalVector(da, &Xloc);
184: return 0;
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;
197: TSGetDM(ts, &da);
198: DMDAGetLocalInfo(da, &info);
199: hx = 1.0 / (PetscReal)(info.mx - 1);
201: /* Get pointers to vector data */
202: DMDAVecGetArrayRead(da, X, &x);
203: 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: DMDAVecRestoreArrayRead(da, X, &x);
214: DMDAVecRestoreArray(da, F, &f);
215: return 0;
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;
232: TSGetDM(ts, &da);
233: DMDAGetLocalInfo(da, &info);
234: hx = 1.0 / (PetscReal)(info.mx - 1);
236: /* Get pointers to vector data */
237: DMDAVecGetArrayRead(da, X, &x);
238: 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: 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: MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
257: }
258: }
260: /* Restore vectors */
261: DMDAVecRestoreArrayRead(da, X, &x);
262: DMDAVecRestoreArrayRead(da, Xdot, &xdot);
264: MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
265: MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
266: if (J != Jpre) {
267: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
268: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
269: }
270: return 0;
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;
283: TSGetDM(ts, &da);
284: DMDAGetLocalInfo(da, &info);
285: hx = 1.0 / (PetscReal)(info.mx - 1);
287: /* Get pointers to vector data */
288: 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: DMDAVecRestoreArray(da, X, &x);
297: return 0;
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*/