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

 70:   *puser = user;
 71:   PetscFunctionReturn(PETSC_SUCCESS);
 72: }

 74: PetscErrorCode Destroy_AppContext(UserCtx *user)
 75: {
 76:   PetscFunctionBegin;
 77:   PetscCall(PetscFree(*user));
 78:   PetscFunctionReturn(PETSC_SUCCESS);
 79: }

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

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

103:   PetscFunctionBeginUser;
104:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

106:   /* Initialize the user context struct */
107:   PetscCall(Initialize_AppContext(&user));

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

115:   /* SetUp the data structures for DMMOAB */
116:   PetscCall(DMSetUp(dm));

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

126:   PetscCall(TSSetRHSFunction(ts, NULL, FormRHSFunction, user));
127:   PetscCall(TSSetIFunction(ts, NULL, FormIFunction, user));
128:   PetscCall(TSSetIJacobian(ts, J, J, FormIJacobian, user));

130:   ftime = 10.0;
131:   PetscCall(TSSetMaxSteps(ts, user->ntsteps));
132:   PetscCall(TSSetMaxTime(ts, ftime));
133:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

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

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

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:      Set runtime options
149:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150:   PetscCall(TSSetFromOptions(ts));

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

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

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

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

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

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

204:   PetscFunctionBegin;
205:   PetscCall(TSGetDM(ts, &dm));

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

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

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

219:     PetscCall(DMMoabGetDofsBlocked(dm, 1, &vhandle, &dof));

221:     /* check if vertex is on the boundary */
222:     PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &vonboundary));

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

241:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
242:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
243:   if (J != Jpre) {
244:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
245:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
246:   }
247:   PetscFunctionReturn(PETSC_SUCCESS);
248: }

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

260:   PetscFunctionBegin;
261:   hx = 1.0 / user->n;
262:   PetscCall(TSGetDM(ts, &dm));

264:   /* Get pointers to vector data */
265:   PetscCall(VecSet(F, 0.0));

267:   PetscCall(DMMoabVecGetArrayRead(dm, X, &x));
268:   PetscCall(DMMoabVecGetArray(dm, F, &f));

270:   PetscCall(DMMoabGetLocalVertices(dm, &ownedvtx, NULL));

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

277:     PetscScalar u = x[dof].u, v = x[dof].v;
278:     f[dof].u = hx * (user->A + u * u * v - (user->B + 1) * u);
279:     f[dof].v = hx * (user->B * u - u * u * v);
280:   }

282:   /* Restore vectors */
283:   PetscCall(DMMoabVecRestoreArrayRead(dm, X, &x));
284:   PetscCall(DMMoabVecRestoreArray(dm, F, &f));
285:   PetscFunctionReturn(PETSC_SUCCESS);
286: }

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

299:   PetscFunctionBegin;
300:   hx = 1.0 / user->n;
301:   PetscCall(TSGetDM(ts, &dm));

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

306:   /* reset the residual vector */
307:   PetscCall(VecSet(F, 0.0));

309:   PetscCall(DMGetLocalVector(dm, &Xloc));
310:   PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, Xloc));
311:   PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, Xloc));

313:   /* get the local representation of the arrays from Vectors */
314:   PetscCall(DMMoabVecGetArrayRead(dm, Xloc, &x));
315:   PetscCall(DMMoabVecGetArrayRead(dm, Xdot, &xdot));
316:   PetscCall(DMMoabVecGetArray(dm, F, &f));

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

322:     PetscCall(DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &i));

324:     /* check if vertex is on the boundary */
325:     PetscCall(DMMoabIsEntityOnBoundary(dm, vhandle, &elem_on_boundary));

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

342:   /* Restore data */
343:   PetscCall(DMMoabVecRestoreArrayRead(dm, Xloc, &x));
344:   PetscCall(DMMoabVecRestoreArrayRead(dm, Xdot, &xdot));
345:   PetscCall(DMMoabVecRestoreArray(dm, F, &f));
346:   PetscCall(DMRestoreLocalVector(dm, &Xloc));
347:   PetscFunctionReturn(PETSC_SUCCESS);
348: }

350: PetscErrorCode FormInitialSolution(TS ts, Vec X, void *ctx)
351: {
352:   UserCtx               user = (UserCtx)ctx;
353:   PetscReal             vpos[3];
354:   DM                    dm;
355:   Field                *x;
356:   const moab::Range    *vowned;
357:   PetscInt              dof;
358:   moab::Range::iterator iter;

360:   PetscFunctionBegin;
361:   PetscCall(TSGetDM(ts, &dm));

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

366:   PetscCall(VecSet(X, 0.0));

368:   /* Get pointers to vector data */
369:   PetscCall(DMMoabVecGetArray(dm, X, &x));

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

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

379:     PetscReal xi = vpos[0];
380:     x[dof].u     = user->leftbc.u * (1. - xi) + user->rightbc.u * xi + PetscSinReal(2. * PETSC_PI * xi);
381:     x[dof].v     = user->leftbc.v * (1. - xi) + user->rightbc.v * xi;
382:   }

384:   /* Restore vectors */
385:   PetscCall(DMMoabVecRestoreArray(dm, X, &x));
386:   PetscFunctionReturn(PETSC_SUCCESS);
387: }

389: /*TEST

391:     build:
392:       requires: moab

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

397:     test:
398:       suffix: 2
399:       nsize: 2
400:       args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io
401:       TODO:

403: TEST*/