Actual source code: ex35.cxx

petsc-3.7.3 2016-08-01
Report Typos and Errors
  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;

 39: PetscErrorCode Initialize_AppContext(UserCtx *puser)
 40: {
 41:   UserCtx           user;
 42:   PetscErrorCode    ierr;

 45:   PetscNew(&user);

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

 74:   *puser = user;
 75:   return(0);
 76: }

 80: PetscErrorCode Destroy_AppContext(UserCtx *user)
 81: {

 85:   PetscFree(*user);
 86:   return(0);
 87: }

 89: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
 90: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);

 92: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 93: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

 95: /****************
 96:  *              *
 97:  *     MAIN     *
 98:  *              *
 99:  ****************/
102: int main(int argc,char **argv)
103: {
104:   TS                ts;         /* nonlinear solver */
105:   Vec               X;          /* solution, residual vectors */
106:   Mat               J;          /* Jacobian matrix */
107:   PetscInt          steps,mx;
108:   PetscErrorCode    ierr;
109:   PetscReal         hx,dt,ftime;
110:   UserCtx           user;       /* user-defined work context */
111:   TSConvergedReason reason;

113:   DM                dm;
114:   moab::ErrorCode   merr;
115:   moab::Interface*  mbImpl;
116:   moab::Tag         solndofs;
117:   const moab::Range *ownedvtx;
118:   const PetscReal   bounds[2] = {0.0,1.0};
119:   const char        *fields[2] = {"U","V"};
120:   PetscScalar       deflt[2]={0.0,0.0};

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

124:   /* Initialize the user context struct */
125:   Initialize_AppContext(&user);

127:   /* Fill in the user defined work context: */
128:   DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, bounds, user->n, 1, &dm);
129:   DMMoabSetBlockSize(dm, user->nvars);
130:   DMMoabSetFieldNames(dm, user->nvars, fields);
131:   DMSetMatType(dm,MATBAIJ);
132:   DMSetFromOptions(dm);

134:   /* SetUp the data structures for DMMOAB */
135:   DMSetUp(dm);

137:   DMMoabGetInterface(dm, &mbImpl);

139:   /*  Create timestepping solver context */
140:   TSCreate(PETSC_COMM_WORLD,&ts);
141:   TSSetDM(ts, dm);
142:   TSSetType(ts,TSARKIMEX);
143:   TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);
144:   TSARKIMEXSetFullyImplicit(ts, PETSC_TRUE);
145:   TSSetRHSFunction(ts,NULL,FormRHSFunction,user);
146:   TSSetIFunction(ts,NULL,FormIFunction,user);
147:   DMCreateMatrix(dm,&J);
148:   TSSetIJacobian(ts,J,J,FormIJacobian,user);

150:   ftime = 10.0;
151:   TSSetDuration(ts,user->ntsteps,ftime);
152:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
153: 
154:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155:      Create the solution vector and set the initial conditions
156:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
157:   /* Use the call to DMMoabCreateVector for creating a named global MOAB Vec object.
158:      Alternately, use the following call to DM for creating an unnamed (anonymous) global 
159:      MOAB Vec object.

161:          DMCreateGlobalVector(dm, &X);
162:   */
163:   DMMoabGetLocalVertices(dm, &ownedvtx, NULL);
164:   /* create a tag to store the unknown fields in the problem */
165:   merr = mbImpl->tag_get_handle("UNKNOWNS",2,moab::MB_TYPE_DOUBLE,solndofs,
166:                                   moab::MB_TAG_DENSE|moab::MB_TAG_CREAT,deflt);MBERRNM(merr);

168:   DMMoabCreateVector(dm, solndofs, ownedvtx, PETSC_TRUE, PETSC_FALSE, &X);

170:   FormInitialSolution(ts,X,user);
171:   TSSetSolution(ts,X);
172:   VecGetSize(X,&mx);
173:   hx = 1.0/(PetscReal)(mx/2-1);
174:   dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */
175:   TSSetInitialTimeStep(ts,0.0,dt);

177:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178:      Set runtime options
179:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
180:   TSSetFromOptions(ts);

182:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183:      Solve nonlinear system
184:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185:   TSSolve(ts,X);
186:   TSGetSolveTime(ts,&ftime);
187:   TSGetTimeStepNumber(ts,&steps);
188:   TSGetConvergedReason(ts,&reason);
189:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],ftime,steps);

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

195:     /* Write out the solution along with the mesh */
196:     DMMoabSetGlobalFieldVector(dm, X);
197: #ifdef MOAB_HDF5_H
198:     DMMoabOutput(dm, "ex35.h5m", "");
199: #else
200:     /* MOAB does not support true parallel writers that aren't HDF5 based
201:        And so if you are using VTK as the output format in parallel,
202:        the data could be jumbled due to the order in which the processors
203:        write out their parts of the mesh and solution tags
204:     */
205:     DMMoabOutput(dm, "ex35.vtk", "");
206: #endif
207:   }

209:   /* Free work space.
210:      Free all PETSc related resources: */
211:   MatDestroy(&J);
212:   VecDestroy(&X);
213:   TSDestroy(&ts);
214:   DMDestroy(&dm);

216:   /* Free all MOAB related resources: */
217:   Destroy_AppContext(&user);

219:   PetscFinalize();
220:   return 0;
221: }

225: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
226: {
227:   UserCtx           user = (UserCtx)ptr;
228:   DM                dm;
229:   PetscReal         hx;
230:   const Field       *x;
231:   Field             *f;
232:   PetscInt          dof_index;
233:   const moab::Range *ownedvtx;
234:   PetscErrorCode    ierr;

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

240:   /* Get pointers to vector data */
241:   VecSet(F,0.0);

243:   DMMoabVecGetArrayRead(dm, X, &x);
244:   DMMoabVecGetArray(dm, F, &f);

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

248:   /* Compute function over the locally owned part of the grid */
249:   for(moab::Range::iterator iter = ownedvtx->begin(); iter != ownedvtx->end(); iter++) {
250:     const moab::EntityHandle vhandle = *iter;
251:     DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index);

253:     const Field& xx = x[dof_index];
254:     f[dof_index].u = hx*(user->A + xx.u*xx.u*xx.v - (user->B+1)*xx.u);
255:     f[dof_index].v = hx*(user->B*xx.u - xx.u*xx.u*xx.v);
256:   }

258:   /* Restore vectors */
259:   DMMoabVecRestoreArrayRead(dm, X, &x);
260:   DMMoabVecRestoreArray(dm, F, &f);
261:   return(0);
262: }


265: /*
266:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
267: */
270: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
271: {
272:   UserCtx             user = (UserCtx)ptr;
273:   PetscErrorCode      ierr;
274:   const moab::EntityHandle *connect;
275:   PetscInt            vpere=2;
276:   PetscReal           hx;
277:   DM                  dm;
278:   moab::Interface*    mbImpl;
279:   const moab::Range   *elocal;
280:   PetscInt            dof_indices[2];
281:   PetscBool           elem_on_boundary;

284:   TSGetDM(ts, &dm);

286:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
287:   DMMoabGetInterface(dm, &mbImpl);
288:   DMMoabGetLocalElements(dm, &elocal);

290:   /* zero out the discrete operator */
291:   MatZeroEntries(Jpre);

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

296:   const int& idl = dof_indices[0];
297:   const int& idr = dof_indices[1];
298:   const PetscScalar dxxL = user->alpha/hx,dxxR = -user->alpha/hx;
299:   const PetscScalar bcvals[2][2] = {{hx,0},{0,hx}};
300:   const PetscScalar e_vals[2][2][2] = {{{a *hx/2+dxxL,0},{dxxR,0}},
301:                                       {{0,a*hx/2+dxxL},{0,dxxR}}};

303:   /* Compute function over the locally owned part of the grid 
304:      Assemble the operator by looping over edges and computing
305:      contribution for each vertex dof                         */
306:   for(moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
307:     const moab::EntityHandle ehandle = *iter;

309:     // Get connectivity information in canonical order
310:     DMMoabGetElementConnectivity(dm, ehandle, &vpere, &connect);
311:     if (vpere != 2) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only EDGE2 element bases are supported in the current example. n(Connectivity)=%D.\n", vpere);

313:     DMMoabGetDofsBlocked(dm, vpere, connect, dof_indices);

315:     const PetscInt lcols[] = {idl,idr}, rcols[] = {idr, idl};

317:     /* check if element is on the boundary */
318:     DMMoabIsEntityOnBoundary(dm,ehandle,&elem_on_boundary);

320:     if (elem_on_boundary) {
321:       if (idl == 0) {
322:         // Left Boundary conditions...
323:         MatSetValuesBlocked(Jpre,1,&idl,1,&idl,&bcvals[0][0],ADD_VALUES);
324:         MatSetValuesBlocked(Jpre,1,&idr,2,rcols,&e_vals[0][0][0],ADD_VALUES);
325:       }
326:       else {
327:         // Right Boundary conditions...
328:         MatSetValuesBlocked(Jpre,1,&idr,1,&idr,&bcvals[0][0],ADD_VALUES);
329:         MatSetValuesBlocked(Jpre,1,&idl,2,lcols,&e_vals[0][0][0],ADD_VALUES);
330:       }
331:     }
332:     else {
333:       MatSetValuesBlocked(Jpre,1,&idr,2,rcols,&e_vals[0][0][0],ADD_VALUES);
334:       MatSetValuesBlocked(Jpre,1,&idl,2,lcols,&e_vals[0][0][0],ADD_VALUES);
335:     }
336:   }

338:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
339:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
340:   if (J != Jpre) {
341:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
342:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
343:   }
344:   return(0);
345: }


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:   PetscErrorCode    ierr;
357:   moab::Interface*  mbImpl;
358:   const moab::Range *vowned;
359:   PetscInt          dof_index;
360:   moab::Range::iterator iter;

363:   TSGetDM(ts, &dm);
364: 
365:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
366:   DMMoabGetInterface(dm, &mbImpl);

368:   DMMoabGetLocalVertices(dm, &vowned, NULL);

370:   VecSet(X, 0.0);

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

375:   /* Compute function over the locally owned part of the grid */
376:   for(moab::Range::iterator iter = vowned->begin(); iter != vowned->end(); iter++) {
377:     const moab::EntityHandle vhandle = *iter;
378:     DMMoabGetDofsBlockedLocal(dm, 1, &vhandle, &dof_index);

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

383:     PetscReal xi = vpos[0];
384:     x[dof_index].u = user->leftbc.u*(1.-xi) + user->rightbc.u*xi + sin(2.*PETSC_PI*xi);
385:     x[dof_index].v = user->leftbc.v*(1.-xi) + user->rightbc.v*xi;
386:   }

388:   /* Restore vectors */
389:   DMMoabVecRestoreArray(dm, X, &x);
390:   return(0);
391: }


396: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
397: {
398:   UserCtx         user = (UserCtx)ctx;
399:   DM              dm;
400:   Field           *x,*xdot,*f;
401:   PetscReal       hx,vpos[2*3];
402:   PetscErrorCode  ierr;
403:   PetscInt        dof_indices[2],bc_indices[2];
404:   const moab::EntityHandle *connect;
405:   PetscInt        vpere=2,nloc,ngh;
406:   PetscBool       elem_on_boundary;
407:   const int& idx_left = dof_indices[0];
408:   const int& idx_right = dof_indices[1];
409:   moab::Interface*  mbImpl;
410:   const moab::Range   *elocal;

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

416:   /* get the essential MOAB mesh related quantities needed for FEM assembly */
417:   DMMoabGetInterface(dm, &mbImpl);
418:   DMMoabGetLocalElements(dm, &elocal);
419:   DMMoabGetLocalSize(dm, NULL, NULL, &nloc, &ngh);

421:   /* reset the residual vector */
422:   VecSet(F,0.0);

424:   /* get the local representation of the arrays from Vectors */
425:   DMMoabVecGetArrayRead(dm, X, &x);
426:   DMMoabVecGetArrayRead(dm, Xdot, &xdot);
427:   DMMoabVecGetArray(dm, F, &f);

429:   /* loop over local elements */
430:   for(moab::Range::iterator iter = elocal->begin(); iter != elocal->end(); iter++) {
431:     const moab::EntityHandle ehandle = *iter;

433:     // Get connectivity information in canonical order
434:     DMMoabGetElementConnectivity(dm, ehandle, &vpere, &connect);
435:     if (vpere != 2) SETERRQ1(PETSC_COMM_WORLD, PETSC_ERR_ARG_WRONG, "Only EDGE2 element bases are supported in the current example. n(Connectivity)=%D.\n", vpere);

437:     /* compute the mid-point of the element and use a 1-point lumped quadrature */
438:     DMMoabGetVertexCoordinates(dm,vpere,connect,vpos);

440:     DMMoabGetDofsBlockedLocal(dm, vpere, connect, dof_indices);

442:     /* check if element is on the boundary */
443:     DMMoabIsEntityOnBoundary(dm,ehandle,&elem_on_boundary);

445:     if (elem_on_boundary) {
446:       DMMoabGetDofsBlocked(dm, vpere, connect, bc_indices);
447:       if (bc_indices[0] == 0) {      /* Apply left BC */
448:         f[idx_left].u = hx * (x[idx_left].u - user->leftbc.u);
449:         f[idx_left].v = hx * (x[idx_left].v - user->leftbc.v);
450:         f[idx_right].u += user->alpha*(x[idx_right].u-x[idx_left].u)/hx;
451:         f[idx_right].v += user->alpha*(x[idx_right].v-x[idx_left].v)/hx;
452:       }
453:       else {                        /* Apply right BC */
454:         f[idx_left].u += hx * xdot[idx_left].u + user->alpha*(x[idx_left].u - x[idx_right].u)/hx;
455:         f[idx_left].v += hx * xdot[idx_left].v + user->alpha*(x[idx_left].v - x[idx_right].v)/hx;
456:         f[idx_right].u = hx * (x[idx_right].u - user->rightbc.u);
457:         f[idx_right].v = hx * (x[idx_right].v - user->rightbc.v);
458:       }
459:     }
460:     else {
461:       f[idx_left].u += hx * xdot[idx_left].u + user->alpha*(x[idx_left].u - x[idx_right].u)/hx;
462:       f[idx_left].v += hx * xdot[idx_left].v + user->alpha*(x[idx_left].v - x[idx_right].v)/hx;
463:       f[idx_right].u += user->alpha*(x[idx_right].u-x[idx_left].u)/hx;
464:       f[idx_right].v += user->alpha*(x[idx_right].v-x[idx_left].v)/hx;
465:     }
466:   }

468:   /* Restore data */
469:   DMMoabVecRestoreArrayRead(dm, X, &x);
470:   DMMoabVecRestoreArrayRead(dm, Xdot, &xdot);
471:   DMMoabVecRestoreArray(dm, F, &f);
472:   return(0);
473: }