Actual source code: ex35.cxx

  1: static const char help[] = "Time-dependent Brusselator reaction-diffusion PDE in 1d. Demonstrates IMEX methods and uses MOAB.\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

  8:    Initial conditions:
  9:    u(x,0) = 1 + sin(2 pi x)
 10:    v(x,0) = 3

 12:    Boundary conditions:
 13:    u(0,t) = u(1,t) = 1
 14:    v(0,t) = v(1,t) = 3
 15: */

 17: // PETSc includes:
 18: #include <petscts.h>
 19: #include <petscdmmoab.h>

 21: typedef struct {
 22:   PetscScalar u, v;
 23: } Field;

 25: struct pUserCtx {
 26:   PetscReal A, B;    /* Reaction coefficients */
 27:   PetscReal alpha;   /* Diffusion coefficient */
 28:   Field     leftbc;  /* Dirichlet boundary conditions at left boundary */
 29:   Field     rightbc; /* Dirichlet boundary conditions at right boundary */
 30:   PetscInt  n, npts; /* Number of mesh points */
 31:   PetscInt  ntsteps; /* Number of time steps */
 32:   PetscInt  nvars;   /* Number of variables in the equation system */
 33:   PetscBool io;
 34: };
 35: typedef pUserCtx *UserCtx;

 37: PetscErrorCode Initialize_AppContext(UserCtx *puser)
 38: {
 39:   UserCtx user;

 41:   PetscNew(&user);
 42:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Advection-reaction options", "ex35.cxx");
 43:   {
 44:     user->nvars     = 2;
 45:     user->A         = 1;
 46:     user->B         = 3;
 47:     user->alpha     = 0.02;
 48:     user->leftbc.u  = 1;
 49:     user->rightbc.u = 1;
 50:     user->leftbc.v  = 3;
 51:     user->rightbc.v = 3;
 52:     user->n         = 10;
 53:     user->ntsteps   = 10000;
 54:     user->io        = PETSC_FALSE;
 55:     PetscOptionsReal("-A", "Reaction rate", "ex35.cxx", user->A, &user->A, NULL);
 56:     PetscOptionsReal("-B", "Reaction rate", "ex35.cxx", user->B, &user->B, NULL);
 57:     PetscOptionsReal("-alpha", "Diffusion coefficient", "ex35.cxx", user->alpha, &user->alpha, NULL);
 58:     PetscOptionsScalar("-uleft", "Dirichlet boundary condition", "ex35.cxx", user->leftbc.u, &user->leftbc.u, NULL);
 59:     PetscOptionsScalar("-uright", "Dirichlet boundary condition", "ex35.cxx", user->rightbc.u, &user->rightbc.u, NULL);
 60:     PetscOptionsScalar("-vleft", "Dirichlet boundary condition", "ex35.cxx", user->leftbc.v, &user->leftbc.v, NULL);
 61:     PetscOptionsScalar("-vright", "Dirichlet boundary condition", "ex35.cxx", user->rightbc.v, &user->rightbc.v, NULL);
 62:     PetscOptionsInt("-n", "Number of 1-D elements", "ex35.cxx", user->n, &user->n, NULL);
 63:     PetscOptionsInt("-ndt", "Number of time steps", "ex35.cxx", user->ntsteps, &user->ntsteps, NULL);
 64:     PetscOptionsBool("-io", "Write the mesh and solution output to a file.", "ex35.cxx", user->io, &user->io, NULL);
 65:     user->npts = user->n + 1;
 66:   }
 67:   PetscOptionsEnd();

 69:   *puser = user;
 70:   return 0;
 71: }

 73: PetscErrorCode Destroy_AppContext(UserCtx *user)
 74: {
 75:   PetscFree(*user);
 76:   return 0;
 77: }

 79: static PetscErrorCode FormInitialSolution(TS, Vec, void *);
 80: static PetscErrorCode FormRHSFunction(TS, PetscReal, Vec, Vec, void *);
 81: static PetscErrorCode FormIFunction(TS, PetscReal, Vec, Vec, Vec, void *);
 82: static PetscErrorCode FormIJacobian(TS, PetscReal, Vec, Vec, PetscReal, Mat, Mat, void *);

 84: /****************
 85:  *              *
 86:  *     MAIN     *
 87:  *              *
 88:  ****************/
 89: int main(int argc, char **argv)
 90: {
 91:   TS                ts; /* nonlinear solver */
 92:   Vec               X;  /* solution, residual vectors */
 93:   Mat               J;  /* Jacobian matrix */
 94:   PetscInt          steps, mx;
 95:   PetscReal         hx, dt, ftime;
 96:   UserCtx           user; /* user-defined work context */
 97:   TSConvergedReason reason;
 98:   DM                dm;
 99:   const char       *fields[2] = {"U", "V"};

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

104:   /* Initialize the user context struct */
105:   Initialize_AppContext(&user);

107:   /* Fill in the user defined work context: */
108:   DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, NULL, user->n, 1, &dm);
109:   DMMoabSetFieldNames(dm, user->nvars, fields);
110:   DMMoabSetBlockSize(dm, user->nvars);
111:   DMSetFromOptions(dm);

113:   /* SetUp the data structures for DMMOAB */
114:   DMSetUp(dm);

116:   /*  Create timestepping solver context */
117:   TSCreate(PETSC_COMM_WORLD, &ts);
118:   TSSetDM(ts, dm);
119:   TSSetType(ts, TSARKIMEX);
120:   TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1);
121:   DMSetMatType(dm, MATBAIJ);
122:   DMCreateMatrix(dm, &J);

124:   TSSetRHSFunction(ts, NULL, FormRHSFunction, user);
125:   TSSetIFunction(ts, NULL, FormIFunction, user);
126:   TSSetIJacobian(ts, J, J, FormIJacobian, user);

128:   ftime = 10.0;
129:   TSSetMaxSteps(ts, user->ntsteps);
130:   TSSetMaxTime(ts, ftime);
131:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

133:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
134:      Create the solution vector and set the initial conditions
135:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136:   DMCreateGlobalVector(dm, &X);

138:   FormInitialSolution(ts, X, user);
139:   TSSetSolution(ts, X);
140:   VecGetSize(X, &mx);
141:   hx = 1.0 / (PetscReal)(mx / 2 - 1);
142:   dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */
143:   TSSetTimeStep(ts, dt);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:      Set runtime options
147:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148:   TSSetFromOptions(ts);

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:      Solve nonlinear system
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:   TSSolve(ts, X);
154:   TSGetSolveTime(ts, &ftime);
155:   TSGetStepNumber(ts, &steps);
156:   TSGetConvergedReason(ts, &reason);
157:   PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, steps);

159:   if (user->io) {
160:     /* Print the numerical solution to screen and then dump to file */
161:     VecView(X, PETSC_VIEWER_STDOUT_WORLD);

163:     /* Write out the solution along with the mesh */
164:     DMMoabSetGlobalFieldVector(dm, X);
165: #ifdef MOAB_HAVE_HDF5
166:     DMMoabOutput(dm, "ex35.h5m", "");
167: #else
168:     /* MOAB does not support true parallel writers that aren't HDF5 based
169:        And so if you are using VTK as the output format in parallel,
170:        the data could be jumbled due to the order in which the processors
171:        write out their parts of the mesh and solution tags
172:     */
173:     DMMoabOutput(dm, "ex35.vtk", "");
174: #endif
175:   }

177:   /* Free work space.
178:      Free all PETSc related resources: */
179:   MatDestroy(&J);
180:   VecDestroy(&X);
181:   TSDestroy(&ts);
182:   DMDestroy(&dm);

184:   /* Free all MOAB related resources: */
185:   Destroy_AppContext(&user);
186:   PetscFinalize();
187:   return 0;
188: }

190: /*
191:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
192: */
193: PetscErrorCode FormIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat J, Mat Jpre, void *ptr)
194: {
195:   UserCtx            user = (UserCtx)ptr;
196:   PetscInt           dof;
197:   PetscReal          hx;
198:   DM                 dm;
199:   const moab::Range *vlocal;
200:   PetscBool          vonboundary;

202:   TSGetDM(ts, &dm);

204:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
205:   DMMoabGetLocalVertices(dm, &vlocal, NULL);

207:   /* compute local element sizes - structured grid */
208:   hx = 1.0 / user->n;

210:   /* Compute function over the locally owned part of the grid
211:      Assemble the operator by looping over edges and computing
212:      contribution for each vertex dof                         */
213:   for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
214:     const moab::EntityHandle vhandle = *iter;

216:     DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof);

218:     /* check if vertex is on the boundary */
219:     DMMoabIsEntityOnBoundary(dm, vhandle, &vonboundary);

221:     if (vonboundary) {
222:       const PetscScalar bcvals[2][2] = {
223:         {hx, 0 },
224:         {0,  hx}
225:       };
226:       MatSetValuesBlocked(Jpre, 1, &dof, 1, &dof, &bcvals[0][0], INSERT_VALUES);
227:     } else {
228:       const PetscInt    row = dof, col[] = {dof - 1, dof, dof + 1};
229:       const PetscScalar dxxL = -user->alpha / hx, dxx0 = 2. * user->alpha / hx, dxxR = -user->alpha / hx;
230:       const PetscScalar vals[2][3][2] = {
231:         {{dxxL, 0}, {a * hx + dxx0, 0}, {dxxR, 0}},
232:         {{0, dxxL}, {0, a * hx + dxx0}, {0, dxxR}}
233:       };
234:       MatSetValuesBlocked(Jpre, 1, &row, 3, col, &vals[0][0][0], INSERT_VALUES);
235:     }
236:   }

238:   MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY);
239:   MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY);
240:   if (J != Jpre) {
241:     MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
242:     MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
243:   }
244:   return 0;
245: }

247: static PetscErrorCode FormRHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr)
248: {
249:   UserCtx            user = (UserCtx)ptr;
250:   DM                 dm;
251:   PetscReal          hx;
252:   const Field       *x;
253:   Field             *f;
254:   PetscInt           dof;
255:   const moab::Range *ownedvtx;

257:   hx = 1.0 / user->n;
258:   TSGetDM(ts, &dm);

260:   /* Get pointers to vector data */
261:   VecSet(F, 0.0);

263:   DMMoabVecGetArrayRead(dm, X, &x);
264:   DMMoabVecGetArray(dm, F, &f);

266:   DMMoabGetLocalVertices(dm, &ownedvtx, NULL);

268:   /* Compute function over the locally owned part of the grid */
269:   for (moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
270:     const moab::EntityHandle vhandle = *iter;
271:     DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof);

273:     PetscScalar u = x[dof].u, v = x[dof].v;
274:     f[dof].u = hx * (user->A + u * u * v - (user->B + 1) * u);
275:     f[dof].v = hx * (user->B * u - u * u * v);
276:   }

278:   /* Restore vectors */
279:   DMMoabVecRestoreArrayRead(dm, X, &x);
280:   DMMoabVecRestoreArray(dm, F, &f);
281:   return 0;
282: }

284: static PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
285: {
286:   UserCtx            user = (UserCtx)ctx;
287:   DM                 dm;
288:   Field             *x, *xdot, *f;
289:   PetscReal          hx;
290:   Vec                Xloc;
291:   PetscInt           i, bcindx;
292:   PetscBool          elem_on_boundary;
293:   const moab::Range *vlocal;

295:   hx = 1.0 / user->n;
296:   TSGetDM(ts, &dm);

298:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
299:   DMMoabGetLocalVertices(dm, &vlocal, NULL);

301:   /* reset the residual vector */
302:   VecSet(F, 0.0);

304:   DMGetLocalVector(dm, &Xloc);
305:   DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc);
306:   DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc);

308:   /* get the local representation of the arrays from Vectors */
309:   DMMoabVecGetArrayRead(dm, Xloc, &x);
310:   DMMoabVecGetArrayRead(dm, Xdot, &xdot);
311:   DMMoabVecGetArray(dm, F, &f);

313:   /* loop over local elements */
314:   for (moab::Range::iterator iter = vlocal->begin(); iter != vlocal->end(); iter++) {
315:     const moab::EntityHandle vhandle = *iter;

317:     DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &i);

319:     /* check if vertex is on the boundary */
320:     DMMoabIsEntityOnBoundary(dm, vhandle, &elem_on_boundary);

322:     if (elem_on_boundary) {
323:       DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx);
324:       if (bcindx == 0) { /* Apply left BC */
325:         f[i].u = hx * (x[i].u - user->leftbc.u);
326:         f[i].v = hx * (x[i].v - user->leftbc.v);
327:       } else { /* Apply right BC */
328:         f[i].u = hx * (x[i].u - user->rightbc.u);
329:         f[i].v = hx * (x[i].v - user->rightbc.v);
330:       }
331:     } else {
332:       f[i].u = hx * xdot[i].u - user->alpha * (x[i - 1].u - 2. * x[i].u + x[i + 1].u) / hx;
333:       f[i].v = hx * xdot[i].v - user->alpha * (x[i - 1].v - 2. * x[i].v + x[i + 1].v) / hx;
334:     }
335:   }

337:   /* Restore data */
338:   DMMoabVecRestoreArrayRead(dm, Xloc, &x);
339:   DMMoabVecRestoreArrayRead(dm, Xdot, &xdot);
340:   DMMoabVecRestoreArray(dm, F, &f);
341:   DMRestoreLocalVector(dm, &Xloc);
342:   return 0;
343: }

345: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
346: {
347:   UserCtx               user = (UserCtx)ctx;
348:   PetscReal             vpos[3];
349:   DM                    dm;
350:   Field                *x;
351:   const moab::Range    *vowned;
352:   PetscInt              dof;
353:   moab::Range::iterator iter;

355:   TSGetDM(ts, &dm);

357:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
358:   DMMoabGetLocalVertices(dm, &vowned, NULL);

360:   VecSet(X, 0.0);

362:   /* Get pointers to vector data */
363:   DMMoabVecGetArray(dm, X, &x);

365:   /* Compute function over the locally owned part of the grid */
366:   for (moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) {
367:     const moab::EntityHandle vhandle = *iter;
368:     DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof);

370:     /* compute the mid-point of the element and use a 1-point lumped quadrature */
371:     DMMoabGetVertexCoordinates(dm, 1, &vhandle, vpos);

373:     PetscReal xi = vpos[0];
374:     x[dof].u     = user->leftbc.u * (1. - xi) + user->rightbc.u * xi + PetscSinReal(2. * PETSC_PI * xi);
375:     x[dof].v     = user->leftbc.v * (1. - xi) + user->rightbc.v * xi;
376:   }

378:   /* Restore vectors */
379:   DMMoabVecRestoreArray(dm, X, &x);
380:   return 0;
381: }

383: /*TEST

385:     build:
386:       requires: moab

388:     test:
389:       args: -n 20 -ts_type rosw -ts_rosw_type 2p -ts_dt 5e-2 -ts_adapt_type none

391:     test:
392:       suffix: 2
393:       nsize: 2
394:       args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io
395:       TODO:

397: TEST*/