Actual source code: ex11.c

petsc-master 2016-06-28
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} AdvectSolType;
132: static const char *const AdvectSolTypes[] = {"TILTED","BUMP","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 Error;
154:   } functional;
155: } Physics_Advect;

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

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

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

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

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

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

205: static PetscErrorCode PhysicsSolution_Advect(Model mod,PetscReal time,const PetscReal *x,PetscScalar *u,void *ctx)
206: {
207:   Physics        phys    = (Physics)ctx;
208:   Physics_Advect *advect = (Physics_Advect*)phys->data;

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

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

252:   PhysicsSolution_Advect(mod,time,x,yexact,phys);
253:   f[advect->functional.Error] = PetscAbsScalar(y[0]-yexact[0]);
254:   return(0);
255: }

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

265:   /* Register "canned" boundary conditions and defaults for where to apply. */
266:   PetscDSAddBoundary(prob, PETSC_TRUE, "inflow",  "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Inflow,  ALEN(inflowids),  inflowids,  phys);
267:   PetscDSAddBoundary(prob, PETSC_FALSE, "outflow", "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_Advect_Outflow, ALEN(outflowids), outflowids, phys);
268:   return(0);
269: }

273: static PetscErrorCode PhysicsCreate_Advect(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)
274: {
275:   Physics_Advect *advect;

279:   phys->field_desc = PhysicsFields_Advect;
280:   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Advect;
281:   PetscNew(&advect);
282:   phys->data       = advect;
283:   mod->setupbc = SetUpBC_Advect;

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

324: /******************* Shallow Water ********************/
325: typedef struct {
326:   PetscReal gravity;
327:   PetscReal boundaryHeight;
328:   struct {
329:     PetscInt Height;
330:     PetscInt Speed;
331:     PetscInt Energy;
332:   } functional;
333: } Physics_SW;
334: typedef struct {
335:   PetscScalar vals[0];
336:   PetscScalar h;
337:   PetscScalar uh[DIM];
338: } SWNode;

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

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

356:   Scale2(1./x->h,x->uh,u);
357:   uhn  = Dot2(x->uh,n);
358:   f->h = uhn;
359:   for (i=0; i<DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
360:   return(0);
361: }

365: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
366: {
368:   xG[0] = xI[0];
369:   xG[1] = -xI[1];
370:   xG[2] = -xI[2];
371:   return(0);
372: }

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

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

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

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

416: static PetscErrorCode PhysicsFunctional_SW(Model mod,PetscReal time,const PetscReal *coord,const PetscScalar *xx,PetscReal *f,void *ctx)
417: {
418:   Physics      phys = (Physics)ctx;
419:   Physics_SW   *sw  = (Physics_SW*)phys->data;
420:   const SWNode *x   = (const SWNode*)xx;
421:   PetscScalar  u[2];
422:   PetscReal    h;

425:   h = PetscRealPart(x->h);
426:   Scale2(1./x->h,x->uh,u);
427:   f[sw->functional.Height] = h;
428:   f[sw->functional.Speed]  = Norm2(u) + PetscSqrtReal(sw->gravity*h);
429:   f[sw->functional.Energy] = 0.5*(Dot2(x->uh,u) + sw->gravity*PetscSqr(h));
430:   return(0);
431: }

435: static PetscErrorCode SetUpBC_SW(PetscDS prob,Physics phys)
436: {
438:   const PetscInt wallids[] = {100,101,200,300};
440:   PetscDSAddBoundary(prob, PETSC_FALSE, "wall", "Face Sets", 0, 0, NULL, (void (*)()) PhysicsBoundary_SW_Wall, ALEN(wallids), wallids, phys);
441:   return(0);
442: }

446: static PetscErrorCode PhysicsCreate_SW(Model mod,Physics phys,PetscOptionItems *PetscOptionsObject)
447: {
448:   Physics_SW     *sw;

452:   phys->field_desc = PhysicsFields_SW;
453:   phys->riema