Actual source code: ex25.c
1: static const char help[] = "Call PetscInitialize multiple times.\n";
2: /*
3: This example is based on the Brusselator tutorial of the same name, but tests multiple calls to PetscInitialize().
4: This is a bad "convergence study" because it only compares min and max values of the solution rather than comparing
5: norms of the errors. For convergence studies, we recommend invoking PetscInitialize() only once and comparing norms
6: of errors (perhaps estimated using an accurate reference solution).
8: Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and multiple solves.
10: u_t - alpha u_xx = A + u^2 v - (B+1) u
11: v_t - alpha v_xx = B u - u^2 v
12: 0 < x < 1;
13: A = 1, B = 3, alpha = 1/10
15: Initial conditions:
16: u(x,0) = 1 + sin(2 pi x)
17: v(x,0) = 3
19: Boundary conditions:
20: u(0,t) = u(1,t) = 1
21: v(0,t) = v(1,t) = 3
22: */
24: #include <petscdm.h>
25: #include <petscdmda.h>
26: #include <petscts.h>
28: typedef struct {
29: PetscScalar u, v;
30: } Field;
32: typedef struct _User *User;
33: struct _User {
34: PetscReal A, B; /* Reaction coefficients */
35: PetscReal alpha; /* Diffusion coefficient */
36: PetscReal uleft, uright; /* Dirichlet boundary conditions */
37: PetscReal vleft, vright; /* Dirichlet boundary conditions */
38: };
40: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
41: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
42: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
43: static PetscErrorCode FormInitialSolution(TS, Vec, void *);
44: static int Brusselator(int, char **, PetscInt);
46: int main(int argc, char **argv)
47: {
48: PetscInt cycle;
51: MPI_Init(&argc, &argv);
52: if (ierr) return ierr;
53: for (cycle = 0; cycle < 4; cycle++) {
54: Brusselator(argc, argv, cycle);
55: if (ierr) return 1;
56: }
57: MPI_Finalize();
58: return ierr;
59: }
61: PetscErrorCode Brusselator(int argc, char **argv, PetscInt cycle)
62: {
63: TS ts; /* nonlinear solver */
64: Vec X; /* solution, residual vectors */
65: Mat J; /* Jacobian matrix */
66: PetscInt steps, mx;
67: DM da;
68: PetscReal ftime, hx, dt, xmax, xmin;
69: struct _User user; /* user-defined work context */
70: TSConvergedReason reason;
73: PetscInitialize(&argc, &argv, (char *)0, help);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create distributed array (DMDA) to manage parallel grid and vectors
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da);
79: DMSetFromOptions(da);
80: DMSetUp(da);
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Extract global vectors from DMDA;
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85: DMCreateGlobalVector(da, &X);
87: /* Initialize user application context */
88: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
89: {
90: user.A = 1;
91: user.B = 3;
92: user.alpha = 0.1;
93: user.uleft = 1;
94: user.uright = 1;
95: user.vleft = 3;
96: user.vright = 3;
97: PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL);
98: PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL);
99: PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL);
100: PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL);
101: PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL);
102: PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL);
103: PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL);
104: }
105: PetscOptionsEnd();
107: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108: Create timestepping solver context
109: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110: TSCreate(PETSC_COMM_WORLD, &ts);
111: TSSetDM(ts, da);
112: TSSetType(ts, TSARKIMEX);
113: TSSetRHSFunction(ts, NULL, FormRHSFunction, &user);
114: TSSetIFunction(ts, NULL, FormIFunction, &user);
115: DMSetMatType(da, MATAIJ);
116: DMCreateMatrix(da, &J);
117: TSSetIJacobian(ts, J, J, FormIJacobian, &user);
119: ftime = 1.0;
120: TSSetMaxTime(ts, ftime);
121: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Set initial conditions
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: FormInitialSolution(ts, X, &user);
127: TSSetSolution(ts, X);
128: VecGetSize(X, &mx);
129: hx = 1.0 / (PetscReal)(mx / 2 - 1);
130: dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
131: dt *= PetscPowRealInt(0.2, cycle); /* Shrink the time step in convergence study. */
132: TSSetTimeStep(ts, dt);
133: TSSetTolerances(ts, 1e-3 * PetscPowRealInt(0.5, cycle), NULL, 1e-3 * PetscPowRealInt(0.5, cycle), NULL);
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Set runtime options
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: TSSetFromOptions(ts);
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Solve nonlinear system
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: TSSolve(ts, X);
144: TSGetSolveTime(ts, &ftime);
145: TSGetStepNumber(ts, &steps);
146: TSGetConvergedReason(ts, &reason);
147: VecMin(X, NULL, &xmin);
148: VecMax(X, NULL, &xmax);
149: PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after % 3" PetscInt_FMT " steps. Range [%6.4f,%6.4f]\n", TSConvergedReasons[reason], (double)ftime, steps, (double)xmin, (double)xmax);
151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: Free work space.
153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154: MatDestroy(&J);
155: VecDestroy(&X);
156: TSDestroy(&ts);
157: DMDestroy(&da);
158: PetscFinalize();
159: return 0;
160: }
162: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
163: {
164: User user = (User)ptr;
165: DM da;
166: DMDALocalInfo info;
167: PetscInt i;
168: Field *x, *xdot, *f;
169: PetscReal hx;
170: Vec Xloc;
173: TSGetDM(ts, &da);
174: DMDAGetLocalInfo(da, &info);
175: hx = 1.0 / (PetscReal)(info.mx - 1);
177: /*
178: Scatter ghost points to local vector,using the 2-step process
179: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
180: By placing code between these two statements, computations can be
181: done while messages are in transition.
182: */
183: DMGetLocalVector(da, &Xloc);
184: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
185: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
187: /* Get pointers to vector data */
188: DMDAVecGetArrayRead(da, Xloc, &x);
189: DMDAVecGetArrayRead(da, Xdot, &xdot);
190: DMDAVecGetArray(da, F, &f);
192: /* Compute function over the locally owned part of the grid */
193: for (i = info.xs; i < info.xs + info.xm; i++) {
194: if (i == 0) {
195: f[i].u = hx * (x[i].u - user->uleft);
196: f[i].v = hx * (x[i].v - user->vleft);
197: } else if (i == info.mx - 1) {
198: f[i].u = hx * (x[i].u - user->uright);
199: f[i].v = hx * (x[i].v - user->vright);
200: } else {
201: f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
202: f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
203: }
204: }
206: /* Restore vectors */
207: DMDAVecRestoreArrayRead(da, Xloc, &x);
208: DMDAVecRestoreArrayRead(da, Xdot, &xdot);
209: DMDAVecRestoreArray(da, F, &f);
210: DMRestoreLocalVector(da, &Xloc);
211: return 0;
212: }
214: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
215: {
216: User user = (User)ptr;
217: DM da;
218: DMDALocalInfo info;
219: PetscInt i;
220: PetscReal hx;
221: Field *x, *f;
224: TSGetDM(ts, &da);
225: DMDAGetLocalInfo(da, &info);
226: hx = 1.0 / (PetscReal)(info.mx - 1);
228: /* Get pointers to vector data */
229: DMDAVecGetArrayRead(da, X, &x);
230: DMDAVecGetArray(da, F, &f);
232: /* Compute function over the locally owned part of the grid */
233: for (i = info.xs; i < info.xs + info.xm; i++) {
234: PetscScalar u = x[i].u, v = x[i].v;
235: f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u);
236: f[i].v = hx * (user->B * u - u * u * v);
237: }
239: /* Restore vectors */
240: DMDAVecRestoreArrayRead(da, X, &x);
241: DMDAVecRestoreArray(da, F, &f);
242: return 0;
243: }
245: /* --------------------------------------------------------------------- */
246: /*
247: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
248: */
249: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
250: {
251: User user = (User)ptr;
252: DMDALocalInfo info;
253: PetscInt i;
254: PetscReal hx;
255: DM da;
256: Field *x, *xdot;
259: TSGetDM(ts, &da);
260: DMDAGetLocalInfo(da, &info);
261: hx = 1.0 / (PetscReal)(info.mx - 1);
263: /* Get pointers to vector data */
264: DMDAVecGetArrayRead(da, X, &x);
265: DMDAVecGetArrayRead(da, Xdot, &xdot);
267: /* Compute function over the locally owned part of the grid */
268: for (i = info.xs; i < info.xs + info.xm; i++) {
269: if (i == 0 || i == info.mx - 1) {
270: const PetscInt row = i, col = i;
271: const PetscScalar vals[2][2] = {
272: {hx, 0 },
273: {0, hx}
274: };
275: MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES);
276: } else {
277: const PetscInt row = i, col[] = {i - 1, i, i + 1};
278: const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
279: const PetscScalar vals[2][3][2] = {
280: {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
281: {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
282: };
283: MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
284: }
285: }
287: /* Restore vectors */
288: DMDAVecRestoreArrayRead(da, X, &x);
289: DMDAVecRestoreArrayRead(da, Xdot, &xdot);
291: MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
292: MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
293: if (J != Jpre) {
294: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
295: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
296: }
297: return 0;
298: }
300: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
301: {
302: User user = (User)ctx;
303: DM da;
304: PetscInt i;
305: DMDALocalInfo info;
306: Field *x;
307: PetscReal hx;
310: TSGetDM(ts, &da);
311: DMDAGetLocalInfo(da, &info);
312: hx = 1.0 / (PetscReal)(info.mx - 1);
314: /* Get pointers to vector data */
315: DMDAVecGetArray(da, X, &x);
317: /* Compute function over the locally owned part of the grid */
318: for (i = info.xs; i < info.xs + info.xm; i++) {
319: PetscReal xi = i * hx;
320: x[i].u = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi);
321: x[i].v = user->vleft * (1. - xi) + user->vright * xi;
322: }
323: DMDAVecRestoreArray(da, X, &x);
324: return 0;
325: }
327: /*TEST
329: test:
330: args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3
332: test:
333: suffix: 2
334: args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3
336: TEST*/