Actual source code: ex11.c

petsc-master 2016-07-26
Report Typos and Errors
  1: static char help[] = "Second Order TVD Finite Volume Example.\n";

We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops
over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values,
\begin{equation}
f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
\end{equation}
and then update the cell values given the cell volume.
\begin{eqnarray}
f^L_i &-=& \frac{f_i}{vol^L} \\
f^R_i &+=& \frac{f_i}{vol^R}
\end{eqnarray}

As an example, we can consider the shallow water wave equation,
\begin{eqnarray}
h_t + \nabla\cdot \left( uh \right) &=& 0 \\
(uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0
\end{eqnarray}
where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity.

A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
\begin{eqnarray}
f^{L,R}_h &=& uh^{L,R} \cdot \hat n \\
f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
c^{L,R} &=& \sqrt{g h^{L,R}} \\
s &=& \max\left( \left|\frac{uh^L \cdot \hat n}{h^L}\right| + c^L, \left|\frac{uh^R \cdot \hat n}{h^R}\right| + c^R \right) \\
f_i &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
\end{eqnarray}
where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.

The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial
over a neighborhood of the given element.

The mesh is read in from an ExodusII file, usually generated by Cubit.
 37: #include <petscdmplex.h>
 38: #include <petscdmforest.h>
 39: #include <petscds.h>
 40: #include <petscts.h>
 41: #include <petscsf.h> /* For SplitFaces() */

 43: #define DIM 2                   /* Geometric dimension */
 44: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

 46: static PetscFunctionList PhysicsList;

 48: /* Represents continuum physical equations. */
 49: typedef struct _n_Physics *Physics;

 51: /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
 52:  * discretization-independent, but its members depend on the scenario being solved. */
 53: typedef struct _n_Model *Model;

 55: /* 'User' implements a discretization of a continuous model. */
 56: typedef struct _n_User *User;
 57: typedef PetscErrorCode (*SolutionFunction)(Model,PetscReal,const PetscReal*,PetscScalar*,void*);
 58: typedef PetscErrorCode (*SetUpBCFunction)(PetscDS,Physics);
 59: typedef PetscErrorCode (*FunctionalFunction)(Model,PetscReal,const PetscReal*,const PetscScalar*,PetscReal*,void*);
 60: typedef PetscErrorCode (*SetupFields)(Physics,PetscSection);
 61: static PetscErrorCode ModelSolutionSetDefault(Model,SolutionFunction,void*);
 62: static PetscErrorCode ModelFunctionalRegister(Model,const char*,PetscInt*,FunctionalFunction,void*);
 63: static PetscErrorCode OutputVTK(DM,const char*,PetscViewer*);

 65: struct FieldDescription {
 66:   const char *name;
 67:   PetscInt dof;
 68: };

 70: typedef struct _n_FunctionalLink *FunctionalLink;
 71: struct _n_FunctionalLink {
 72:   char               *name;
 73:   FunctionalFunction func;
 74:   void               *ctx;
 75:   PetscInt           offset;
 76:   FunctionalLink     next;
 77: };

 79: struct _n_Physics {
 80:   PetscRiemannFunc riemann;
 81:   PetscInt         dof;          /* number of degrees of freedom per cell */
 82:   PetscReal        maxspeed;     /* kludge to pick initial time step, need to add monitoring and step control */
 83:   void             *data;
 84:   PetscInt         nfields;
 85:   const struct FieldDescription *field_desc;
 86: };

 88: struct _n_Model {
 89:   MPI_Comm         comm;        /* Does not do collective communicaton, but some error conditions can be collective */
 90:   Physics          physics;
 91:   FunctionalLink   functionalRegistry;
 92:   PetscInt         maxComputed;
 93:   PetscInt         numMonitored;
 94:   FunctionalLink   *functionalMonitored;
 95:   PetscInt         numCall;
 96:   FunctionalLink   *functionalCall;
 97:   SolutionFunction solution;
 98:   SetUpBCFunction  setupbc;
 99:   void             *solutionctx;
100:   PetscReal        maxspeed;    /* estimate of global maximum speed (for CFL calculation) */
101:   PetscReal        bounds[2*DIM];
102:   DMBoundaryType   bcs[3];
103:   PetscErrorCode   (*errorIndicator)(PetscInt, PetscReal, PetscInt, const PetscScalar[], const PetscScalar[], PetscReal *, void *);
104:   void             *errorCtx;
105: };

107: struct _n_User {
108:   PetscInt numSplitFaces;
109:   PetscInt vtkInterval;   /* For monitor */
110:   PetscInt monitorStepOffset;
111:   Model    model;
112: };

114: PETSC_STATIC_INLINE PetscScalar DotDIM(const PetscScalar *x,const PetscScalar *y)
115: {
116:   PetscInt    i;
117:   PetscScalar prod=0.0;

119:   for (i=0; i<DIM; i++) prod += x[i]*y[i];
120:   return prod;
121: }
122: PETSC_STATIC_INLINE PetscReal NormDIM(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(DotDIM(x,x))); }

124: PETSC_STATIC_INLINE PetscScalar Dot2(const PetscScalar *x,const PetscScalar *y) { return x[0]*y[0] + x[1]*y[1];}
125: PETSC_STATIC_INLINE PetscReal Norm2(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(Dot2(x,x)));}
126: PETSC_STATIC_INLINE void Normalize2(PetscScalar *x) { PetscReal a = 1./Norm2(x); x[0] *= a; x[1] *= a; }
127: PETSC_STATIC_INLINE void Waxpy2(PetscScalar a,const PetscScalar *x,const PetscScalar *y,PetscScalar *w) { w[0] = a*x[0] + y[0]; w[1] = a*x[1] + y[1]; }
128: PETSC_STATIC_INLINE void Scale2(PetscScalar a,const PetscScalar *x,PetscScalar *y) { y[0] = a*x[0]; y[1] = a*x[1]; }

130: /******************* Advect ********************/
131: typedef enum {ADVECT_SOL_TILTED,ADVECT_SOL_BUMP,ADVECT_SOL_BUMP_CAVITY} AdvectSolType;
132: static const char *const AdvectSolTypes[] = {"TILTED","BUMP","BUMP_CAVITY","AdvectSolType","ADVECT_SOL_",0};
133: typedef enum {ADVECT_SOL_BUMP_CONE,ADVECT_SOL_BUMP_COS} AdvectSolBumpType;
134: static const char *const AdvectSolBumpTypes[] = {"CONE","COS","AdvectSolBumpType","ADVECT_SOL_BUMP_",0};

136: typedef struct {
137:   PetscReal wind[DIM];
138: } Physics_Advect_Tilted;
139: typedef struct {
140:   PetscReal         center[DIM];
141:   PetscReal         radius;
142:   AdvectSolBumpType type;
143: } Physics_Advect_Bump;

145: typedef struct {
146:   PetscReal     inflowState;
147:   AdvectSolType soltype;
148:   union {
149:     Physics_Advect_Tilted tilted;
150:     Physics_Advect_Bump   bump;
151:   } sol;
152:   struct {
153:     PetscInt Solution;
154:     PetscInt Error;
155:   } functional;
156: } Physics_Advect;

158: static const struct FieldDescription PhysicsFields_Advect[] = {{"U",1},{NULL,0}};

162: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
163: {
164:   Physics        phys    = (Physics)ctx;
165:   Physics_Advect *advect = (Physics_Advect*)phys->data;

168:   xG[0] = advect->inflowState;
169:   return(0);
170: }

174: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
175: {
177:   xG[0] = xI[0];
178:   return(0);
179: }

183: static void PhysicsRiemann_Advect(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
184: {
185:   Physics_Advect *advect = (Physics_Advect*)phys->data;
186:   PetscReal      wind[DIM],wn;

188:   switch (advect->soltype) {
189:   case ADVECT_SOL_TILTED: {
190:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
191:     wind[0] = tilted->wind[0];
192:     wind[1] = tilted->wind[1];
193:   } break;
194:   case ADVECT_SOL_BUMP:
195:     wind[0] = -qp[1];
196:     wind[1] = qp[0];
197:     break;
198:   case ADVECT_SOL_BUMP_CAVITY:
199:     {
200:       PetscInt  i;
201:       PetscReal comp2[3], rad2;

203:       rad2 = 0.;
204:       for (i = 0; i < dim; i++) {
205:         comp2[i] = qp[i] * qp[i];
206:         rad2    += comp2[i];
207:       }

209:       wind[0] = -qp[1];
210:       wind[1] = qp[0];
211:       if (rad2 > 1.) {
212:         PetscInt  maxI = 0;
213:         PetscReal maxComp2 = comp2[0];

215:         for (i = 1; i < dim; i++) {
216:           if (comp2[i] > maxComp2) {
217:             maxI     = i;
218:             maxComp2 = comp2[i];
219:           }
220:         }
221:         wind[maxI] = 0.;
222:       }
223:     }
224:     break;
225:     /* default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
226:   }
227:   wn      = Dot2(wind, n);
228:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
229: }

233: static PetscErrorCode PhysicsSolution_Advect(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
234: {
235:   Physics        phys    = (Physics)ctx;
236:   Physics_Advect *advect = (Physics_Advect*)phys->data;

239:   switch (advect->soltype) {
240:   case ADVECT_SOL_TILTED: {
241:     PetscReal             x0[DIM];
242:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
243:     Waxpy2(-time,tilted->wind,x,x0);
244:     if (x0[1] > 0) u[0] = 1.*x[0] + 3.*x[1];
245:     else u[0] = advect->inflowState;
246:   } break;
247:   case ADVECT_SOL_BUMP_CAVITY:
248:   case ADVECT_SOL_BUMP: {
249:     Physics_Advect_Bump *bump = &advect->sol.bump;
250:     PetscReal           x0[DIM],v[DIM],r,cost,sint;
251:     cost  = PetscCosReal(time);
252:     sint  = PetscSinReal(time);
253:     x0[0] = cost*x[0] + sint*x[1];
254:     x0[1] = -sint*x[0] + cost*x[1];
255:     Waxpy2(-1,bump->center,x0,v);
256:     r = Norm2(v);
257:     switch (bump->type) {
258:     case ADVECT_SOL_BUMP_CONE:
259:       u[0] = PetscMax(1 - r/bump->radius,0);
260:       break;
261:     case ADVECT_SOL_BUMP_COS:
262:       u[0] = 0.5 + 0.5*PetscCosReal(PetscMin(r/bump->radius,1)*PETSC_PI);
263:       break;
264:     }
265:   } break;
266:   default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown solution type");
267:   }
268:   return(0);
269: }

273: static PetscErrorCode PhysicsFunctional_Advect(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *y,PetscReal *f,void *ctx)
274: {
275:   Physics        phys    = (Physics)ctx;
276:   Physics_Advect *advect = (Physics_Advect*)phys->data;
277:   PetscScalar    yexact[1];

281:   PhysicsSolution_Advect(mod,time,x,yexact,phys);
282:   f[advect->functional.Solution] = PetscRealPart(y[0]);
283:   f[advect->functional.Error] = PetscAbsScalar(y[0]-yexact[0]);
284:   return(0);
285: }

289: static PetscErrorCode SetUpBC_Advect(PetscDS prob, Physics phys)
290: {
292:   const PetscInt inflowids[] = {100,200,300},outflowids[] = {101};

295:   /* Register "canned" boundary conditions and defaults for where to apply. */
296:   PetscDSAddBoundary(prob, PETSC_TRUE, "inflow",  "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Inflow,  ALEN(inflowids),   inflowids,  phys);
297:   PetscDSAddBoundary(prob, PETSC_FALSE, "outflow", "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Outflow, ALEN(outflowids), outflowids, phys);
298:   return(0);
299: }

303: static PetscErrorCode PhysicsCreate_Advect(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)
304: {
305:   Physics_Advect *advect;

309:   phys->field_desc = PhysicsFields_Advect;
310:   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Advect;
311:   PetscNew(&advect);
312:   phys->data       = advect;
313:   mod->setupbc = SetUpBC_Advect;

315:   PetscOptionsHead(PetscOptionsObject,"Advect options");
316:   {
317:     PetscInt two = 2,dof = 1;
318:     advect->soltype = ADVECT_SOL_TILTED;
319:     PetscOptionsEnum("-advect_sol_type","solution type","",AdvectSolTypes,(PetscEnum)advect->soltype,(PetscEnum*)&advect->soltype,NULL);
320:     switch (advect->soltype) {
321:     case ADVECT_SOL_TILTED: {
322:       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
323:       two = 2;
324:       tilted->wind[0] = 0.0;
325:       tilted->wind[1] = 1.0;
326:       PetscOptionsRealArray("-advect_tilted_wind","background wind vx,vy","",tilted->wind,&two,NULL);
327:       advect->inflowState = -2.0;
328:       PetscOptionsRealArray("-advect_tilted_inflow","Inflow state","",&advect->inflowState,&dof,NULL);
329:       phys->maxspeed = Norm2(tilted->wind);
330:     } break;
331:     case ADVECT_SOL_BUMP_CAVITY:
332:     case ADVECT_SOL_BUMP: {
333:       Physics_Advect_Bump *bump = &advect->sol.bump;
334:       two = 2;
335:       bump->center[0] = 2.;
336:       bump->center[1] = 0.;
337:       PetscOptionsRealArray("-advect_bump_center","location of center of bump x,y","",bump->center,&two,NULL);
338:       bump->radius = 0.9;
339:       PetscOptionsReal("-advect_bump_radius","radius of bump","",bump->radius,&bump->radius,NULL);
340:       bump->type = ADVECT_SOL_BUMP_CONE;
341:       PetscOptionsEnum("-advect_bump_type","type of bump","",AdvectSolBumpTypes,(PetscEnum)bump->type,(PetscEnum*)&bump->type,NULL);
342:       phys->maxspeed = 3.;       /* radius of mesh, kludge */
343:     } break;
344:     }
345:   }
346:   PetscOptionsTail();
347:   /* Initial/transient solution with default boundary conditions */
348:   ModelSolutionSetDefault(mod,PhysicsSolution_Advect,phys);
349:   /* Register "canned" functionals */
350:   ModelFunctionalRegister(mod,"Solution",&advect->functional.Solution,PhysicsFunctional_Advect,phys);
351:   ModelFunctionalRegister(mod,"Error",&advect->functional.Error,PhysicsFunctional_Advect,phys);
352:   mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_GHOSTED;
353:   return(0);
354: }

356: /******************* Shallow Water ********************/
357: typedef struct {
358:   PetscReal gravity;
359:   PetscReal boundaryHeight;
360:   struct {
361:     PetscInt Height;
362:     PetscInt Speed;
363:     PetscInt Energy;
364:   } functional;
365: } Physics_SW;
366: typedef struct {
367:   PetscScalar vals[0];
368:   PetscScalar h;
369:   PetscScalar uh[DIM];
370: } SWNode;

372: static const struct FieldDescription PhysicsFields_SW[] = {{"Height",1},{"Momentum",DIM},{NULL,0}};

376: /*
377:  * h_t + div(uh) = 0
378:  * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
379:  *
380:  * */
381: static PetscErrorCode SWFlux(Physics phys,const PetscReal *n,const SWNode *x,SWNode *f)
382: {
383:   Physics_SW  *sw = (Physics_SW*)phys->data;
384:   PetscScalar uhn,u[DIM];
385:   PetscInt    i;

388:   Scale2(1./x->h,x->uh,u);
389:   uhn  = Dot2(x->uh,n);
390:   f->h = uhn;
391:   for (i=0; i<DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
392:   return(0);
393: }

397: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
398: {
400:   xG[0] = xI[0];
401:   xG[1] = -xI[1];
402:   xG[2] = -xI[2];
403:   return(0);
404: }

408: static void PhysicsRiemann_SW(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
409: {
410:   Physics_SW   *sw = (Physics_SW*)phys->data;
411:   PetscReal    cL,cR,speed,nn[DIM];
412:   const SWNode *uL = (const SWNode*)xL,*uR = (const SWNode*)xR;
413:   SWNode       fL,fR;
414:   PetscInt     i;

416:   if (uL->h < 0 || uR->h < 0) {for (i=0; i<1+dim; i++) flux[i] = NAN; return;} /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); */
417:   nn[0] = n[0];
418:   nn[1] = n[1];
419:   Normalize2(nn);
420:   SWFlux(phys,nn,uL,&fL);
421:   SWFlux(phys,nn,uR,&fR);
422:   cL    = PetscSqrtReal(sw->gravity*PetscRealPart(uL->h));
423:   cR    = PetscSqrtReal(sw->gravity*PetscRealPart(uR->h)); /* gravity wave speed */
424:   speed = PetscMax(PetscAbsScalar(Dot2(uL->uh,nn)/uL->h) + cL,PetscAbsScalar(Dot2(uR->uh,nn)/uR->h) + cR);
425:   for (i=0; i<1+dim; i++) flux[i] = (0.5*(fL.vals[i] + fR.vals[i]) + 0.5*speed*(xL[i] - xR[i])) * Norm2(n);
426: }

430: static PetscErrorCode PhysicsSolution_SW(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
431: {
432:   PetscReal dx[2],r,sigma;

435:   if (time != 0.0) SETERRQ1(mod->comm,PETSC_ERR_SUP,"No solution known for time %g",(double)time);
436:   dx[0] = x[0] - 1.5;
437:   dx[1] = x[1] - 1.0;
438:   r     = Norm2(dx);
439:   sigma = 0.5;
440:   u[0]  = 1 + 2*PetscExpScalar(-PetscSqr(r)/(2*PetscSqr(sigma)));
441:   u[1]  = 0.0;
442:   u[2]  = 0.0;
443:   return(0);
444: }

448: static PetscErrorCode PhysicsFunctional_SW(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)
449: {
450:   Physics      phys = (Physics)ctx;
451:   Physics_SW   *sw  = (Physics_SW*)phys->data;
452:   const SWNode *x   = (const SWNode*)xx;
453:   PetscScalar  u[2];
454:   PetscReal    h;

457:   h = PetscRealPart(x->h);
458:   Scale2(1./x->h,x->uh,u);
459:   f[sw->functional.Height] = h;
460:   f[sw->functional.Speed]  = Norm2(u) + PetscSqrtReal(sw->gravity*h);
461:   f[sw->functional.Energy] = 0.5*(Dot2(x->uh,u) + sw->gravity*PetscSqr(h));
462:   return(0);
463: }

467: static PetscErrorCode SetUpBC_SW(PetscDS prob,Physics phys)
468: {
470:   const PetscInt wallids[] = {100,101,200,300};
472:   PetscDSAddBoundary(prob, PETSC_FALSE, "wall", "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_SW_Wall, ALEN(wallids), wallids, phys);
473:   return(0);
474: }

478: static PetscErrorCode PhysicsCreate_SW(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)
479: {
480:   Physics_SW     *sw;

484:   phys->field_desc = PhysicsFields_SW;
485:   phys->riemann = (PetscRiemannFunc) PhysicsRiemann_SW;
486:   PetscNew(&sw);
487:   phys->data    = sw;
488:   mod->setupbc  = SetUpBC_SW;

490:   PetscOptionsHead(PetscOptionsObject,"SW options");
491:   {
492:     sw->gravity = 1.0;
493:     PetscOptionsReal("-sw_gravity","Gravitational constant","",sw->gravity,&sw->gravity,NULL);
494:   }
495:   PetscOptionsTail();
496:   phys->maxspeed = PetscSqrtReal(2.0*sw->gravity); /* Mach 1 for depth of 2 */

498:   ModelSolutionSetDefault(mod,PhysicsSolution_SW,phys);
499:   ModelFunctionalRegister(mod,"Height",&sw->functional.Height,PhysicsFunctional_SW,phys);
500:   ModelFunctionalRegister(mod,"Speed",&sw->functional.Speed,PhysicsFunctional_SW,phys);
501:   ModelFunctionalRegister(mod,"Energy",&sw->functional.Energy,PhysicsFunctional_SW,phys);

503:   mod->bcs[0] = mod->bcs[1] = mod->bcs[2] =