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;
 40:   PetscErrorCode    ierr;

 43:   PetscNew(&user);

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

 72:   *puser = user;
 73:   return(0);
 74: }

 76: PetscErrorCode Destroy_AppContext(UserCtx *user)
 77: {

 81:   PetscFree(*user);
 82:   return(0);
 83: }

 85: static PetscErrorCode FormInitialSolution(TS,Vec,void*);
 86: static PetscErrorCode FormRHSFunction(TS,PetscReal,Vec,Vec,void*);
 87: static PetscErrorCode FormIFunction(TS,PetscReal,Vec,Vec,Vec,void*);
 88: static PetscErrorCode FormIJacobian(TS,PetscReal,Vec,Vec,PetscReal,Mat,Mat,void*);

 90: /****************
 91:  *              *
 92:  *     MAIN     *
 93:  *              *
 94:  ****************/
 95: int main(int argc,char **argv)
 96: {
 97:   TS                ts;         /* nonlinear solver */
 98:   Vec               X;          /* solution, residual vectors */
 99:   Mat               J;          /* Jacobian matrix */
100:   PetscInt          steps,mx;
101:   PetscErrorCode    ierr;
102:   PetscReal         hx,dt,ftime;
103:   UserCtx           user;       /* user-defined work context */
104:   TSConvergedReason reason;

106:   DM                dm;
107:   const char        *fields[2] = {"U","V"};

109:   PetscInitialize(&argc,&argv,(char *)0,help);if (ierr) return ierr;

111:   /* Initialize the user context struct */
112:   Initialize_AppContext(&user);

114:   /* Fill in the user defined work context: */
115:   DMMoabCreateBoxMesh(PETSC_COMM_WORLD, 1, PETSC_FALSE, NULL, user->n, 1, &dm);
116:   DMMoabSetFieldNames(dm, user->nvars, fields);
117:   DMMoabSetBlockSize(dm, user->nvars);
118:   DMSetFromOptions(dm);

120:   /* SetUp the data structures for DMMOAB */
121:   DMSetUp(dm);

123:   /*  Create timestepping solver context */
124:   TSCreate(PETSC_COMM_WORLD,&ts);
125:   TSSetDM(ts, dm);
126:   TSSetType(ts,TSARKIMEX);
127:   TSSetEquationType(ts,TS_EQ_DAE_IMPLICIT_INDEX1);
128:   DMSetMatType(dm,MATBAIJ);
129:   DMCreateMatrix(dm,&J);

131:   TSSetRHSFunction(ts,NULL,FormRHSFunction,user);
132:   TSSetIFunction(ts,NULL,FormIFunction,user);
133:   TSSetIJacobian(ts,J,J,FormIJacobian,user);

135:   ftime = 10.0;
136:   TSSetMaxSteps(ts,user->ntsteps);
137:   TSSetMaxTime(ts,ftime);
138:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

140:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141:      Create the solution vector and set the initial conditions
142:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143:   DMCreateGlobalVector(dm, &X);

145:   FormInitialSolution(ts,X,user);
146:   TSSetSolution(ts,X);
147:   VecGetSize(X,&mx);
148:   hx = 1.0/(PetscReal)(mx/2-1);
149:   dt = 0.4 * PetscSqr(hx) / user->alpha; /* Diffusive stability limit */
150:   TSSetTimeStep(ts,dt);

152:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153:      Set runtime options
154:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
155:   TSSetFromOptions(ts);

157:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158:      Solve nonlinear system
159:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160:   TSSolve(ts,X);
161:   TSGetSolveTime(ts,&ftime);
162:   TSGetStepNumber(ts,&steps);
163:   TSGetConvergedReason(ts,&reason);
164:   PetscPrintf(PETSC_COMM_WORLD,"%s at time %g after %D steps\n",TSConvergedReasons[reason],ftime,steps);

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

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

184:   /* Free work space.
185:      Free all PETSc related resources: */
186:   MatDestroy(&J);
187:   VecDestroy(&X);
188:   TSDestroy(&ts);
189:   DMDestroy(&dm);

191:   /* Free all MOAB related resources: */
192:   Destroy_AppContext(&user);
193:   PetscFinalize();
194:   return ierr;
195: }


198: /*
199:   IJacobian - Compute IJacobian = dF/dU + a dF/dUdot
200: */
201: PetscErrorCode FormIJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat J,Mat Jpre,void *ptr)
202: {
203:   UserCtx             user = (UserCtx)ptr;
204:   PetscErrorCode      ierr;
205:   PetscInt            dof;
206:   PetscReal           hx;
207:   DM                  dm;
208:   const moab::Range   *vlocal;
209:   PetscBool           vonboundary;

212:   TSGetDM(ts, &dm);

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

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

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

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

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

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

244:   MatAssemblyBegin(Jpre,MAT_FINAL_ASSEMBLY);
245:   MatAssemblyEnd(Jpre,MAT_FINAL_ASSEMBLY);
246:   if (J != Jpre) {
247:     MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
248:     MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
249:   }
250:   return(0);
251: }


254: static PetscErrorCode FormRHSFunction(TS ts,PetscReal t,Vec X,Vec F,void *ptr)
255: {
256:   UserCtx           user = (UserCtx)ptr;
257:   DM                dm;
258:   PetscReal         hx;
259:   const Field       *x;
260:   Field             *f;
261:   PetscInt          dof;
262:   const moab::Range *ownedvtx;
263:   PetscErrorCode    ierr;

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

269:   /* Get pointers to vector data */
270:   VecSet(F,0.0);

272:   DMMoabVecGetArrayRead(dm, X, &x);
273:   DMMoabVecGetArray(dm, F, &f);

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

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

282:     PetscScalar u = x[dof].u,v = x[dof].v;
283:     f[dof].u = hx*(user->A + u*u*v - (user->B+1)*u);
284:     f[dof].v = hx*(user->B*u - u*u*v);
285:   }

287:   /* Restore vectors */
288:   DMMoabVecRestoreArrayRead(dm, X, &x);
289:   DMMoabVecRestoreArray(dm, F, &f);
290:   return(0);
291: }


294: static PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,void *ctx)
295: {
296:   UserCtx         user = (UserCtx)ctx;
297:   DM              dm;
298:   Field           *x,*xdot,*f;
299:   PetscReal       hx;
300:   Vec             Xloc;
301:   PetscErrorCode  ierr;
302:   PetscInt        i,bcindx;
303:   PetscBool       elem_on_boundary;
304:   const moab::Range   *vlocal;

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

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

313:   /* reset the residual vector */
314:   VecSet(F,0.0);

316:   DMGetLocalVector(dm,&Xloc);
317:   DMGlobalToLocalBegin(dm,X,INSERT_VALUES,Xloc);
318:   DMGlobalToLocalEnd(dm,X,INSERT_VALUES,Xloc);

320:   /* get the local representation of the arrays from Vectors */
321:   DMMoabVecGetArrayRead(dm, Xloc, &x);
322:   DMMoabVecGetArrayRead(dm, Xdot, &xdot);
323:   DMMoabVecGetArray(dm, F, &f);

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

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

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

334:     if (elem_on_boundary) {
335:       DMMoabGetDofsBlocked(dm, 1, &vhandle, &bcindx);
336:       if (bcindx == 0) {  /* Apply left BC */
337:         f[i].u = hx * (x[i].u - user->leftbc.u);
338:         f[i].v = hx * (x[i].v - user->leftbc.v);
339:       } else {       /* Apply right BC */
340:         f[i].u = hx * (x[i].u - user->rightbc.u);
341:         f[i].v = hx * (x[i].v - user->rightbc.v);
342:       }
343:     }
344:     else {
345:       f[i].u = hx * xdot[i].u - user->alpha * (x[i-1].u - 2.*x[i].u + x[i+1].u) / hx;
346:       f[i].v = hx * xdot[i].v - user->alpha * (x[i-1].v - 2.*x[i].v + x[i+1].v) / hx;
347:     }
348:   }

350:   /* Restore data */
351:   DMMoabVecRestoreArrayRead(dm, Xloc, &x);
352:   DMMoabVecRestoreArrayRead(dm, Xdot, &xdot);
353:   DMMoabVecRestoreArray(dm, F, &f);
354:   DMRestoreLocalVector(dm, &Xloc);
355:   return(0);
356: }


359: PetscErrorCode FormInitialSolution(TS ts,Vec X,void *ctx)
360: {
361:   UserCtx           user = (UserCtx)ctx;
362:   PetscReal         vpos[3];
363:   DM                dm;
364:   Field             *x;
365:   PetscErrorCode    ierr;
366:   const moab::Range *vowned;
367:   PetscInt          dof;
368:   moab::Range::iterator iter;

371:   TSGetDM(ts, &dm);

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

376:   VecSet(X, 0.0);

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

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

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

389:     PetscReal xi = vpos[0];
390:     x[dof].u = user->leftbc.u*(1.-xi) + user->rightbc.u*xi + PetscSinReal(2.*PETSC_PI*xi);
391:     x[dof].v = user->leftbc.v*(1.-xi) + user->rightbc.v*xi;
392:   }

394:   /* Restore vectors */
395:   DMMoabVecRestoreArray(dm, X, &x);
396:   return(0);
397: }

399: /*TEST

401:     build:
402:       requires: moab

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

407:     test:
408:       suffix: 2
409:       nsize: 2
410:       args: -n 50 -ts_type glee -ts_adapt_type none -ts_dt 0.1 -io
411:       TODO:

413: TEST*/