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 PetscErrorCode Brusselator(int, char **, PetscInt);
46: int main(int argc, char **argv)
47: {
48: int ierr = MPI_Init(&argc, &argv);
49: if (ierr) return ierr;
50: for (PetscInt cycle = 0; cycle < 4; cycle++) {
51: if (Brusselator(argc, argv, cycle)) return 1;
52: }
53: return MPI_Finalize();
54: }
56: PetscErrorCode Brusselator(int argc, char **argv, PetscInt cycle)
57: {
58: TS ts; /* nonlinear solver */
59: Vec X; /* solution, residual vectors */
60: Mat J; /* Jacobian matrix */
61: PetscInt steps, mx;
62: DM da;
63: PetscReal ftime, hx, dt, xmax, xmin;
64: struct _User user; /* user-defined work context */
65: TSConvergedReason reason;
67: PetscFunctionBeginUser;
68: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
70: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71: Create distributed array (DMDA) to manage parallel grid and vectors
72: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da));
74: PetscCall(DMSetFromOptions(da));
75: PetscCall(DMSetUp(da));
77: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78: Extract global vectors from DMDA;
79: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80: PetscCall(DMCreateGlobalVector(da, &X));
82: /* Initialize user application context */
83: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
84: {
85: user.A = 1;
86: user.B = 3;
87: user.alpha = 0.1;
88: user.uleft = 1;
89: user.uright = 1;
90: user.vleft = 3;
91: user.vright = 3;
92: PetscCall(PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL));
93: PetscCall(PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL));
94: PetscCall(PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL));
95: PetscCall(PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL));
96: PetscCall(PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL));
97: PetscCall(PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL));
98: PetscCall(PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL));
99: }
100: PetscOptionsEnd();
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create timestepping solver context
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
106: PetscCall(TSSetDM(ts, da));
107: PetscCall(TSSetType(ts, TSARKIMEX));
108: PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
109: PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
110: PetscCall(DMSetMatType(da, MATAIJ));
111: PetscCall(DMCreateMatrix(da, &J));
112: PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));
114: ftime = 1.0;
115: PetscCall(TSSetMaxTime(ts, ftime));
116: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
118: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: Set initial conditions
120: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: PetscCall(FormInitialSolution(ts, X, &user));
122: PetscCall(TSSetSolution(ts, X));
123: PetscCall(VecGetSize(X, &mx));
124: hx = 1.0 / (PetscReal)(mx / 2 - 1);
125: dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
126: dt *= PetscPowRealInt(0.2, cycle); /* Shrink the time step in convergence study. */
127: PetscCall(TSSetTimeStep(ts, dt));
128: PetscCall(TSSetTolerances(ts, 1e-3 * PetscPowRealInt(0.5, cycle), NULL, 1e-3 * PetscPowRealInt(0.5, cycle), NULL));
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Set runtime options
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: PetscCall(TSSetFromOptions(ts));
135: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136: Solve nonlinear system
137: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscCall(TSSolve(ts, X));
139: PetscCall(TSGetSolveTime(ts, &ftime));
140: PetscCall(TSGetStepNumber(ts, &steps));
141: PetscCall(TSGetConvergedReason(ts, &reason));
142: PetscCall(VecMin(X, NULL, &xmin));
143: PetscCall(VecMax(X, NULL, &xmax));
144: PetscCall(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));
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Free work space.
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: PetscCall(MatDestroy(&J));
150: PetscCall(VecDestroy(&X));
151: PetscCall(TSDestroy(&ts));
152: PetscCall(DMDestroy(&da));
153: PetscCall(PetscFinalize());
154: return PETSC_SUCCESS;
155: }
157: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
158: {
159: User user = (User)ptr;
160: DM da;
161: DMDALocalInfo info;
162: PetscInt i;
163: Field *x, *xdot, *f;
164: PetscReal hx;
165: Vec Xloc;
167: PetscFunctionBeginUser;
168: PetscCall(TSGetDM(ts, &da));
169: PetscCall(DMDAGetLocalInfo(da, &info));
170: hx = 1.0 / (PetscReal)(info.mx - 1);
172: /*
173: Scatter ghost points to local vector,using the 2-step process
174: DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
175: By placing code between these two statements, computations can be
176: done while messages are in transition.
177: */
178: PetscCall(DMGetLocalVector(da, &Xloc));
179: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
180: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
182: /* Get pointers to vector data */
183: PetscCall(DMDAVecGetArrayRead(da, Xloc, &x));
184: PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
185: PetscCall(DMDAVecGetArray(da, F, &f));
187: /* Compute function over the locally owned part of the grid */
188: for (i = info.xs; i < info.xs + info.xm; i++) {
189: if (i == 0) {
190: f[i].u = hx * (x[i].u - user->uleft);
191: f[i].v = hx * (x[i].v - user->vleft);
192: } else if (i == info.mx - 1) {
193: f[i].u = hx * (x[i].u - user->uright);
194: f[i].v = hx * (x[i].v - user->vright);
195: } else {
196: f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
197: f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
198: }
199: }
201: /* Restore vectors */
202: PetscCall(DMDAVecRestoreArrayRead(da, Xloc, &x));
203: PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
204: PetscCall(DMDAVecRestoreArray(da, F, &f));
205: PetscCall(DMRestoreLocalVector(da, &Xloc));
206: PetscFunctionReturn(PETSC_SUCCESS);
207: }
209: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
210: {
211: User user = (User)ptr;
212: DM da;
213: DMDALocalInfo info;
214: PetscInt i;
215: PetscReal hx;
216: Field *x, *f;
218: PetscFunctionBeginUser;
219: PetscCall(TSGetDM(ts, &da));
220: PetscCall(DMDAGetLocalInfo(da, &info));
221: hx = 1.0 / (PetscReal)(info.mx - 1);
223: /* Get pointers to vector data */
224: PetscCall(DMDAVecGetArrayRead(da, X, &x));
225: PetscCall(DMDAVecGetArray(da, F, &f));
227: /* Compute function over the locally owned part of the grid */
228: for (i = info.xs; i < info.xs + info.xm; i++) {
229: PetscScalar u = x[i].u, v = x[i].v;
230: f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u);
231: f[i].v = hx * (user->B * u - u * u * v);
232: }
234: /* Restore vectors */
235: PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
236: PetscCall(DMDAVecRestoreArray(da, F, &f));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: /* --------------------------------------------------------------------- */
241: /*
242: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
243: */
244: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
245: {
246: User user = (User)ptr;
247: DMDALocalInfo info;
248: PetscInt i;
249: PetscReal hx;
250: DM da;
251: Field *x, *xdot;
253: PetscFunctionBeginUser;
254: PetscCall(TSGetDM(ts, &da));
255: PetscCall(DMDAGetLocalInfo(da, &info));
256: hx = 1.0 / (PetscReal)(info.mx - 1);
258: /* Get pointers to vector data */
259: PetscCall(DMDAVecGetArrayRead(da, X, &x));
260: PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
262: /* Compute function over the locally owned part of the grid */
263: for (i = info.xs; i < info.xs + info.xm; i++) {
264: if (i == 0 || i == info.mx - 1) {
265: const PetscInt row = i, col = i;
266: const PetscScalar vals[2][2] = {
267: {hx, 0 },
268: {0, hx}
269: };
270: PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES));
271: } else {
272: const PetscInt row = i, col[] = {i - 1, i, i + 1};
273: const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
274: const PetscScalar vals[2][3][2] = {
275: {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
276: {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
277: };
278: PetscCall(MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES));
279: }
280: }
282: /* Restore vectors */
283: PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
284: PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
286: PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
287: PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
288: if (J != Jpre) {
289: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
290: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
291: }
292: PetscFunctionReturn(PETSC_SUCCESS);
293: }
295: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
296: {
297: User user = (User)ctx;
298: DM da;
299: PetscInt i;
300: DMDALocalInfo info;
301: Field *x;
302: PetscReal hx;
304: PetscFunctionBeginUser;
305: PetscCall(TSGetDM(ts, &da));
306: PetscCall(DMDAGetLocalInfo(da, &info));
307: hx = 1.0 / (PetscReal)(info.mx - 1);
309: /* Get pointers to vector data */
310: PetscCall(DMDAVecGetArray(da, X, &x));
312: /* Compute function over the locally owned part of the grid */
313: for (i = info.xs; i < info.xs + info.xm; i++) {
314: PetscReal xi = i * hx;
315: x[i].u = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi);
316: x[i].v = user->vleft * (1. - xi) + user->vright * xi;
317: }
318: PetscCall(DMDAVecRestoreArray(da, X, &x));
319: PetscFunctionReturn(PETSC_SUCCESS);
320: }
322: /*TEST
324: test:
325: args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3
327: test:
328: suffix: 2
329: args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3
331: TEST*/