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*/