1 static const char help[] = "Time-dependent advection-reaction PDE in 1d, demonstrates IMEX methods.\n";
2 /*
3 u_t + a1*u_x = -k1*u + k2*v + s1
4 v_t + a2*v_x = k1*u - k2*v + s2
5 0 < x < 1;
6 a1 = 1, k1 = 10^6, s1 = 0,
7 a2 = 0, k2 = 2*k1, s2 = 1
8
9 Initial conditions:
10 u(x,0) = 1 + s2*x
11 v(x,0) = k0/k1*u(x,0) + s1/k1
12
13 Upstream boundary conditions:
14 u(0,t) = 1-sin(12*t)^4
15
16 Useful command-line parameters:
17 -ts_arkimex_fully_implicit - treats advection implicitly, only relevant with -ts_type arkimex (default)
18 -ts_type rosw - use Rosenbrock-W method (linearly implicit IMEX, amortizes assembly and preconditioner setup)
19 */
20
21 #include <petscdm.h>
22 #include <petscdmda.h>
23 #include <petscts.h>
24
25 typedef PetscScalar Field[2];
26
27 typedef struct _User *User;
28 struct _User {
29 PetscReal a[2]; /* Advection speeds */
30 PetscReal k[2]; /* Reaction coefficients */
31 PetscReal s[2]; /* Source terms */
32 };
33
34 static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
35 static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
36 static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);
37 static PetscErrorCode FormInitialSolution(TS, Vec, void *);
38
main(int argc,char ** argv)39 int main(int argc, char **argv)
40 {
41 TS ts; /* time integrator */
42 SNES snes; /* nonlinear solver */
43 SNESLineSearch linesearch; /* line search */
44 Vec X; /* solution, residual vectors */
45 Mat J; /* Jacobian matrix */
46 PetscInt steps, mx;
47 DM da;
48 PetscReal ftime, dt;
49 struct _User user; /* user-defined work context */
50 TSConvergedReason reason;
51
52 PetscFunctionBeginUser;
53 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
54
55 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
56 Create distributed array (DMDA) to manage parallel grid and vectors
57 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
58 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da));
59 PetscCall(DMSetFromOptions(da));
60 PetscCall(DMSetUp(da));
61
62 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63 Extract global vectors from DMDA;
64 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65 PetscCall(DMCreateGlobalVector(da, &X));
66
67 /* Initialize user application context */
68 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
69 {
70 user.a[0] = 1;
71 PetscCall(PetscOptionsReal("-a0", "Advection rate 0", "", user.a[0], &user.a[0], NULL));
72 user.a[1] = 0;
73 PetscCall(PetscOptionsReal("-a1", "Advection rate 1", "", user.a[1], &user.a[1], NULL));
74 user.k[0] = 1e6;
75 PetscCall(PetscOptionsReal("-k0", "Reaction rate 0", "", user.k[0], &user.k[0], NULL));
76 user.k[1] = 2 * user.k[0];
77 PetscCall(PetscOptionsReal("-k1", "Reaction rate 1", "", user.k[1], &user.k[1], NULL));
78 user.s[0] = 0;
79 PetscCall(PetscOptionsReal("-s0", "Source 0", "", user.s[0], &user.s[0], NULL));
80 user.s[1] = 1;
81 PetscCall(PetscOptionsReal("-s1", "Source 1", "", user.s[1], &user.s[1], NULL));
82 }
83 PetscOptionsEnd();
84
85 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
86 Create timestepping solver context
87 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
89 PetscCall(TSSetDM(ts, da));
90 PetscCall(TSSetType(ts, TSARKIMEX));
91 PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, &user));
92 PetscCall(TSSetIFunction(ts, NULL, FormIFunction, &user));
93 PetscCall(DMSetMatType(da, MATAIJ));
94 PetscCall(DMCreateMatrix(da, &J));
95 PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, &user));
96
97 /* A line search in the nonlinear solve can fail due to ill-conditioning unless an absolute tolerance is set. Since
98 * this problem is linear, we deactivate the line search. For a linear problem, it is usually recommended to also use
99 * SNESSetType(snes,SNESKSPONLY). */
100 PetscCall(TSGetSNES(ts, &snes));
101 PetscCall(SNESGetLineSearch(snes, &linesearch));
102 PetscCall(SNESLineSearchSetType(linesearch, SNESLINESEARCHBASIC));
103
104 ftime = .1;
105 PetscCall(TSSetMaxTime(ts, ftime));
106 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
107
108 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
109 Set initial conditions
110 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111 PetscCall(FormInitialSolution(ts, X, &user));
112 PetscCall(TSSetSolution(ts, X));
113 PetscCall(VecGetSize(X, &mx));
114 dt = .1 * PetscMax(user.a[0], user.a[1]) / mx; /* Advective CFL, I don't know why it needs so much safety factor. */
115 PetscCall(TSSetTimeStep(ts, dt));
116
117 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118 Set runtime options
119 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120 PetscCall(TSSetFromOptions(ts));
121
122 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123 Solve nonlinear system
124 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
125 PetscCall(TSSolve(ts, X));
126 PetscCall(TSGetSolveTime(ts, &ftime));
127 PetscCall(TSGetStepNumber(ts, &steps));
128 PetscCall(TSGetConvergedReason(ts, &reason));
129 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps));
130
131 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 Free work space.
133 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134 PetscCall(MatDestroy(&J));
135 PetscCall(VecDestroy(&X));
136 PetscCall(TSDestroy(&ts));
137 PetscCall(DMDestroy(&da));
138 PetscCall(PetscFinalize());
139 return 0;
140 }
141
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void * ptr)142 static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
143 {
144 User user = (User)ptr;
145 DM da;
146 DMDALocalInfo info;
147 PetscInt i;
148 Field *f;
149 const Field *x, *xdot;
150
151 PetscFunctionBeginUser;
152 PetscCall(TSGetDM(ts, &da));
153 PetscCall(DMDAGetLocalInfo(da, &info));
154
155 /* Get pointers to vector data */
156 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
157 PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot));
158 PetscCall(DMDAVecGetArray(da, F, &f));
159
160 /* Compute function over the locally owned part of the grid */
161 for (i = info.xs; i < info.xs + info.xm; i++) {
162 f[i][0] = xdot[i][0] + user->k[0] * x[i][0] - user->k[1] * x[i][1] - user->s[0];
163 f[i][1] = xdot[i][1] - user->k[0] * x[i][0] + user->k[1] * x[i][1] - user->s[1];
164 }
165
166 /* Restore vectors */
167 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));
168 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot));
169 PetscCall(DMDAVecRestoreArray(da, F, &f));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void * ptr)173 static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
174 {
175 User user = (User)ptr;
176 DM da;
177 Vec Xloc;
178 DMDALocalInfo info;
179 PetscInt i, j;
180 PetscReal hx;
181 Field *f;
182 const Field *x;
183
184 PetscFunctionBeginUser;
185 PetscCall(TSGetDM(ts, &da));
186 PetscCall(DMDAGetLocalInfo(da, &info));
187 hx = 1.0 / (PetscReal)info.mx;
188
189 /*
190 Scatter ghost points to local vector,using the 2-step process
191 DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
192 By placing code between these two statements, computations can be
193 done while messages are in transition.
194 */
195 PetscCall(DMGetLocalVector(da, &Xloc));
196 PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
197 PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
198
199 /* Get pointers to vector data */
200 PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
201 PetscCall(DMDAVecGetArray(da, F, &f));
202
203 /* Compute function over the locally owned part of the grid */
204 for (i = info.xs; i < info.xs + info.xm; i++) {
205 const PetscReal *a = user->a;
206 PetscReal u0t[2];
207 u0t[0] = 1.0 - PetscPowRealInt(PetscSinReal(12 * t), 4);
208 u0t[1] = 0.0;
209 for (j = 0; j < 2; j++) {
210 if (i == 0) f[i][j] = a[j] / hx * (1. / 3 * u0t[j] + 0.5 * x[i][j] - x[i + 1][j] + 1. / 6 * x[i + 2][j]);
211 else if (i == 1) f[i][j] = a[j] / hx * (-1. / 12 * u0t[j] + 2. / 3 * x[i - 1][j] - 2. / 3 * x[i + 1][j] + 1. / 12 * x[i + 2][j]);
212 else if (i == info.mx - 2) f[i][j] = a[j] / hx * (-1. / 6 * x[i - 2][j] + x[i - 1][j] - 0.5 * x[i][j] - 1. / 3 * x[i + 1][j]);
213 else if (i == info.mx - 1) f[i][j] = a[j] / hx * (-x[i][j] + x[i - 1][j]);
214 else f[i][j] = a[j] / hx * (-1. / 12 * x[i - 2][j] + 2. / 3 * x[i - 1][j] - 2. / 3 * x[i + 1][j] + 1. / 12 * x[i + 2][j]);
215 }
216 }
217
218 /* Restore vectors */
219 PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
220 PetscCall(DMDAVecRestoreArray(da, F, &f));
221 PetscCall(DMRestoreLocalVector(da, &Xloc));
222 PetscFunctionReturn(PETSC_SUCCESS);
223 }
224
225 /* --------------------------------------------------------------------- */
226 /*
227 IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
228 */
FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void * ptr)229 PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
230 {
231 User user = (User)ptr;
232 DMDALocalInfo info;
233 PetscInt i;
234 DM da;
235 const Field *x, *xdot;
236
237 PetscFunctionBeginUser;
238 PetscCall(TSGetDM(ts, &da));
239 PetscCall(DMDAGetLocalInfo(da, &info));
240
241 /* Get pointers to vector data */
242 PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
243 PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot));
244
245 /* Compute function over the locally owned part of the grid */
246 for (i = info.xs; i < info.xs + info.xm; i++) {
247 const PetscReal *k = user->k;
248 PetscScalar v[2][2];
249
250 v[0][0] = a + k[0];
251 v[0][1] = -k[1];
252 v[1][0] = -k[0];
253 v[1][1] = a + k[1];
254 PetscCall(MatSetValuesBlocked(Jpre, 1, &i, 1, &i, &v[0][0], INSERT_VALUES));
255 }
256
257 /* Restore vectors */
258 PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));
259 PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot));
260
261 PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
262 PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
263 if (J != Jpre) {
264 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
265 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
266 }
267 PetscFunctionReturn(PETSC_SUCCESS);
268 }
269
FormInitialSolution(TS ts,Vec X,PetscCtx ctx)270 PetscErrorCode FormInitialSolution(TS ts, Vec X, PetscCtx ctx)
271 {
272 User user = (User)ctx;
273 DM da;
274 PetscInt i;
275 DMDALocalInfo info;
276 Field *x;
277 PetscReal hx;
278
279 PetscFunctionBeginUser;
280 PetscCall(TSGetDM(ts, &da));
281 PetscCall(DMDAGetLocalInfo(da, &info));
282 hx = 1.0 / (PetscReal)info.mx;
283
284 /* Get pointers to vector data */
285 PetscCall(DMDAVecGetArray(da, X, &x));
286
287 /* Compute function over the locally owned part of the grid */
288 for (i = info.xs; i < info.xs + info.xm; i++) {
289 PetscReal r = (i + 1) * hx, ik = user->k[1] != 0.0 ? 1.0 / user->k[1] : 1.0;
290 x[i][0] = 1 + user->s[1] * r;
291 x[i][1] = user->k[0] * ik * x[i][0] + user->s[1] * ik;
292 }
293 PetscCall(DMDAVecRestoreArray(da, X, &x));
294 PetscFunctionReturn(PETSC_SUCCESS);
295 }
296
297 /*TEST
298
299 test:
300 args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_arkimex_type 4 -ts_adapt_type none -ts_time_step .005 -ts_max_time .1
301 requires: !single
302
303 test:
304 suffix: 2
305 args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_time_step 1e-3 -ts_adapt_type none -ts_time_step .005 -ts_max_time .1
306 nsize: 2
307
308 test:
309 suffix: 3
310 args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type ra34pw2 -ts_time_step 5e-3 -ts_max_time .1 -ts_adapt_type none
311 nsize: 2
312
313 test:
314 suffix: 4
315 args: -ts_type eimex -da_grid_x 200 -ts_eimex_order_adapt true -ts_time_step 0.001 -ts_monitor -ts_max_steps 100
316 filter: sed "s/ITS/TIME/g"
317 nsize: 2
318
319 TEST*/
320