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
7
8 Initial conditions:
9 u(x,0) = 1 + sin(2 pi x)
10 v(x,0) = 3
11
12 Boundary conditions:
13 u(0,t) = u(1,t) = 1
14 v(0,t) = v(1,t) = 3
15 */
16
17 #include <petscdm.h>
18 #include <petscdmda.h>
19 #include <petscts.h>
20
21 typedef struct {
22 PetscScalar u, v;
23 } Field;
24
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 };
32
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 *);
37
main(int argc,char ** argv)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;
48
49 PetscFunctionBeginUser;
50 PetscCall(PetscInitialize(&argc, &argv, NULL, 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));
57
58 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59 Extract global vectors from DMDA;
60 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
61 PetscCall(DMCreateGlobalVector(da, &X));
62
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();
82
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));
95
96 ftime = 10.0;
97 PetscCall(TSSetMaxTime(ts, ftime));
98 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
99
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));
109
110 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
111 Set runtime options
112 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
113 PetscCall(TSSetFromOptions(ts));
114
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));
123
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 }
134
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ptr)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;
144
145 PetscFunctionBeginUser;
146 PetscCall(TSGetDM(ts, &da));
147 PetscCall(DMDAGetLocalInfo(da, &info));
148 hx = 1.0 / (PetscReal)(info.mx - 1);
149
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));
159
160 /* Get pointers to vector data */
161 PetscCall(DMDAVecGetArrayRead(da, Xloc, &x));
162 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
163 PetscCall(DMDAVecGetArray(da, F, &f));
164
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 }
178
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 }
186
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)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;
195
196 PetscFunctionBeginUser;
197 PetscCall(TSGetDM(ts, &da));
198 PetscCall(DMDAGetLocalInfo(da, &info));
199 hx = 1.0 / (PetscReal)(info.mx - 1);
200
201 /* Get pointers to vector data */
202 PetscCall(DMDAVecGetArrayRead(da, X, &x));
203 PetscCall(DMDAVecGetArray(da, F, &f));
204
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 }
211
212 /* Restore vectors */
213 PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
214 PetscCall(DMDAVecRestoreArray(da, F, &f));
215 PetscFunctionReturn(PETSC_SUCCESS);
216 }
217
218 /* --------------------------------------------------------------------- */
219 /*
220 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
221 */
FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void * ptr)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;
230
231 PetscFunctionBeginUser;
232 PetscCall(TSGetDM(ts, &da));
233 PetscCall(DMDAGetLocalInfo(da, &info));
234 hx = 1.0 / (PetscReal)(info.mx - 1);
235
236 /* Get pointers to vector data */
237 PetscCall(DMDAVecGetArrayRead(da, X, &x));
238 PetscCall(DMDAVecGetArrayRead(da, Xdot, &xdot));
239
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 }
259
260 /* Restore vectors */
261 PetscCall(DMDAVecRestoreArrayRead(da, X, &x));
262 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, &xdot));
263
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 }
272
FormInitialSolution(TS ts,Vec X,PetscCtx ctx)273 PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
274 {
275 User user = (User)ctx;
276 DM da;
277 PetscInt i;
278 DMDALocalInfo info;
279 Field *x;
280 PetscReal hx;
281
282 PetscFunctionBeginUser;
283 PetscCall(TSGetDM(ts, &da));
284 PetscCall(DMDAGetLocalInfo(da, &info));
285 hx = 1.0 / (PetscReal)(info.mx - 1);
286
287 /* Get pointers to vector data */
288 PetscCall(DMDAVecGetArray(da, X, &x));
289
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 }
299
300 /*TEST
301
302 test:
303 args: -nox -da_grid_x 20 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type 2p -ts_time_step 5e-2 -ts_adapt_type none
304
305 TEST*/
306