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).
7
8 Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and multiple solves.
9
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
14
15 Initial conditions:
16 u(x,0) = 1 + sin(2 pi x)
17 v(x,0) = 3
18
19 Boundary conditions:
20 u(0,t) = u(1,t) = 1
21 v(0,t) = v(1,t) = 3
22 */
23
24 #include <petscdm.h>
25 #include <petscdmda.h>
26 #include <petscts.h>
27
28 typedef struct {
29 PetscScalar u, v;
30 } Field;
31
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 };
39
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);
45
main(int argc,char ** argv)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 }
55
Brusselator(int argc,char ** argv,PetscInt cycle)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;
66
67 PetscFunctionBeginUser;
68 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
69
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));
76
77 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
78 Extract global vectors from DMDA;
79 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
80 PetscCall(DMCreateGlobalVector(da, &X));
81
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();
101
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));
113
114 ftime = 1.0;
115 PetscCall(TSSetMaxTime(ts, ftime));
116 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
117
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));
129
130 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131 Set runtime options
132 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133 PetscCall(TSSetFromOptions(ts));
134
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));
145
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 }
156
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ptr)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;
166
167 PetscFunctionBeginUser;
168 PetscCall(TSGetDM(ts, &da));
169 PetscCall(DMDAGetLocalInfo(da, &info));
170 hx = 1.0 / (PetscReal)(info.mx - 1);
171
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));
181
182 /* Get pointers to vector data */
183 PetscCall(DMDAVecGetArrayRead(da, Xloc, &x));
184 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
185 PetscCall(DMDAVecGetArray(da, F, &f));
186
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 }
200
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 }
208
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
217
218 PetscFunctionBeginUser;
219 PetscCall(TSGetDM(ts, &da));
220 PetscCall(DMDAGetLocalInfo(da, &info));
221 hx = 1.0 / (PetscReal)(info.mx - 1);
222
223 /* Get pointers to vector data */
224 PetscCall(DMDAVecGetArrayRead(da, X, &x));
225 PetscCall(DMDAVecGetArray(da, F, &f));
226
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 }
233
234 /* Restore vectors */
235 PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
236 PetscCall(DMDAVecRestoreArray(da, F, &f));
237 PetscFunctionReturn(PETSC_SUCCESS);
238 }
239
240 /* --------------------------------------------------------------------- */
241 /*
242 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
243 */
FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void * ptr)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;
252
253 PetscFunctionBeginUser;
254 PetscCall(TSGetDM(ts, &da));
255 PetscCall(DMDAGetLocalInfo(da, &info));
256 hx = 1.0 / (PetscReal)(info.mx - 1);
257
258 /* Get pointers to vector data */
259 PetscCall(DMDAVecGetArrayRead(da, X, &x));
260 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
261
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 }
281
282 /* Restore vectors */
283 PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
284 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
285
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 }
294
FormInitialSolution(TS ts,Vec X,PetscCtx ctx)295 PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
296 {
297 User user = (User)ctx;
298 DM da;
299 PetscInt i;
300 DMDALocalInfo info;
301 Field *x;
302 PetscReal hx;
303
304 PetscFunctionBeginUser;
305 PetscCall(TSGetDM(ts, &da));
306 PetscCall(DMDAGetLocalInfo(da, &info));
307 hx = 1.0 / (PetscReal)(info.mx - 1);
308
309 /* Get pointers to vector data */
310 PetscCall(DMDAVecGetArray(da, X, &x));
311
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 }
321
322 /*TEST
323
324 test:
325 args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3
326
327 test:
328 suffix: 2
329 args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3
330
331 TEST*/
332