Actual source code: ex22.c
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
9: Initial conditions:
10: u(x,0) = 1 + s2*x
11: v(x,0) = k0/k1*u(x,0) + s1/k1
13: Upstream boundary conditions:
14: u(0,t) = 1-sin(12*t)^4
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: */
21: #include <petscdm.h>
22: #include <petscdmda.h>
23: #include <petscts.h>
25: typedef PetscScalar Field[2];
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: };
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 *);
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;
52: PetscFunctionBeginUser;
53: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
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));
62: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
63: Extract global vectors from DMDA;
64: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: PetscCall(DMCreateGlobalVector(da, &X));
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();
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));
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));
104: ftime = .1;
105: PetscCall(TSSetMaxTime(ts, ftime));
106: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
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));
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Set runtime options
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
120: PetscCall(TSSetFromOptions(ts));
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));
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: }
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;
151: PetscFunctionBeginUser;
152: PetscCall(TSGetDM(ts, &da));
153: PetscCall(DMDAGetLocalInfo(da, &info));
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));
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: }
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: }
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;
184: PetscFunctionBeginUser;
185: PetscCall(TSGetDM(ts, &da));
186: PetscCall(DMDAGetLocalInfo(da, &info));
187: hx = 1.0 / (PetscReal)info.mx;
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));
199: /* Get pointers to vector data */
200: PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
201: PetscCall(DMDAVecGetArray(da, F, &f));
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: }
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: }
225: /* --------------------------------------------------------------------- */
226: /*
227: IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
228: */
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;
237: PetscFunctionBeginUser;
238: PetscCall(TSGetDM(ts, &da));
239: PetscCall(DMDAGetLocalInfo(da, &info));
241: /* Get pointers to vector data */
242: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
243: PetscCall(DMDAVecGetArrayRead(da, Xdot, (void *)&xdot));
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];
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: }
257: /* Restore vectors */
258: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));
259: PetscCall(DMDAVecRestoreArrayRead(da, Xdot, (void *)&xdot));
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: }
270: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
271: {
272: User user = (User)ctx;
273: DM da;
274: PetscInt i;
275: DMDALocalInfo info;
276: Field *x;
277: PetscReal hx;
279: PetscFunctionBeginUser;
280: PetscCall(TSGetDM(ts, &da));
281: PetscCall(DMDAGetLocalInfo(da, &info));
282: hx = 1.0 / (PetscReal)info.mx;
284: /* Get pointers to vector data */
285: PetscCall(DMDAVecGetArray(da, X, &x));
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: }
297: /*TEST
299: test:
300: args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_arkimex_type 4 -ts_adapt_type none -ts_dt .005 -ts_max_time .1
301: requires: !single
303: test:
304: suffix: 2
305: args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_dt 1e-3 -ts_adapt_type none -ts_dt .005 -ts_max_time .1
306: nsize: 2
308: test:
309: suffix: 3
310: args: -nox -da_grid_x 200 -ts_monitor_draw_solution -ts_type rosw -ts_rosw_type ra34pw2 -ts_dt 5e-3 -ts_max_time .1 -ts_adapt_type none
311: nsize: 2
313: test:
314: suffix: 4
315: args: -ts_type eimex -da_grid_x 200 -ts_eimex_order_adapt true -ts_dt 0.001 -ts_monitor -ts_max_steps 100
316: filter: sed "s/ITS/TIME/g"
317: nsize: 2
319: TEST*/