Actual source code: ex25.c

  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).

  8:    Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and multiple solves.

 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

 15:    Initial conditions:
 16:    u(x,0) = 1 + sin(2 pi x)
 17:    v(x,0) = 3

 19:    Boundary conditions:
 20:    u(0,t) = u(1,t) = 1
 21:    v(0,t) = v(1,t) = 3
 22: */

 24: #include <petscdm.h>
 25: #include <petscdmda.h>
 26: #include <petscts.h>

 28: typedef struct {
 29:   PetscScalar u, v;
 30: } Field;

 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: };

 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 int            Brusselator(int, char **, PetscInt);

 46: int main(int argc, char **argv)
 47: {
 48:   PetscInt       cycle;

 51:   MPI_Init(&argc, &argv);
 52:   if (ierr) return ierr;
 53:   for (cycle = 0; cycle < 4; cycle++) {
 54:     Brusselator(argc, argv, cycle);
 55:     if (ierr) return 1;
 56:   }
 57:   MPI_Finalize();
 58:   return ierr;
 59: }

 61: PetscErrorCode Brusselator(int argc, char **argv, PetscInt cycle)
 62: {
 63:   TS                ts; /* nonlinear solver */
 64:   Vec               X;  /* solution, residual vectors */
 65:   Mat               J;  /* Jacobian matrix */
 66:   PetscInt          steps, mx;
 67:   DM                da;
 68:   PetscReal         ftime, hx, dt, xmax, xmin;
 69:   struct _User      user; /* user-defined work context */
 70:   TSConvergedReason reason;

 73:   PetscInitialize(&argc, &argv, (char *)0, help);

 75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 76:      Create distributed array (DMDA) to manage parallel grid and vectors
 77:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 78:   DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, 11, 2, 2, NULL, &da);
 79:   DMSetFromOptions(da);
 80:   DMSetUp(da);

 82:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 83:      Extract global vectors from DMDA;
 84:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 85:   DMCreateGlobalVector(da, &X);

 87:   /* Initialize user application context */
 88:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "");
 89:   {
 90:     user.A      = 1;
 91:     user.B      = 3;
 92:     user.alpha  = 0.1;
 93:     user.uleft  = 1;
 94:     user.uright = 1;
 95:     user.vleft  = 3;
 96:     user.vright = 3;
 97:     PetscOptionsReal("-A", "Reaction rate", "", user.A, &user.A, NULL);
 98:     PetscOptionsReal("-B", "Reaction rate", "", user.B, &user.B, NULL);
 99:     PetscOptionsReal("-alpha", "Diffusion coefficient", "", user.alpha, &user.alpha, NULL);
100:     PetscOptionsReal("-uleft", "Dirichlet boundary condition", "", user.uleft, &user.uleft, NULL);
101:     PetscOptionsReal("-uright", "Dirichlet boundary condition", "", user.uright, &user.uright, NULL);
102:     PetscOptionsReal("-vleft", "Dirichlet boundary condition", "", user.vleft, &user.vleft, NULL);
103:     PetscOptionsReal("-vright", "Dirichlet boundary condition", "", user.vright, &user.vright, NULL);
104:   }
105:   PetscOptionsEnd();

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:      Create timestepping solver context
109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
110:   TSCreate(PETSC_COMM_WORLD, &ts);
111:   TSSetDM(ts, da);
112:   TSSetType(ts, TSARKIMEX);
113:   TSSetRHSFunction(ts, NULL, FormRHSFunction, &user);
114:   TSSetIFunction(ts, NULL, FormIFunction, &user);
115:   DMSetMatType(da, MATAIJ);
116:   DMCreateMatrix(da, &J);
117:   TSSetIJacobian(ts, J, J, FormIJacobian, &user);

119:   ftime = 1.0;
120:   TSSetMaxTime(ts, ftime);
121:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

123:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124:      Set initial conditions
125:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126:   FormInitialSolution(ts, X, &user);
127:   TSSetSolution(ts, X);
128:   VecGetSize(X, &mx);
129:   hx = 1.0 / (PetscReal)(mx / 2 - 1);
130:   dt = 0.4 * PetscSqr(hx) / user.alpha; /* Diffusive stability limit */
131:   dt *= PetscPowRealInt(0.2, cycle);    /* Shrink the time step in convergence study. */
132:   TSSetTimeStep(ts, dt);
133:   TSSetTolerances(ts, 1e-3 * PetscPowRealInt(0.5, cycle), NULL, 1e-3 * PetscPowRealInt(0.5, cycle), NULL);

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:      Set runtime options
137:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138:   TSSetFromOptions(ts);

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:      Solve nonlinear system
142:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143:   TSSolve(ts, X);
144:   TSGetSolveTime(ts, &ftime);
145:   TSGetStepNumber(ts, &steps);
146:   TSGetConvergedReason(ts, &reason);
147:   VecMin(X, NULL, &xmin);
148:   VecMax(X, NULL, &xmax);
149:   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);

151:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152:      Free work space.
153:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154:   MatDestroy(&J);
155:   VecDestroy(&X);
156:   TSDestroy(&ts);
157:   DMDestroy(&da);
158:   PetscFinalize();
159:   return 0;
160: }

162: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ptr)
163: {
164:   User          user = (User)ptr;
165:   DM            da;
166:   DMDALocalInfo info;
167:   PetscInt      i;
168:   Field        *x, *xdot, *f;
169:   PetscReal     hx;
170:   Vec           Xloc;

173:   TSGetDM(ts, &da);
174:   DMDAGetLocalInfo(da, &info);
175:   hx = 1.0 / (PetscReal)(info.mx - 1);

177:   /*
178:      Scatter ghost points to local vector,using the 2-step process
179:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
180:      By placing code between these two statements, computations can be
181:      done while messages are in transition.
182:   */
183:   DMGetLocalVector(da, &Xloc);
184:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
185:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);

187:   /* Get pointers to vector data */
188:   DMDAVecGetArrayRead(da, Xloc, &x);
189:   DMDAVecGetArrayRead(da, Xdot, &xdot);
190:   DMDAVecGetArray(da, F, &f);

192:   /* Compute function over the locally owned part of the grid */
193:   for (i = info.xs; i < info.xs + info.xm; i++) {
194:     if (i == 0) {
195:       f[i].u = hx * (x[i].u - user->uleft);
196:       f[i].v = hx * (x[i].v - user->vleft);
197:     } else if (i == info.mx - 1) {
198:       f[i].u = hx * (x[i].u - user->uright);
199:       f[i].v = hx * (x[i].v - user->vright);
200:     } else {
201:       f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
202:       f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
203:     }
204:   }

206:   /* Restore vectors */
207:   DMDAVecRestoreArrayRead(da, Xloc, &x);
208:   DMDAVecRestoreArrayRead(da, Xdot, &xdot);
209:   DMDAVecRestoreArray(da, F, &f);
210:   DMRestoreLocalVector(da, &Xloc);
211:   return 0;
212: }

214: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
215: {
216:   User          user = (User)ptr;
217:   DM            da;
218:   DMDALocalInfo info;
219:   PetscInt      i;
220:   PetscReal     hx;
221:   Field        *x, *f;

224:   TSGetDM(ts, &da);
225:   DMDAGetLocalInfo(da, &info);
226:   hx = 1.0 / (PetscReal)(info.mx - 1);

228:   /* Get pointers to vector data */
229:   DMDAVecGetArrayRead(da, X, &x);
230:   DMDAVecGetArray(da, F, &f);

232:   /* Compute function over the locally owned part of the grid */
233:   for (i = info.xs; i < info.xs + info.xm; i++) {
234:     PetscScalar u = x[i].u, v = x[i].v;
235:     f[i].u = hx * (user->A + u * u * v - (user->B + 1) * u);
236:     f[i].v = hx * (user->B * u - u * u * v);
237:   }

239:   /* Restore vectors */
240:   DMDAVecRestoreArrayRead(da, X, &x);
241:   DMDAVecRestoreArray(da, F, &f);
242:   return 0;
243: }

245: /* --------------------------------------------------------------------- */
246: /*
247:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
248: */
249: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
250: {
251:   User          user = (User)ptr;
252:   DMDALocalInfo info;
253:   PetscInt      i;
254:   PetscReal     hx;
255:   DM            da;
256:   Field        *x, *xdot;

259:   TSGetDM(ts, &da);
260:   DMDAGetLocalInfo(da, &info);
261:   hx = 1.0 / (PetscReal)(info.mx - 1);

263:   /* Get pointers to vector data */
264:   DMDAVecGetArrayRead(da, X, &x);
265:   DMDAVecGetArrayRead(da, Xdot, &xdot);

267:   /* Compute function over the locally owned part of the grid */
268:   for (i = info.xs; i < info.xs + info.xm; i++) {
269:     if (i == 0 || i == info.mx - 1) {
270:       const PetscInt    row = i, col = i;
271:       const PetscScalar vals[2][2] = {
272:         {hx, 0 },
273:         {0,  hx}
274:       };
275:       MatSetValuesBlocked(Jpre, 1, &row, 1, &col, &vals[0][0], INSERT_VALUES);
276:     } else {
277:       const PetscInt    row = i, col[] = {i - 1, i, i + 1};
278:       const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
279:       const PetscScalar vals[2][3][2] = {
280:         {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
281:         {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
282:       };
283:       MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
284:     }
285:   }

287:   /* Restore vectors */
288:   DMDAVecRestoreArrayRead(da, X, &x);
289:   DMDAVecRestoreArrayRead(da, Xdot, &xdot);

291:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
292:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
293:   if (J != Jpre) {
294:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
295:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
296:   }
297:   return 0;
298: }

300: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
301: {
302:   User          user = (User)ctx;
303:   DM            da;
304:   PetscInt      i;
305:   DMDALocalInfo info;
306:   Field        *x;
307:   PetscReal     hx;

310:   TSGetDM(ts, &da);
311:   DMDAGetLocalInfo(da, &info);
312:   hx = 1.0 / (PetscReal)(info.mx - 1);

314:   /* Get pointers to vector data */
315:   DMDAVecGetArray(da, X, &x);

317:   /* Compute function over the locally owned part of the grid */
318:   for (i = info.xs; i < info.xs + info.xm; i++) {
319:     PetscReal xi = i * hx;
320:     x[i].u       = user->uleft * (1. - xi) + user->uright * xi + PetscSinReal(2. * PETSC_PI * xi);
321:     x[i].v       = user->vleft * (1. - xi) + user->vright * xi;
322:   }
323:   DMDAVecRestoreArray(da, X, &x);
324:   return 0;
325: }

327: /*TEST

329:     test:
330:       args: -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3

332:     test:
333:       suffix: 2
334:       args:   -ts_exact_final_time INTERPOLATE -snes_rtol 1.e-3

336: TEST*/