Actual source code: ex11.c

petsc-3.7.3 2016-08-01

  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 <petscds.h>
 39: #include <petscts.h>
 40: #include <petscsf.h> /* For SplitFaces() */

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

 45: static PetscFunctionList PhysicsList;

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

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

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

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

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

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

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

104: struct _n_User {
105:   PetscInt numSplitFaces;
106:   PetscInt vtkInterval;   /* For monitor */
107:   Model    model;
108: };

110: PETSC_STATIC_INLINE PetscScalar DotDIM(const PetscScalar *x,const PetscScalar *y)
111: {
112:   PetscInt    i;
113:   PetscScalar prod=0.0;

115:   for (i=0; i<DIM; i++) prod += x[i]*y[i];
116:   return prod;
117: }
118: PETSC_STATIC_INLINE PetscReal NormDIM(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(DotDIM(x,x))); }

120: PETSC_STATIC_INLINE PetscScalar Dot2(const PetscScalar *x,const PetscScalar *y) { return x[0]*y[0] + x[1]*y[1];}
121: PETSC_STATIC_INLINE PetscReal Norm2(const PetscScalar *x) { return PetscSqrtReal(PetscAbsScalar(Dot2(x,x)));}
122: PETSC_STATIC_INLINE void Normalize2(PetscScalar *x) { PetscReal a = 1./Norm2(x); x[0] *= a; x[1] *= a; }
123: 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]; }
124: PETSC_STATIC_INLINE void Scale2(PetscScalar a,const PetscScalar *x,PetscScalar *y) { y[0] = a*x[0]; y[1] = a*x[1]; }

126: /******************* Advect ********************/
127: typedef enum {ADVECT_SOL_TILTED,ADVECT_SOL_BUMP} AdvectSolType;
128: static const char *const AdvectSolTypes[] = {"TILTED","BUMP","AdvectSolType","ADVECT_SOL_",0};
129: typedef enum {ADVECT_SOL_BUMP_CONE,ADVECT_SOL_BUMP_COS} AdvectSolBumpType;
130: static const char *const AdvectSolBumpTypes[] = {"CONE","COS","AdvectSolBumpType","ADVECT_SOL_BUMP_",0};

132: typedef struct {
133:   PetscReal wind[DIM];
134: } Physics_Advect_Tilted;
135: typedef struct {
136:   PetscReal         center[DIM];
137:   PetscReal         radius;
138:   AdvectSolBumpType type;
139: } Physics_Advect_Bump;

141: typedef struct {
142:   PetscReal     inflowState;
143:   AdvectSolType soltype;
144:   union {
145:     Physics_Advect_Tilted tilted;
146:     Physics_Advect_Bump   bump;
147:   } sol;
148:   struct {
149:     PetscInt Error;
150:   } functional;
151: } Physics_Advect;

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

157: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
158: {
159:   Physics        phys    = (Physics)ctx;
160:   Physics_Advect *advect = (Physics_Advect*)phys->data;

163:   xG[0] = advect->inflowState;
164:   return(0);
165: }

169: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
170: {
172:   xG[0] = xI[0];
173:   return(0);
174: }

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

183:   switch (advect->soltype) {
184:   case ADVECT_SOL_TILTED: {
185:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
186:     wind[0] = tilted->wind[0];
187:     wind[1] = tilted->wind[1];
188:   } break;
189:   case ADVECT_SOL_BUMP:
190:     wind[0] = -qp[1];
191:     wind[1] = qp[0];
192:     break;
193:     /* default: SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
194:   }
195:   wn      = Dot2(wind, n);
196:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
197: }

201: static PetscErrorCode PhysicsSolution_Advect(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
202: {
203:   Physics        phys    = (Physics)ctx;
204:   Physics_Advect *advect = (Physics_Advect*)phys->data;

207:   switch (advect->soltype) {
208:   case ADVECT_SOL_TILTED: {
209:     PetscReal             x0[DIM];
210:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
211:     Waxpy2(-time,tilted->wind,x,x0);
212:     if (x0[1] > 0) u[0] = 1.*x[0] + 3.*x[1];
213:     else u[0] = advect->inflowState;
214:   } break;
215:   case ADVECT_SOL_BUMP: {
216:     Physics_Advect_Bump *bump = &advect->sol.bump;
217:     PetscReal           x0[DIM],v[DIM],r,cost,sint;
218:     cost  = PetscCosReal(time);
219:     sint  = PetscSinReal(time);
220:     x0[0] = cost*x[0] + sint*x[1];
221:     x0[1] = -sint*x[0] + cost*x[1];
222:     Waxpy2(-1,bump->center,x0,v);
223:     r = Norm2(v);
224:     switch (bump->type) {
225:     case ADVECT_SOL_BUMP_CONE:
226:       u[0] = PetscMax(1 - r/bump->radius,0);
227:       break;
228:     case ADVECT_SOL_BUMP_COS:
229:       u[0] = 0.5 + 0.5*PetscCosReal(PetscMin(r/bump->radius,1)*PETSC_PI);
230:       break;
231:     }
232:   } break;
233:   default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Unknown solution type");
234:   }
235:   return(0);
236: }

240: static PetscErrorCode PhysicsFunctional_Advect(Model mod,PetscReal time,const PetscScalar *x,const PetscScalar *y,PetscReal *f,void *ctx)
241: {
242:   Physics        phys    = (Physics)ctx;
243:   Physics_Advect *advect = (Physics_Advect*)phys->data;
244:   PetscScalar    yexact[1];

248:   PhysicsSolution_Advect(mod,time,x,yexact,phys);
249:   f[advect->functional.Error] = PetscAbsScalar(y[0]-yexact[0]);
250:   return(0);
251: }

255: static PetscErrorCode SetUpBC_Advect(DM dm, Physics phys)
256: {
258:   const PetscInt inflowids[] = {100,200,300},outflowids[] = {101};
260:   /* Register "canned" boundary conditions and defaults for where to apply. */
261:   DMAddBoundary(dm, PETSC_TRUE, "inflow",  "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Inflow,  ALEN(inflowids),  inflowids,  phys);
262:   DMAddBoundary(dm, PETSC_FALSE, "outflow", "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Outflow, ALEN(outflowids), outflowids, phys);
263:   return(0);
264: }

268: static PetscErrorCode PhysicsCreate_Advect(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)
269: {
270:   Physics_Advect *advect;

274:   phys->field_desc = PhysicsFields_Advect;
275:   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Advect;
276:   PetscNew(&advect);
277:   phys->data       = advect;
278:   mod->setupbc = SetUpBC_Advect;

280:   PetscOptionsHead(PetscOptionsObject,"Advect options");
281:   {
282:     PetscInt two = 2,dof = 1;
283:     advect->soltype = ADVECT_SOL_TILTED;
284:     PetscOptionsEnum("-advect_sol_type","solution type","",AdvectSolTypes,(PetscEnum)advect->soltype,(PetscEnum*)&advect->soltype,NULL);
285:     switch (advect->soltype) {
286:     case ADVECT_SOL_TILTED: {
287:       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
288:       two = 2;
289:       tilted->wind[0] = 0.0;
290:       tilted->wind[1] = 1.0;
291:       PetscOptionsRealArray("-advect_tilted_wind","background wind vx,vy","",tilted->wind,&two,NULL);
292:       advect->inflowState = -2.0;
293:       PetscOptionsRealArray("-advect_tilted_inflow","Inflow state","",&advect->inflowState,&dof,NULL);
294:       phys->maxspeed = Norm2(tilted->wind);
295:     } break;
296:     case ADVECT_SOL_BUMP: {
297:       Physics_Advect_Bump *bump = &advect->sol.bump;
298:       two = 2;
299:       bump->center[0] = 2.;
300:       bump->center[1] = 0.;
301:       PetscOptionsRealArray("-advect_bump_center","location of center of bump x,y","",bump->center,&two,NULL);
302:       bump->radius = 0.9;
303:       PetscOptionsReal("-advect_bump_radius","radius of bump","",bump->radius,&bump->radius,NULL);
304:       bump->type = ADVECT_SOL_BUMP_CONE;
305:       PetscOptionsEnum("-advect_bump_type","type of bump","",AdvectSolBumpTypes,(PetscEnum)bump->type,(PetscEnum*)&bump->type,NULL);
306:       phys->maxspeed = 3.;       /* radius of mesh, kludge */
307:     } break;
308:     }
309:   }
310:   PetscOptionsTail();
311:   /* Initial/transient solution with default boundary conditions */
312:   ModelSolutionSetDefault(mod,PhysicsSolution_Advect,phys);
313:   /* Register "canned" functionals */
314:   ModelFunctionalRegister(mod,"Error",&advect->functional.Error,PhysicsFunctional_Advect,phys);
315:   mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_GHOSTED;
316:   return(0);
317: }

319: /******************* Shallow Water ********************/
320: typedef struct {
321:   PetscReal gravity;
322:   PetscReal boundaryHeight;
323:   struct {
324:     PetscInt Height;
325:     PetscInt Speed;
326:     PetscInt Energy;
327:   } functional;
328: } Physics_SW;
329: typedef struct {
330:   PetscScalar vals[0];
331:   PetscScalar h;
332:   PetscScalar uh[DIM];
333: } SWNode;

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

339: /*
340:  * h_t + div(uh) = 0
341:  * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
342:  *
343:  * */
344: static PetscErrorCode SWFlux(Physics phys,const PetscReal *n,const SWNode *x,SWNode *f)
345: {
346:   Physics_SW  *sw = (Physics_SW*)phys->data;
347:   PetscScalar uhn,u[DIM];
348:   PetscInt    i;

351:   Scale2(1./x->h,x->uh,u);
352:   uhn  = Dot2(x->uh,n);
353:   f->h = uhn;
354:   for (i=0; i<DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
355:   return(0);
356: }

360: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
361: {
363:   xG[0] = xI[0];
364:   xG[1] = -xI[1];
365:   xG[2] = -xI[2];
366:   return(0);
367: }

371: static void PhysicsRiemann_SW(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
372: {
373:   Physics_SW   *sw = (Physics_SW*)phys->data;
374:   PetscReal    cL,cR,speed,nn[DIM];
375:   const SWNode *uL = (const SWNode*)xL,*uR = (const SWNode*)xR;
376:   SWNode       fL,fR;
377:   PetscInt     i;

379:   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"); */
380:   nn[0] = n[0];
381:   nn[1] = n[1];
382:   Normalize2(nn);
383:   SWFlux(phys,nn,uL,&fL);
384:   SWFlux(phys,nn,uR,&fR);
385:   cL    = PetscSqrtReal(sw->gravity*PetscRealPart(uL->h));
386:   cR    = PetscSqrtReal(sw->gravity*PetscRealPart(uR->h)); /* gravity wave speed */
387:   speed = PetscMax(PetscAbsScalar(Dot2(uL->uh,nn)/uL->h) + cL,PetscAbsScalar(Dot2(uR->uh,nn)/uR->h) + cR);
388:   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);
389: }

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

398:   if (time != 0.0) SETERRQ1(mod->comm,PETSC_ERR_SUP,"No solution known for time %g",(double)time);
399:   dx[0] = x[0] - 1.5;
400:   dx[1] = x[1] - 1.0;
401:   r     = Norm2(dx);
402:   sigma = 0.5;
403:   u[0]  = 1 + 2*PetscExpScalar(-PetscSqr(r)/(2*PetscSqr(sigma)));
404:   u[1]  = 0.0;
405:   u[2]  = 0.0;
406:   return(0);
407: }

411: static PetscErrorCode PhysicsFunctional_SW(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)
412: {
413:   Physics      phys = (Physics)ctx;
414:   Physics_SW   *sw  = (Physics_SW*)phys->data;
415:   const SWNode *x   = (const SWNode*)xx;
416:   PetscScalar  u[2];
417:   PetscReal    h;

420:   h = PetscRealPart(x->h);
421:   Scale2(1./x->h,x->uh,u);
422:   f[sw->functional.Height] = h;
423:   f[sw->functional.Speed]  = Norm2(u) + PetscSqrtReal(sw->gravity*h);
424:   f[sw->functional.Energy] = 0.5*(Dot2(x->uh,u) + sw->gravity*PetscSqr(h));
425:   return(0);
426: }

430: static PetscErrorCode SetUpBC_SW(DM dm,Physics phys)
431: {
433:   const PetscInt wallids[] = {100,101,200,300};
435:   DMAddBoundary(dm, PETSC_FALSE, "wall", "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_SW_Wall, ALEN(wallids), wallids, phys);
436:   return(0);
437: }

441: static PetscErrorCode PhysicsCreate_SW(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)
442: {
443:   Physics_SW     *sw;

447:   phys->field_desc = PhysicsFields_SW;
448:   phys->riemann = (PetscRiemannFunc) PhysicsRiemann_SW;
449:   PetscNew(&sw);
450:   phys->data    = sw;
451:   mod->setupbc  = SetUpBC_SW;

453:   PetscOptionsHead(PetscOptionsObject,"SW options");
454:   {
455:     sw->gravity = 1.0;
456:     PetscOptionsReal("-sw_gravity","Gravitational constant","",sw->gravity,&sw->gravity,NULL);
457:   }
458:   PetscOptionsTail();
459:   phys->maxspeed = PetscSqrtReal(2.0*sw->gravity); /* Mach 1 for depth of 2 */

461:   ModelSolutionSetDefault(mod,PhysicsSolution_SW,phys);
462:   ModelFunctionalRegister(mod,"Height",&sw->functional.Height,PhysicsFunctional_SW,phys);
463:   ModelFunctionalRegister(mod,"Speed",&sw->functional.Speed,PhysicsFunctional_SW,phys);
464:   ModelFunctionalRegister(mod,"Energy",&sw->functional.Energy,PhysicsFunctional_SW,phys);

466:   mod->bcs[0] = mod->bcs[1] = mod->bcs[2] = DM_BOUNDARY_GHOSTED;

468:   return(0);
469: }

471: /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
Binary file (standard input) matches