Actual source code: ex11.c

  1: static char help[] = "Second Order TVD Finite Volume Example.\n";
  2: /*F

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

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

 22: A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
 23: \begin{eqnarray}
 24:   f^{L,R}_h    &=& uh^{L,R} \cdot \hat n \\
 25:   f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
 26:   c^{L,R}      &=& \sqrt{g h^{L,R}} \\
 27:   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) \\
 28:   f_i          &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
 29: \end{eqnarray}
 30: where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.

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

 35: The mesh is read in from an ExodusII file, usually generated by Cubit.
 36: F*/
 37: #include <petscdmplex.h>
 38: #include <petscdmforest.h>
 39: #include <petscds.h>
 40: #include <petscts.h>

 42: #define DIM 2 /* Geometric dimension */

 44: static PetscFunctionList PhysicsList, PhysicsRiemannList_SW;

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

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

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

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

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

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

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

104: struct _n_User {
105:   PetscInt  vtkInterval;                        /* For monitor */
106:   char      outputBasename[PETSC_MAX_PATH_LEN]; /* Basename for output files */
107:   PetscInt  monitorStepOffset;
108:   Model     model;
109:   PetscBool vtkmon;
110: };

112: static inline PetscReal DotDIMReal(const PetscReal *x, const PetscReal *y)
113: {
114:   PetscInt  i;
115:   PetscReal prod = 0.0;

117:   for (i = 0; i < DIM; i++) prod += x[i] * y[i];
118:   return prod;
119: }
120: static inline PetscReal NormDIM(const PetscReal *x)
121: {
122:   return PetscSqrtReal(PetscAbsReal(DotDIMReal(x, x)));
123: }

125: static inline PetscReal Dot2Real(const PetscReal *x, const PetscReal *y)
126: {
127:   return x[0] * y[0] + x[1] * y[1];
128: }
129: static inline PetscReal Norm2Real(const PetscReal *x)
130: {
131:   return PetscSqrtReal(PetscAbsReal(Dot2Real(x, x)));
132: }
133: static inline void Normalize2Real(PetscReal *x)
134: {
135:   PetscReal a = 1. / Norm2Real(x);
136:   x[0] *= a;
137:   x[1] *= a;
138: }
139: static inline void Waxpy2Real(PetscReal a, const PetscReal *x, const PetscReal *y, PetscReal *w)
140: {
141:   w[0] = a * x[0] + y[0];
142:   w[1] = a * x[1] + y[1];
143: }
144: static inline void Scale2Real(PetscReal a, const PetscReal *x, PetscReal *y)
145: {
146:   y[0] = a * x[0];
147:   y[1] = a * x[1];
148: }

150: /******************* Advect ********************/
151: typedef enum {
152:   ADVECT_SOL_TILTED,
153:   ADVECT_SOL_BUMP,
154:   ADVECT_SOL_BUMP_CAVITY
155: } AdvectSolType;
156: static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "BUMP_CAVITY", "AdvectSolType", "ADVECT_SOL_", 0};
157: typedef enum {
158:   ADVECT_SOL_BUMP_CONE,
159:   ADVECT_SOL_BUMP_COS
160: } AdvectSolBumpType;
161: static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};

163: typedef struct {
164:   PetscReal wind[DIM];
165: } Physics_Advect_Tilted;
166: typedef struct {
167:   PetscReal         center[DIM];
168:   PetscReal         radius;
169:   AdvectSolBumpType type;
170: } Physics_Advect_Bump;

172: typedef struct {
173:   PetscReal     inflowState;
174:   AdvectSolType soltype;
175:   union
176:   {
177:     Physics_Advect_Tilted tilted;
178:     Physics_Advect_Bump   bump;
179:   } sol;
180:   struct {
181:     PetscInt Solution;
182:     PetscInt Error;
183:   } functional;
184: } Physics_Advect;

186: static const struct FieldDescription PhysicsFields_Advect[] = {
187:   {"U",  1},
188:   {NULL, 0}
189: };

191: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
192: {
193:   Physics         phys   = (Physics)ctx;
194:   Physics_Advect *advect = (Physics_Advect *)phys->data;

197:   xG[0] = advect->inflowState;
198:   return 0;
199: }

201: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
202: {
204:   xG[0] = xI[0];
205:   return 0;
206: }

208: static void PhysicsRiemann_Advect(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
209: {
210:   Physics_Advect *advect = (Physics_Advect *)phys->data;
211:   PetscReal       wind[DIM], wn;

213:   switch (advect->soltype) {
214:   case ADVECT_SOL_TILTED: {
215:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
216:     wind[0]                       = tilted->wind[0];
217:     wind[1]                       = tilted->wind[1];
218:   } break;
219:   case ADVECT_SOL_BUMP:
220:     wind[0] = -qp[1];
221:     wind[1] = qp[0];
222:     break;
223:   case ADVECT_SOL_BUMP_CAVITY: {
224:     PetscInt  i;
225:     PetscReal comp2[3] = {0., 0., 0.}, rad2;

227:     rad2 = 0.;
228:     for (i = 0; i < dim; i++) {
229:       comp2[i] = qp[i] * qp[i];
230:       rad2 += comp2[i];
231:     }

233:     wind[0] = -qp[1];
234:     wind[1] = qp[0];
235:     if (rad2 > 1.) {
236:       PetscInt  maxI     = 0;
237:       PetscReal maxComp2 = comp2[0];

239:       for (i = 1; i < dim; i++) {
240:         if (comp2[i] > maxComp2) {
241:           maxI     = i;
242:           maxComp2 = comp2[i];
243:         }
244:       }
245:       wind[maxI] = 0.;
246:     }
247:   } break;
248:   default: {
249:     PetscInt i;
250:     for (i = 0; i < DIM; ++i) wind[i] = 0.0;
251:   }
252:     /* default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
253:   }
254:   wn      = Dot2Real(wind, n);
255:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
256: }

258: static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
259: {
260:   Physics         phys   = (Physics)ctx;
261:   Physics_Advect *advect = (Physics_Advect *)phys->data;

264:   switch (advect->soltype) {
265:   case ADVECT_SOL_TILTED: {
266:     PetscReal              x0[DIM];
267:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
268:     Waxpy2Real(-time, tilted->wind, x, x0);
269:     if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
270:     else u[0] = advect->inflowState;
271:   } break;
272:   case ADVECT_SOL_BUMP_CAVITY:
273:   case ADVECT_SOL_BUMP: {
274:     Physics_Advect_Bump *bump = &advect->sol.bump;
275:     PetscReal            x0[DIM], v[DIM], r, cost, sint;
276:     cost  = PetscCosReal(time);
277:     sint  = PetscSinReal(time);
278:     x0[0] = cost * x[0] + sint * x[1];
279:     x0[1] = -sint * x[0] + cost * x[1];
280:     Waxpy2Real(-1, bump->center, x0, v);
281:     r = Norm2Real(v);
282:     switch (bump->type) {
283:     case ADVECT_SOL_BUMP_CONE:
284:       u[0] = PetscMax(1 - r / bump->radius, 0);
285:       break;
286:     case ADVECT_SOL_BUMP_COS:
287:       u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
288:       break;
289:     }
290:   } break;
291:   default:
292:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
293:   }
294:   return 0;
295: }

297: static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx)
298: {
299:   Physics         phys      = (Physics)ctx;
300:   Physics_Advect *advect    = (Physics_Advect *)phys->data;
301:   PetscScalar     yexact[1] = {0.0};

304:   PhysicsSolution_Advect(mod, time, x, yexact, phys);
305:   f[advect->functional.Solution] = PetscRealPart(y[0]);
306:   f[advect->functional.Error]    = PetscAbsScalar(y[0] - yexact[0]);
307:   return 0;
308: }

310: static PetscErrorCode SetUpBC_Advect(DM dm, PetscDS prob, Physics phys)
311: {
312:   const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
313:   DMLabel        label;

316:   /* Register "canned" boundary conditions and defaults for where to apply. */
317:   DMGetLabel(dm, "Face Sets", &label);
318:   PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Inflow, NULL, phys, NULL);
319:   PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Outflow, NULL, phys, NULL);
320:   return 0;
321: }

323: static PetscErrorCode PhysicsCreate_Advect(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
324: {
325:   Physics_Advect *advect;

328:   phys->field_desc = PhysicsFields_Advect;
329:   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Advect;
330:   PetscNew(&advect);
331:   phys->data   = advect;
332:   mod->setupbc = SetUpBC_Advect;

334:   PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
335:   {
336:     PetscInt two = 2, dof = 1;
337:     advect->soltype = ADVECT_SOL_TILTED;
338:     PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL);
339:     switch (advect->soltype) {
340:     case ADVECT_SOL_TILTED: {
341:       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
342:       two                           = 2;
343:       tilted->wind[0]               = 0.0;
344:       tilted->wind[1]               = 1.0;
345:       PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL);
346:       advect->inflowState = -2.0;
347:       PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL);
348:       phys->maxspeed = Norm2Real(tilted->wind);
349:     } break;
350:     case ADVECT_SOL_BUMP_CAVITY:
351:     case ADVECT_SOL_BUMP: {
352:       Physics_Advect_Bump *bump = &advect->sol.bump;
353:       two                       = 2;
354:       bump->center[0]           = 2.;
355:       bump->center[1]           = 0.;
356:       PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL);
357:       bump->radius = 0.9;
358:       PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL);
359:       bump->type = ADVECT_SOL_BUMP_CONE;
360:       PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL);
361:       phys->maxspeed = 3.; /* radius of mesh, kludge */
362:     } break;
363:     }
364:   }
365:   PetscOptionsHeadEnd();
366:   /* Initial/transient solution with default boundary conditions */
367:   ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys);
368:   /* Register "canned" functionals */
369:   ModelFunctionalRegister(mod, "Solution", &advect->functional.Solution, PhysicsFunctional_Advect, phys);
370:   ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys);
371:   return 0;
372: }

374: /******************* Shallow Water ********************/
375: typedef struct {
376:   PetscReal gravity;
377:   PetscReal boundaryHeight;
378:   struct {
379:     PetscInt Height;
380:     PetscInt Speed;
381:     PetscInt Energy;
382:   } functional;
383: } Physics_SW;
384: typedef struct {
385:   PetscReal h;
386:   PetscReal uh[DIM];
387: } SWNode;
388: typedef union
389: {
390:   SWNode    swnode;
391:   PetscReal vals[DIM + 1];
392: } SWNodeUnion;

394: static const struct FieldDescription PhysicsFields_SW[] = {
395:   {"Height",   1  },
396:   {"Momentum", DIM},
397:   {NULL,       0  }
398: };

400: /*
401:  * h_t + div(uh) = 0
402:  * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
403:  *
404:  * */
405: static PetscErrorCode SWFlux(Physics phys, const PetscReal *n, const SWNode *x, SWNode *f)
406: {
407:   Physics_SW *sw = (Physics_SW *)phys->data;
408:   PetscReal   uhn, u[DIM];
409:   PetscInt    i;

412:   Scale2Real(1. / x->h, x->uh, u);
413:   uhn  = x->uh[0] * n[0] + x->uh[1] * n[1];
414:   f->h = uhn;
415:   for (i = 0; i < DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
416:   return 0;
417: }

419: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
420: {
422:   xG[0] = xI[0];
423:   xG[1] = -xI[1];
424:   xG[2] = -xI[2];
425:   return 0;
426: }

428: static void PhysicsRiemann_SW_HLL(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
429: {
430:   Physics_SW *sw = (Physics_SW *)phys->data;
431:   PetscReal   aL, aR;
432:   PetscReal   nn[DIM];
433: #if !defined(PETSC_USE_COMPLEX)
434:   const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
435: #else
436:   SWNodeUnion   uLreal, uRreal;
437:   const SWNode *uL = &uLreal.swnode;
438:   const SWNode *uR = &uRreal.swnode;
439: #endif
440:   SWNodeUnion fL, fR;
441:   PetscInt    i;
442:   PetscReal   zero = 0.;

444: #if defined(PETSC_USE_COMPLEX)
445:   uLreal.swnode.h = 0;
446:   uRreal.swnode.h = 0;
447:   for (i = 0; i < 1 + dim; i++) uLreal.vals[i] = PetscRealPart(xL[i]);
448:   for (i = 0; i < 1 + dim; i++) uRreal.vals[i] = PetscRealPart(xR[i]);
449: #endif
450:   if (uL->h <= 0 || uR->h <= 0) {
451:     for (i = 0; i < 1 + dim; i++) flux[i] = zero;
452:     return;
453:   } /* SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Reconstructed thickness is negative"); */
454:   nn[0] = n[0];
455:   nn[1] = n[1];
456:   Normalize2Real(nn);
457:   SWFlux(phys, nn, uL, &(fL.swnode));
458:   SWFlux(phys, nn, uR, &(fR.swnode));
459:   /* gravity wave speed */
460:   aL = PetscSqrtReal(sw->gravity * uL->h);
461:   aR = PetscSqrtReal(sw->gravity * uR->h);
462:   // Defining u_tilda and v_tilda as u and v
463:   PetscReal u_L, u_R;
464:   u_L = Dot2Real(uL->uh, nn) / uL->h;
465:   u_R = Dot2Real(uR->uh, nn) / uR->h;
466:   PetscReal sL, sR;
467:   sL = PetscMin(u_L - aL, u_R - aR);
468:   sR = PetscMax(u_L + aL, u_R + aR);
469:   if (sL > zero) {
470:     for (i = 0; i < dim + 1; i++) flux[i] = fL.vals[i] * Norm2Real(n);
471:   } else if (sR < zero) {
472:     for (i = 0; i < dim + 1; i++) flux[i] = fR.vals[i] * Norm2Real(n);
473:   } else {
474:     for (i = 0; i < dim + 1; i++) flux[i] = ((sR * fL.vals[i] - sL * fR.vals[i] + sR * sL * (xR[i] - xL[i])) / (sR - sL)) * Norm2Real(n);
475:   }
476: }

478: static void PhysicsRiemann_SW_Rusanov(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
479: {
480:   Physics_SW *sw = (Physics_SW *)phys->data;
481:   PetscReal   cL, cR, speed;
482:   PetscReal   nn[DIM];
483: #if !defined(PETSC_USE_COMPLEX)
484:   const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
485: #else
486:   SWNodeUnion   uLreal, uRreal;
487:   const SWNode *uL = &uLreal.swnode;
488:   const SWNode *uR = &uRreal.swnode;
489: #endif
490:   SWNodeUnion fL, fR;
491:   PetscInt    i;
492:   PetscReal   zero = 0.;

494: #if defined(PETSC_USE_COMPLEX)
495:   uLreal.swnode.h = 0;
496:   uRreal.swnode.h = 0;
497:   for (i = 0; i < 1 + dim; i++) uLreal.vals[i] = PetscRealPart(xL[i]);
498:   for (i = 0; i < 1 + dim; i++) uRreal.vals[i] = PetscRealPart(xR[i]);
499: #endif
500:   if (uL->h < 0 || uR->h < 0) {
501:     for (i = 0; i < 1 + dim; i++) flux[i] = zero / zero;
502:     return;
503:   } /* reconstructed thickness is negative */
504:   nn[0] = n[0];
505:   nn[1] = n[1];
506:   Normalize2Real(nn);
507:   SWFlux(phys, nn, uL, &(fL.swnode));
508:   SWFlux(phys, nn, uR, &(fR.swnode));
509:   cL    = PetscSqrtReal(sw->gravity * uL->h);
510:   cR    = PetscSqrtReal(sw->gravity * uR->h); /* gravity wave speed */
511:   speed = PetscMax(PetscAbsReal(Dot2Real(uL->uh, nn) / uL->h) + cL, PetscAbsReal(Dot2Real(uR->uh, nn) / uR->h) + cR);
512:   for (i = 0; i < 1 + dim; i++) flux[i] = (0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i])) * Norm2Real(n);
513: }

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

521:   dx[0] = x[0] - 1.5;
522:   dx[1] = x[1] - 1.0;
523:   r     = Norm2Real(dx);
524:   sigma = 0.5;
525:   u[0]  = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma)));
526:   u[1]  = 0.0;
527:   u[2]  = 0.0;
528:   return 0;
529: }

531: static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
532: {
533:   Physics       phys = (Physics)ctx;
534:   Physics_SW   *sw   = (Physics_SW *)phys->data;
535:   const SWNode *x    = (const SWNode *)xx;
536:   PetscReal     u[2];
537:   PetscReal     h;

540:   h = x->h;
541:   Scale2Real(1. / x->h, x->uh, u);
542:   f[sw->functional.Height] = h;
543:   f[sw->functional.Speed]  = Norm2Real(u) + PetscSqrtReal(sw->gravity * h);
544:   f[sw->functional.Energy] = 0.5 * (Dot2Real(x->uh, u) + sw->gravity * PetscSqr(h));
545:   return 0;
546: }

548: static PetscErrorCode SetUpBC_SW(DM dm, PetscDS prob, Physics phys)
549: {
550:   const PetscInt wallids[] = {100, 101, 200, 300};
551:   DMLabel        label;

554:   DMGetLabel(dm, "Face Sets", &label);
555:   PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_SW_Wall, NULL, phys, NULL);
556:   return 0;
557: }

559: static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
560: {
561:   Physics_SW *sw;
562:   char        sw_riemann[64] = "rusanov";

565:   phys->field_desc = PhysicsFields_SW;
566:   PetscNew(&sw);
567:   phys->data   = sw;
568:   mod->setupbc = SetUpBC_SW;

570:   PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov);
571:   PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL);

573:   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
574:   {
575:     void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
576:     sw->gravity = 1.0;
577:     PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL);
578:     PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL);
579:     PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW);
580:     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_SW;
581:   }
582:   PetscOptionsHeadEnd();
583:   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */

585:   ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys);
586:   ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys);
587:   ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys);
588:   ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys);

590:   return 0;
591: }

593: /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
594: /* An initial-value and self-similar solutions of the compressible Euler equations */
595: /* Ravi Samtaney and D. I. Pullin */
596: /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
597: typedef enum {
598:   EULER_PAR_GAMMA,
599:   EULER_PAR_RHOR,
600:   EULER_PAR_AMACH,
601:   EULER_PAR_ITANA,
602:   EULER_PAR_SIZE
603: } EulerParamIdx;
604: typedef enum {
605:   EULER_IV_SHOCK,
606:   EULER_SS_SHOCK,
607:   EULER_SHOCK_TUBE,
608:   EULER_LINEAR_WAVE
609: } EulerType;
610: typedef struct {
611:   PetscReal r;
612:   PetscReal ru[DIM];
613:   PetscReal E;
614: } EulerNode;
615: typedef union
616: {
617:   EulerNode eulernode;
618:   PetscReal vals[DIM + 2];
619: } EulerNodeUnion;
620: typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode *, PetscReal *);
621: typedef struct {
622:   EulerType       type;
623:   PetscReal       pars[EULER_PAR_SIZE];
624:   EquationOfState sound;
625:   struct {
626:     PetscInt Density;
627:     PetscInt Momentum;
628:     PetscInt Energy;
629:     PetscInt Pressure;
630:     PetscInt Speed;
631:   } monitor;
632: } Physics_Euler;

634: static const struct FieldDescription PhysicsFields_Euler[] = {
635:   {"Density",  1  },
636:   {"Momentum", DIM},
637:   {"Energy",   1  },
638:   {NULL,       0  }
639: };

641: /* initial condition */
642: int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
643: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
644: {
645:   PetscInt       i;
646:   Physics        phys = (Physics)ctx;
647:   Physics_Euler *eu   = (Physics_Euler *)phys->data;
648:   EulerNode     *uu   = (EulerNode *)u;
649:   PetscReal      p0, gamma, c;

653:   for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
654:   /* set E and rho */
655:   gamma = eu->pars[EULER_PAR_GAMMA];

657:   if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) {
658:     /******************* Euler Density Shock ********************/
659:     /* On initial-value and self-similar solutions of the compressible Euler equations */
660:     /* Ravi Samtaney and D. I. Pullin */
661:     /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
662:     /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity,  */
663:     p0 = 1.;
664:     if (x[0] < 0.0 + x[1] * eu->pars[EULER_PAR_ITANA]) {
665:       if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */
666:         PetscReal amach, rho, press, gas1, p1;
667:         amach     = eu->pars[EULER_PAR_AMACH];
668:         rho       = 1.;
669:         press     = p0;
670:         p1        = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0));
671:         gas1      = (gamma - 1.0) / (gamma + 1.0);
672:         uu->r     = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0);
673:         uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach);
674:         uu->E     = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0];
675:       } else {      /* left of discontinuity (0) */
676:         uu->r = 1.; /* rho = 1 */
677:         uu->E = p0 / (gamma - 1.0);
678:       }
679:     } else { /* right of discontinuity (2) */
680:       uu->r = eu->pars[EULER_PAR_RHOR];
681:       uu->E = p0 / (gamma - 1.0);
682:     }
683:   } else if (eu->type == EULER_SHOCK_TUBE) {
684:     /* For (x<x0) set (rho,u,p)=(8,0,10) and for (x>x0) set (rho,u,p)=(1,0,1). Choose x0 to the midpoint of the domain in the x-direction. */
685:     if (x[0] < 0.0) {
686:       uu->r = 8.;
687:       uu->E = 10. / (gamma - 1.);
688:     } else {
689:       uu->r = 1.;
690:       uu->E = 1. / (gamma - 1.);
691:     }
692:   } else if (eu->type == EULER_LINEAR_WAVE) {
693:     initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
694:   } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type);

696:   /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */
697:   eu->sound(&gamma, uu, &c);
698:   c = (uu->ru[0] / uu->r) + c;
699:   if (c > phys->maxspeed) phys->maxspeed = c;

701:   return 0;
702: }

704: static PetscErrorCode Pressure_PG(const PetscReal gamma, const EulerNode *x, PetscReal *p)
705: {
706:   PetscReal ru2;

709:   ru2  = DotDIMReal(x->ru, x->ru);
710:   (*p) = (x->E - 0.5 * ru2 / x->r) * (gamma - 1.0); /* (E - rho V^2/2)(gamma-1) = e rho (gamma-1) */
711:   return 0;
712: }

714: static PetscErrorCode SpeedOfSound_PG(const PetscReal *gamma, const EulerNode *x, PetscReal *c)
715: {
716:   PetscReal p;

719:   Pressure_PG(*gamma, x, &p);
721:   /* pars[EULER_PAR_GAMMA] = heat capacity ratio */
722:   (*c) = PetscSqrtReal(*gamma * p / x->r);
723:   return 0;
724: }

726: /*
727:  * x = (rho,rho*(u_1),...,rho*e)^T
728:  * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
729:  *
730:  * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
731:  *
732:  */
733: static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
734: {
735:   Physics_Euler *eu = (Physics_Euler *)phys->data;
736:   PetscReal      nu, p;
737:   PetscInt       i;

740:   Pressure_PG(eu->pars[EULER_PAR_GAMMA], x, &p);
741:   nu   = DotDIMReal(x->ru, n);
742:   f->r = nu;                                                     /* A rho u */
743:   nu /= x->r;                                                    /* A u */
744:   for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p; /* r u^2 + p */
745:   f->E = nu * (x->E + p);                                        /* u(e+p) */
746:   return 0;
747: }

749: /* PetscReal* => EulerNode* conversion */
750: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
751: {
752:   PetscInt         i;
753:   const EulerNode *xI   = (const EulerNode *)a_xI;
754:   EulerNode       *xG   = (EulerNode *)a_xG;
755:   Physics          phys = (Physics)ctx;
756:   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
758:   xG->r = xI->r;                                     /* ghost cell density - same */
759:   xG->E = xI->E;                                     /* ghost cell energy - same */
760:   if (n[1] != 0.) {                                  /* top and bottom */
761:     xG->ru[0] = xI->ru[0];                           /* copy tang to wall */
762:     xG->ru[1] = -xI->ru[1];                          /* reflect perp to t/b wall */
763:   } else {                                           /* sides */
764:     for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
765:   }
766:   if (eu->type == EULER_LINEAR_WAVE) { /* debug */
767: #if 0
768:     PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]);
769: #endif
770:   }
771:   return 0;
772: }
773: int godunovflux(const PetscScalar *ul, const PetscScalar *ur, PetscScalar *flux, const PetscReal *nn, const int *ndim, const PetscReal *gamma);
774: /* PetscReal* => EulerNode* conversion */
775: static void PhysicsRiemann_Euler_Godunov(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
776: {
777:   Physics_Euler *eu = (Physics_Euler *)phys->data;
778:   PetscReal      cL, cR, speed, velL, velR, nn[DIM], s2;
779:   PetscInt       i;

783:   for (i = 0, s2 = 0.; i < DIM; i++) {
784:     nn[i] = n[i];
785:     s2 += nn[i] * nn[i];
786:   }
787:   s2 = PetscSqrtReal(s2); /* |n|_2 = sum(n^2)^1/2 */
788:   for (i = 0.; i < DIM; i++) nn[i] /= s2;
789:   if (0) { /* Rusanov */
790:     const EulerNode *uL = (const EulerNode *)xL, *uR = (const EulerNode *)xR;
791:     EulerNodeUnion   fL, fR;
792:     EulerFlux(phys, nn, uL, &(fL.eulernode));
793:     EulerFlux(phys, nn, uR, &(fR.eulernode));
794:     eu->sound(&eu->pars[EULER_PAR_GAMMA], uL, &cL);
795:     if (ierr) exit(13);
796:     eu->sound(&eu->pars[EULER_PAR_GAMMA], uR, &cR);
797:     if (ierr) exit(14);
798:     velL  = DotDIMReal(uL->ru, nn) / uL->r;
799:     velR  = DotDIMReal(uR->ru, nn) / uR->r;
800:     speed = PetscMax(velR + cR, velL + cL);
801:     for (i = 0; i < 2 + dim; i++) flux[i] = 0.5 * ((fL.vals[i] + fR.vals[i]) + speed * (xL[i] - xR[i])) * s2;
802:   } else {
803:     int dim = DIM;
804:     /* int iwave =  */
805:     godunovflux(xL, xR, flux, nn, &dim, &eu->pars[EULER_PAR_GAMMA]);
806:     for (i = 0; i < 2 + dim; i++) flux[i] *= s2;
807:   }
808:   return;
809: }

811: static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
812: {
813:   Physics          phys = (Physics)ctx;
814:   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
815:   const EulerNode *x    = (const EulerNode *)xx;
816:   PetscReal        p;

819:   f[eu->monitor.Density]  = x->r;
820:   f[eu->monitor.Momentum] = NormDIM(x->ru);
821:   f[eu->monitor.Energy]   = x->E;
822:   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
823:   Pressure_PG(eu->pars[EULER_PAR_GAMMA], x, &p);
824:   f[eu->monitor.Pressure] = p;
825:   return 0;
826: }

828: static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys)
829: {
830:   Physics_Euler *eu = (Physics_Euler *)phys->data;
831:   DMLabel        label;

834:   DMGetLabel(dm, "Face Sets", &label);
835:   if (eu->type == EULER_LINEAR_WAVE) {
836:     const PetscInt wallids[] = {100, 101};
837:     PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL);
838:   } else {
839:     const PetscInt wallids[] = {100, 101, 200, 300};
840:     PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL);
841:   }
842:   return 0;
843: }

845: static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
846: {
847:   Physics_Euler *eu;

850:   phys->field_desc = PhysicsFields_Euler;
851:   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Euler_Godunov;
852:   PetscNew(&eu);
853:   phys->data   = eu;
854:   mod->setupbc = SetUpBC_Euler;
855:   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
856:   {
857:     PetscReal alpha;
858:     char      type[64] = "linear_wave";
859:     PetscBool is;
860:     eu->pars[EULER_PAR_GAMMA] = 1.4;
861:     eu->pars[EULER_PAR_AMACH] = 2.02;
862:     eu->pars[EULER_PAR_RHOR]  = 3.0;
863:     eu->pars[EULER_PAR_ITANA] = 0.57735026918963; /* angle of Euler self similar (SS) shock */
864:     PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->pars[EULER_PAR_GAMMA], &eu->pars[EULER_PAR_GAMMA], NULL);
865:     PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->pars[EULER_PAR_AMACH], &eu->pars[EULER_PAR_AMACH], NULL);
866:     PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->pars[EULER_PAR_RHOR], &eu->pars[EULER_PAR_RHOR], NULL);
867:     alpha = 60.;
868:     PetscOptionsReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL);
870:     eu->pars[EULER_PAR_ITANA] = 1. / PetscTanReal(alpha * PETSC_PI / 180.0);
871:     PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL);
872:     PetscStrcmp(type, "linear_wave", &is);
873:     if (is) {
874:       /* Remember this should be periodic */
875:       eu->type = EULER_LINEAR_WAVE;
876:       PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave");
877:     } else {
879:       PetscStrcmp(type, "iv_shock", &is);
880:       if (is) {
881:         eu->type = EULER_IV_SHOCK;
882:         PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock");
883:       } else {
884:         PetscStrcmp(type, "ss_shock", &is);
885:         if (is) {
886:           eu->type = EULER_SS_SHOCK;
887:           PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock");
888:         } else {
889:           PetscStrcmp(type, "shock_tube", &is);
890:           if (is) eu->type = EULER_SHOCK_TUBE;
891:           else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type);
892:           PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube");
893:         }
894:       }
895:     }
896:   }
897:   PetscOptionsHeadEnd();
898:   eu->sound      = SpeedOfSound_PG;
899:   phys->maxspeed = 0.; /* will get set in solution */
900:   ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys);
901:   ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys);
902:   ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys);
903:   ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys);
904:   ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys);
905:   ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys);

907:   return 0;
908: }

910: static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx)
911: {
912:   PetscReal err = 0.;
913:   PetscInt  i, j;

916:   for (i = 0; i < numComps; i++) {
917:     for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
918:   }
919:   *error = volume * err;
920:   return 0;
921: }

923: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
924: {
925:   PetscSF      sfPoint;
926:   PetscSection coordSection;
927:   Vec          coordinates;
928:   PetscSection sectionCell;
929:   PetscScalar *part;
930:   PetscInt     cStart, cEnd, c;
931:   PetscMPIInt  rank;

934:   DMGetCoordinateSection(dm, &coordSection);
935:   DMGetCoordinatesLocal(dm, &coordinates);
936:   DMClone(dm, dmCell);
937:   DMGetPointSF(dm, &sfPoint);
938:   DMSetPointSF(*dmCell, sfPoint);
939:   DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection);
940:   DMSetCoordinatesLocal(*dmCell, coordinates);
941:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
942:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);
943:   DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd);
944:   PetscSectionSetChart(sectionCell, cStart, cEnd);
945:   for (c = cStart; c < cEnd; ++c) PetscSectionSetDof(sectionCell, c, 1);
946:   PetscSectionSetUp(sectionCell);
947:   DMSetLocalSection(*dmCell, sectionCell);
948:   PetscSectionDestroy(&sectionCell);
949:   DMCreateLocalVector(*dmCell, partition);
950:   PetscObjectSetName((PetscObject)*partition, "partition");
951:   VecGetArray(*partition, &part);
952:   for (c = cStart; c < cEnd; ++c) {
953:     PetscScalar *p;

955:     DMPlexPointLocalRef(*dmCell, c, part, &p);
956:     p[0] = rank;
957:   }
958:   VecRestoreArray(*partition, &part);
959:   return 0;
960: }

962: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
963: {
964:   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
965:   PetscSection       coordSection;
966:   Vec                coordinates, facegeom, cellgeom;
967:   PetscSection       sectionMass;
968:   PetscScalar       *m;
969:   const PetscScalar *fgeom, *cgeom, *coords;
970:   PetscInt           vStart, vEnd, v;

973:   DMConvert(dm, DMPLEX, &plex);
974:   DMGetCoordinateSection(dm, &coordSection);
975:   DMGetCoordinatesLocal(dm, &coordinates);
976:   DMClone(dm, &dmMass);
977:   DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection);
978:   DMSetCoordinatesLocal(dmMass, coordinates);
979:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass);
980:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
981:   PetscSectionSetChart(sectionMass, vStart, vEnd);
982:   for (v = vStart; v < vEnd; ++v) {
983:     PetscInt numFaces;

985:     DMPlexGetSupportSize(dmMass, v, &numFaces);
986:     PetscSectionSetDof(sectionMass, v, numFaces * numFaces);
987:   }
988:   PetscSectionSetUp(sectionMass);
989:   DMSetLocalSection(dmMass, sectionMass);
990:   PetscSectionDestroy(&sectionMass);
991:   DMGetLocalVector(dmMass, massMatrix);
992:   VecGetArray(*massMatrix, &m);
993:   DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL);
994:   VecGetDM(facegeom, &dmFace);
995:   VecGetArrayRead(facegeom, &fgeom);
996:   VecGetDM(cellgeom, &dmCell);
997:   VecGetArrayRead(cellgeom, &cgeom);
998:   DMGetCoordinateDM(dm, &dmCoord);
999:   VecGetArrayRead(coordinates, &coords);
1000:   for (v = vStart; v < vEnd; ++v) {
1001:     const PetscInt  *faces;
1002:     PetscFVFaceGeom *fgA, *fgB, *cg;
1003:     PetscScalar     *vertex;
1004:     PetscInt         numFaces, sides[2], f, g;

1006:     DMPlexPointLocalRead(dmCoord, v, coords, &vertex);
1007:     DMPlexGetSupportSize(dmMass, v, &numFaces);
1008:     DMPlexGetSupport(dmMass, v, &faces);
1009:     for (f = 0; f < numFaces; ++f) {
1010:       sides[0] = faces[f];
1011:       DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA);
1012:       for (g = 0; g < numFaces; ++g) {
1013:         const PetscInt *cells = NULL;
1014:         PetscReal       area  = 0.0;
1015:         PetscInt        numCells;

1017:         sides[1] = faces[g];
1018:         DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB);
1019:         DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells);
1021:         DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg);
1022:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
1023:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
1024:         m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
1025:         DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells);
1026:       }
1027:     }
1028:   }
1029:   VecRestoreArrayRead(facegeom, &fgeom);
1030:   VecRestoreArrayRead(cellgeom, &cgeom);
1031:   VecRestoreArrayRead(coordinates, &coords);
1032:   VecRestoreArray(*massMatrix, &m);
1033:   DMDestroy(&dmMass);
1034:   DMDestroy(&plex);
1035:   return 0;
1036: }

1038: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
1039: static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
1040: {
1042:   mod->solution    = func;
1043:   mod->solutionctx = ctx;
1044:   return 0;
1045: }

1047: static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
1048: {
1049:   FunctionalLink link, *ptr;
1050:   PetscInt       lastoffset = -1;

1053:   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1054:   PetscNew(&link);
1055:   PetscStrallocpy(name, &link->name);
1056:   link->offset = lastoffset + 1;
1057:   link->func   = func;
1058:   link->ctx    = ctx;
1059:   link->next   = NULL;
1060:   *ptr         = link;
1061:   *offset      = link->offset;
1062:   return 0;
1063: }

1065: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject)
1066: {
1067:   PetscInt       i, j;
1068:   FunctionalLink link;
1069:   char          *names[256];

1072:   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
1073:   PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL);
1074:   /* Create list of functionals that will be computed somehow */
1075:   PetscMalloc1(mod->numMonitored, &mod->functionalMonitored);
1076:   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1077:   PetscMalloc1(mod->numMonitored, &mod->functionalCall);
1078:   mod->numCall = 0;
1079:   for (i = 0; i < mod->numMonitored; i++) {
1080:     for (link = mod->functionalRegistry; link; link = link->next) {
1081:       PetscBool match;
1082:       PetscStrcasecmp(names[i], link->name, &match);
1083:       if (match) break;
1084:     }
1086:     mod->functionalMonitored[i] = link;
1087:     for (j = 0; j < i; j++) {
1088:       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
1089:     }
1090:     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1091:   next_name:
1092:     PetscFree(names[i]);
1093:   }

1095:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1096:   mod->maxComputed = -1;
1097:   for (link = mod->functionalRegistry; link; link = link->next) {
1098:     for (i = 0; i < mod->numCall; i++) {
1099:       FunctionalLink call = mod->functionalCall[i];
1100:       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
1101:     }
1102:   }
1103:   return 0;
1104: }

1106: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
1107: {
1108:   FunctionalLink l, next;

1111:   if (!link) return 0;
1112:   l     = *link;
1113:   *link = NULL;
1114:   for (; l; l = next) {
1115:     next = l->next;
1116:     PetscFree(l->name);
1117:     PetscFree(l);
1118:   }
1119:   return 0;
1120: }

1122: /* put the solution callback into a functional callback */
1123: static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
1124: {
1125:   Model mod;
1127:   mod = (Model)modctx;
1128:   (*mod->solution)(mod, time, x, u, mod->solutionctx);
1129:   return 0;
1130: }

1132: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
1133: {
1134:   PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
1135:   void *ctx[1];
1136:   Model mod = user->model;

1139:   func[0] = SolutionFunctional;
1140:   ctx[0]  = (void *)mod;
1141:   DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X);
1142:   return 0;
1143: }

1145: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1146: {
1148:   PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);
1149:   PetscViewerSetType(*viewer, PETSCVIEWERVTK);
1150:   PetscViewerFileSetName(*viewer, filename);
1151:   return 0;
1152: }

1154: static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
1155: {
1156:   User        user = (User)ctx;
1157:   DM          dm, plex;
1158:   PetscViewer viewer;
1159:   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
1160:   PetscReal   xnorm;

1163:   PetscObjectSetName((PetscObject)X, "u");
1164:   VecGetDM(X, &dm);
1165:   VecNorm(X, NORM_INFINITY, &xnorm);

1167:   if (stepnum >= 0) stepnum += user->monitorStepOffset;
1168:   if (stepnum >= 0) { /* No summary for final time */
1169:     Model              mod = user->model;
1170:     Vec                cellgeom;
1171:     PetscInt           c, cStart, cEnd, fcount, i;
1172:     size_t             ftableused, ftablealloc;
1173:     const PetscScalar *cgeom, *x;
1174:     DM                 dmCell;
1175:     DMLabel            vtkLabel;
1176:     PetscReal         *fmin, *fmax, *fintegral, *ftmp;

1178:     DMConvert(dm, DMPLEX, &plex);
1179:     DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL);
1180:     fcount = mod->maxComputed + 1;
1181:     PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp);
1182:     for (i = 0; i < fcount; i++) {
1183:       fmin[i]      = PETSC_MAX_REAL;
1184:       fmax[i]      = PETSC_MIN_REAL;
1185:       fintegral[i] = 0;
1186:     }
1187:     VecGetDM(cellgeom, &dmCell);
1188:     DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd);
1189:     VecGetArrayRead(cellgeom, &cgeom);
1190:     VecGetArrayRead(X, &x);
1191:     DMGetLabel(dm, "vtk", &vtkLabel);
1192:     for (c = cStart; c < cEnd; ++c) {
1193:       PetscFVCellGeom   *cg;
1194:       const PetscScalar *cx     = NULL;
1195:       PetscInt           vtkVal = 0;

1197:       /* not that these two routines as currently implemented work for any dm with a
1198:        * localSection/globalSection */
1199:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1200:       DMPlexPointGlobalRead(dm, c, x, &cx);
1201:       if (vtkLabel) DMLabelGetValue(vtkLabel, c, &vtkVal);
1202:       if (!vtkVal || !cx) continue; /* ghost, or not a global cell */
1203:       for (i = 0; i < mod->numCall; i++) {
1204:         FunctionalLink flink = mod->functionalCall[i];
1205:         (*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx);
1206:       }
1207:       for (i = 0; i < fcount; i++) {
1208:         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1209:         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1210:         fintegral[i] += cg->volume * ftmp[i];
1211:       }
1212:     }
1213:     VecRestoreArrayRead(cellgeom, &cgeom);
1214:     VecRestoreArrayRead(X, &x);
1215:     DMDestroy(&plex);
1216:     MPI_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts));
1217:     MPI_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
1218:     MPI_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));

1220:     ftablealloc = fcount * 100;
1221:     ftableused  = 0;
1222:     PetscMalloc1(ftablealloc, &ftable);
1223:     for (i = 0; i < mod->numMonitored; i++) {
1224:       size_t         countused;
1225:       char           buffer[256], *p;
1226:       FunctionalLink flink = mod->functionalMonitored[i];
1227:       PetscInt       id    = flink->offset;
1228:       if (i % 3) {
1229:         PetscArraycpy(buffer, "  ", 2);
1230:         p = buffer + 2;
1231:       } else if (i) {
1232:         char newline[] = "\n";
1233:         PetscMemcpy(buffer, newline, sizeof(newline) - 1);
1234:         p = buffer + sizeof(newline) - 1;
1235:       } else {
1236:         p = buffer;
1237:       }
1238:       PetscSNPrintfCount(p, sizeof buffer - (p - buffer), "%12s [%10.7g,%10.7g] int %10.7g", &countused, flink->name, (double)fmin[id], (double)fmax[id], (double)fintegral[id]);
1239:       countused--;
1240:       countused += p - buffer;
1241:       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1242:         char *ftablenew;
1243:         ftablealloc = 2 * ftablealloc + countused;
1244:         PetscMalloc(ftablealloc, &ftablenew);
1245:         PetscArraycpy(ftablenew, ftable, ftableused);
1246:         PetscFree(ftable);
1247:         ftable = ftablenew;
1248:       }
1249:       PetscArraycpy(ftable + ftableused, buffer, countused);
1250:       ftableused += countused;
1251:       ftable[ftableused] = 0;
1252:     }
1253:     PetscFree4(fmin, fmax, fintegral, ftmp);

1255:     PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : "");
1256:     PetscFree(ftable);
1257:   }
1258:   if (user->vtkInterval < 1) return 0;
1259:   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1260:     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
1261:       TSGetStepNumber(ts, &stepnum);
1262:     }
1263:     PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum);
1264:     OutputVTK(dm, filename, &viewer);
1265:     VecView(X, viewer);
1266:     PetscViewerDestroy(&viewer);
1267:   }
1268:   return 0;
1269: }

1271: static PetscErrorCode initializeTS(DM dm, User user, TS *ts)
1272: {
1274:   TSCreate(PetscObjectComm((PetscObject)dm), ts);
1275:   TSSetType(*ts, TSSSP);
1276:   TSSetDM(*ts, dm);
1277:   if (user->vtkmon) TSMonitorSet(*ts, MonitorVTK, user, NULL);
1278:   DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user);
1279:   DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user);
1280:   TSSetMaxTime(*ts, 2.0);
1281:   TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER);
1282:   return 0;
1283: }

1285: static PetscErrorCode adaptToleranceFVM(PetscFV fvm, TS ts, Vec sol, VecTagger refineTag, VecTagger coarsenTag, User user, TS *tsNew, Vec *solNew)
1286: {
1287:   DM                 dm, gradDM, plex, cellDM, adaptedDM = NULL;
1288:   Vec                cellGeom, faceGeom;
1289:   PetscBool          isForest, computeGradient;
1290:   Vec                grad, locGrad, locX, errVec;
1291:   PetscInt           cStart, cEnd, c, dim, nRefine, nCoarsen;
1292:   PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2], minInd, maxInd, time;
1293:   PetscScalar       *errArray;
1294:   const PetscScalar *pointVals;
1295:   const PetscScalar *pointGrads;
1296:   const PetscScalar *pointGeom;
1297:   DMLabel            adaptLabel = NULL;
1298:   IS                 refineIS, coarsenIS;

1301:   TSGetTime(ts, &time);
1302:   VecGetDM(sol, &dm);
1303:   DMGetDimension(dm, &dim);
1304:   PetscFVGetComputeGradients(fvm, &computeGradient);
1305:   PetscFVSetComputeGradients(fvm, PETSC_TRUE);
1306:   DMIsForest(dm, &isForest);
1307:   DMConvert(dm, DMPLEX, &plex);
1308:   DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM);
1309:   DMCreateLocalVector(plex, &locX);
1310:   DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL);
1311:   DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX);
1312:   DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX);
1313:   DMCreateGlobalVector(gradDM, &grad);
1314:   DMPlexReconstructGradientsFVM(plex, locX, grad);
1315:   DMCreateLocalVector(gradDM, &locGrad);
1316:   DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad);
1317:   DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad);
1318:   VecDestroy(&grad);
1319:   DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd);
1320:   VecGetArrayRead(locGrad, &pointGrads);
1321:   VecGetArrayRead(cellGeom, &pointGeom);
1322:   VecGetArrayRead(locX, &pointVals);
1323:   VecGetDM(cellGeom, &cellDM);
1324:   DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel);
1325:   VecCreateMPI(PetscObjectComm((PetscObject)plex), cEnd - cStart, PETSC_DETERMINE, &errVec);
1326:   VecSetUp(errVec);
1327:   VecGetArray(errVec, &errArray);
1328:   for (c = cStart; c < cEnd; c++) {
1329:     PetscReal        errInd = 0.;
1330:     PetscScalar     *pointGrad;
1331:     PetscScalar     *pointVal;
1332:     PetscFVCellGeom *cg;

1334:     DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad);
1335:     DMPlexPointLocalRead(cellDM, c, pointGeom, &cg);
1336:     DMPlexPointLocalRead(plex, c, pointVals, &pointVal);

1338:     (user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx);
1339:     errArray[c - cStart] = errInd;
1340:     minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
1341:     minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
1342:   }
1343:   VecRestoreArray(errVec, &errArray);
1344:   VecRestoreArrayRead(locX, &pointVals);
1345:   VecRestoreArrayRead(cellGeom, &pointGeom);
1346:   VecRestoreArrayRead(locGrad, &pointGrads);
1347:   VecDestroy(&locGrad);
1348:   VecDestroy(&locX);
1349:   DMDestroy(&plex);

1351:   VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL);
1352:   VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL);
1353:   ISGetSize(refineIS, &nRefine);
1354:   ISGetSize(coarsenIS, &nCoarsen);
1355:   if (nRefine) DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS);
1356:   if (nCoarsen) DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS);
1357:   ISDestroy(&coarsenIS);
1358:   ISDestroy(&refineIS);
1359:   VecDestroy(&errVec);

1361:   PetscFVSetComputeGradients(fvm, computeGradient);
1362:   minMaxInd[1] = -minMaxInd[1];
1363:   MPI_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm));
1364:   minInd = minMaxIndGlobal[0];
1365:   maxInd = -minMaxIndGlobal[1];
1366:   PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minInd, (double)maxInd);
1367:   if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
1368:     DMAdaptLabel(dm, adaptLabel, &adaptedDM);
1369:   }
1370:   DMLabelDestroy(&adaptLabel);
1371:   if (adaptedDM) {
1372:     PetscInfo(ts, "Adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nRefine, nCoarsen);
1373:     if (tsNew) initializeTS(adaptedDM, user, tsNew);
1374:     if (solNew) {
1375:       DMCreateGlobalVector(adaptedDM, solNew);
1376:       PetscObjectSetName((PetscObject)*solNew, "solution");
1377:       DMForestTransferVec(dm, sol, adaptedDM, *solNew, PETSC_TRUE, time);
1378:     }
1379:     if (isForest) DMForestSetAdaptivityForest(adaptedDM, NULL); /* clear internal references to the previous dm */
1380:     DMDestroy(&adaptedDM);
1381:   } else {
1382:     if (tsNew) *tsNew = NULL;
1383:     if (solNew) *solNew = NULL;
1384:   }
1385:   return 0;
1386: }

1388: int main(int argc, char **argv)
1389: {
1390:   MPI_Comm          comm;
1391:   PetscDS           prob;
1392:   PetscFV           fvm;
1393:   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1394:   User              user;
1395:   Model             mod;
1396:   Physics           phys;
1397:   DM                dm, plex;
1398:   PetscReal         ftime, cfl, dt, minRadius;
1399:   PetscInt          dim, nsteps;
1400:   TS                ts;
1401:   TSConvergedReason reason;
1402:   Vec               X;
1403:   PetscViewer       viewer;
1404:   PetscBool         vtkCellGeom, useAMR;
1405:   PetscInt          adaptInterval;
1406:   char              physname[256] = "advect";
1407:   VecTagger         refineTag = NULL, coarsenTag = NULL;

1410:   PetscInitialize(&argc, &argv, (char *)0, help);
1411:   comm = PETSC_COMM_WORLD;

1413:   PetscNew(&user);
1414:   PetscNew(&user->model);
1415:   PetscNew(&user->model->physics);
1416:   mod           = user->model;
1417:   phys          = mod->physics;
1418:   mod->comm     = comm;
1419:   useAMR        = PETSC_FALSE;
1420:   adaptInterval = 1;

1422:   /* Register physical models to be available on the command line */
1423:   PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect);
1424:   PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW);
1425:   PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler);

1427:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1428:   {
1429:     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1430:     PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL);
1431:     user->vtkInterval = 1;
1432:     PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL);
1433:     user->vtkmon = PETSC_TRUE;
1434:     PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL);
1435:     vtkCellGeom = PETSC_FALSE;
1436:     PetscStrcpy(user->outputBasename, "ex11");
1437:     PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL);
1438:     PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL);
1439:     PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL);
1440:     PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL);
1441:   }
1442:   PetscOptionsEnd();

1444:   if (useAMR) {
1445:     VecTaggerBox refineBox, coarsenBox;

1447:     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1448:     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;

1450:     VecTaggerCreate(comm, &refineTag);
1451:     PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_");
1452:     VecTaggerSetType(refineTag, VECTAGGERABSOLUTE);
1453:     VecTaggerAbsoluteSetBox(refineTag, &refineBox);
1454:     VecTaggerSetFromOptions(refineTag);
1455:     VecTaggerSetUp(refineTag);
1456:     PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view");

1458:     VecTaggerCreate(comm, &coarsenTag);
1459:     PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_");
1460:     VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE);
1461:     VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox);
1462:     VecTaggerSetFromOptions(coarsenTag);
1463:     VecTaggerSetUp(coarsenTag);
1464:     PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view");
1465:   }

1467:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1468:   {
1469:     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
1470:     PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL);
1471:     PetscFunctionListFind(PhysicsList, physname, &physcreate);
1472:     PetscMemzero(phys, sizeof(struct _n_Physics));
1473:     (*physcreate)(mod, phys, PetscOptionsObject);
1474:     /* Count number of fields and dofs */
1475:     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1477:     ModelFunctionalSetFromOptions(mod, PetscOptionsObject);
1478:   }
1479:   PetscOptionsEnd();

1481:   /* Create mesh */
1482:   {
1483:     PetscInt i;

1485:     DMCreate(comm, &dm);
1486:     DMSetType(dm, DMPLEX);
1487:     DMSetFromOptions(dm);
1488:     for (i = 0; i < DIM; i++) {
1489:       mod->bounds[2 * i]     = 0.;
1490:       mod->bounds[2 * i + 1] = 1.;
1491:     };
1492:     dim = DIM;
1493:     { /* a null name means just do a hex box */
1494:       PetscInt  cells[3] = {1, 1, 1}, n = 3;
1495:       PetscBool flg2, skew              = PETSC_FALSE;
1496:       PetscInt  nret2 = 2 * DIM;
1497:       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
1498:       PetscOptionsRealArray("-grid_bounds", "bounds of the mesh in each direction (i.e., x_min,x_max,y_min,y_max", "", mod->bounds, &nret2, &flg2);
1499:       PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL);
1500:       PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL);
1501:       PetscOptionsEnd();
1502:       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1503:       if (flg2) {
1504:         PetscInt     dimEmbed, i;
1505:         PetscInt     nCoords;
1506:         PetscScalar *coords;
1507:         Vec          coordinates;

1509:         DMGetCoordinatesLocal(dm, &coordinates);
1510:         DMGetCoordinateDim(dm, &dimEmbed);
1511:         VecGetLocalSize(coordinates, &nCoords);
1513:         VecGetArray(coordinates, &coords);
1514:         for (i = 0; i < nCoords; i += dimEmbed) {
1515:           PetscInt j;

1517:           PetscScalar *coord = &coords[i];
1518:           for (j = 0; j < dimEmbed; j++) {
1519:             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1520:             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1521:               if (cells[0] == 2 && i == 8) {
1522:                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
1523:               } else if (cells[0] == 3) {
1524:                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1525:                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1526:                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1527:               }
1528:             }
1529:           }
1530:         }
1531:         VecRestoreArray(coordinates, &coords);
1532:         DMSetCoordinatesLocal(dm, coordinates);
1533:       }
1534:     }
1535:   }
1536:   DMViewFromOptions(dm, NULL, "-orig_dm_view");
1537:   DMGetDimension(dm, &dim);

1539:   /* set up BCs, functions, tags */
1540:   DMCreateLabel(dm, "Face Sets");
1541:   mod->errorIndicator = ErrorIndicator_Simple;

1543:   {
1544:     DM gdm;

1546:     DMPlexConstructGhostCells(dm, NULL, NULL, &gdm);
1547:     DMDestroy(&dm);
1548:     dm = gdm;
1549:     DMViewFromOptions(dm, NULL, "-dm_view");
1550:   }

1552:   PetscFVCreate(comm, &fvm);
1553:   PetscFVSetFromOptions(fvm);
1554:   PetscFVSetNumComponents(fvm, phys->dof);
1555:   PetscFVSetSpatialDimension(fvm, dim);
1556:   PetscObjectSetName((PetscObject)fvm, "");
1557:   {
1558:     PetscInt f, dof;
1559:     for (f = 0, dof = 0; f < phys->nfields; f++) {
1560:       PetscInt newDof = phys->field_desc[f].dof;

1562:       if (newDof == 1) {
1563:         PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name);
1564:       } else {
1565:         PetscInt j;

1567:         for (j = 0; j < newDof; j++) {
1568:           char compName[256] = "Unknown";

1570:           PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j);
1571:           PetscFVSetComponentName(fvm, dof + j, compName);
1572:         }
1573:       }
1574:       dof += newDof;
1575:     }
1576:   }
1577:   /* FV is now structured with one field having all physics as components */
1578:   DMAddField(dm, NULL, (PetscObject)fvm);
1579:   DMCreateDS(dm);
1580:   DMGetDS(dm, &prob);
1581:   PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann);
1582:   PetscDSSetContext(prob, 0, user->model->physics);
1583:   (*mod->setupbc)(dm, prob, phys);
1584:   PetscDSSetFromOptions(prob);
1585:   {
1586:     char      convType[256];
1587:     PetscBool flg;

1589:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
1590:     PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg);
1591:     PetscOptionsEnd();
1592:     if (flg) {
1593:       DM dmConv;

1595:       DMConvert(dm, convType, &dmConv);
1596:       if (dmConv) {
1597:         DMViewFromOptions(dmConv, NULL, "-dm_conv_view");
1598:         DMDestroy(&dm);
1599:         dm = dmConv;
1600:         DMSetFromOptions(dm);
1601:       }
1602:     }
1603:   }

1605:   initializeTS(dm, user, &ts);

1607:   DMCreateGlobalVector(dm, &X);
1608:   PetscObjectSetName((PetscObject)X, "solution");
1609:   SetInitialCondition(dm, X, user);
1610:   if (useAMR) {
1611:     PetscInt adaptIter;

1613:     /* use no limiting when reconstructing gradients for adaptivity */
1614:     PetscFVGetLimiter(fvm, &limiter);
1615:     PetscObjectReference((PetscObject)limiter);
1616:     PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter);
1617:     PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE);

1619:     PetscFVSetLimiter(fvm, noneLimiter);
1620:     for (adaptIter = 0;; ++adaptIter) {
1621:       PetscLogDouble bytes;
1622:       TS             tsNew = NULL;

1624:       PetscMemoryGetCurrentUsage(&bytes);
1625:       PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes);
1626:       DMViewFromOptions(dm, NULL, "-initial_dm_view");
1627:       VecViewFromOptions(X, NULL, "-initial_vec_view");
1628: #if 0
1629:       if (viewInitial) {
1630:         PetscViewer viewer;
1631:         char        buf[256];
1632:         PetscBool   isHDF5, isVTK;

1634:         PetscViewerCreate(comm,&viewer);
1635:         PetscViewerSetType(viewer,PETSCVIEWERVTK);
1636:         PetscViewerSetOptionsPrefix(viewer,"initial_");
1637:         PetscViewerSetFromOptions(viewer);
1638:         PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&isHDF5);
1639:         PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERVTK,&isVTK);
1640:         if (isHDF5) {
1641:           PetscSNPrintf(buf, 256, "ex11-initial-%" PetscInt_FMT ".h5", adaptIter);
1642:         } else if (isVTK) {
1643:           PetscSNPrintf(buf, 256, "ex11-initial-%" PetscInt_FMT ".vtu", adaptIter);
1644:           PetscViewerPushFormat(viewer,PETSC_VIEWER_VTK_VTU);
1645:         }
1646:         PetscViewerFileSetMode(viewer,FILE_MODE_WRITE);
1647:         PetscViewerFileSetName(viewer,buf);
1648:         if (isHDF5) {
1649:           DMView(dm,viewer);
1650:           PetscViewerFileSetMode(viewer,FILE_MODE_UPDATE);
1651:         }
1652:         VecView(X,viewer);
1653:         PetscViewerDestroy(&viewer);
1654:       }
1655: #endif

1657:       adaptToleranceFVM(fvm, ts, X, refineTag, coarsenTag, user, &tsNew, NULL);
1658:       if (!tsNew) {
1659:         break;
1660:       } else {
1661:         DMDestroy(&dm);
1662:         VecDestroy(&X);
1663:         TSDestroy(&ts);
1664:         ts = tsNew;
1665:         TSGetDM(ts, &dm);
1666:         PetscObjectReference((PetscObject)dm);
1667:         DMCreateGlobalVector(dm, &X);
1668:         PetscObjectSetName((PetscObject)X, "solution");
1669:         SetInitialCondition(dm, X, user);
1670:       }
1671:     }
1672:     /* restore original limiter */
1673:     PetscFVSetLimiter(fvm, limiter);
1674:   }

1676:   DMConvert(dm, DMPLEX, &plex);
1677:   if (vtkCellGeom) {
1678:     DM  dmCell;
1679:     Vec cellgeom, partition;

1681:     DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL);
1682:     OutputVTK(dm, "ex11-cellgeom.vtk", &viewer);
1683:     VecView(cellgeom, viewer);
1684:     PetscViewerDestroy(&viewer);
1685:     CreatePartitionVec(dm, &dmCell, &partition);
1686:     OutputVTK(dmCell, "ex11-partition.vtk", &viewer);
1687:     VecView(partition, viewer);
1688:     PetscViewerDestroy(&viewer);
1689:     VecDestroy(&partition);
1690:     DMDestroy(&dmCell);
1691:   }
1692:   /* collect max maxspeed from all processes -- todo */
1693:   DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius);
1694:   DMDestroy(&plex);
1695:   MPI_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
1697:   dt = cfl * minRadius / mod->maxspeed;
1698:   TSSetTimeStep(ts, dt);
1699:   TSSetFromOptions(ts);
1700:   if (!useAMR) {
1701:     TSSolve(ts, X);
1702:     TSGetSolveTime(ts, &ftime);
1703:     TSGetStepNumber(ts, &nsteps);
1704:   } else {
1705:     PetscReal finalTime;
1706:     PetscInt  adaptIter;
1707:     TS        tsNew  = NULL;
1708:     Vec       solNew = NULL;

1710:     TSGetMaxTime(ts, &finalTime);
1711:     TSSetMaxSteps(ts, adaptInterval);
1712:     TSSolve(ts, X);
1713:     TSGetSolveTime(ts, &ftime);
1714:     TSGetStepNumber(ts, &nsteps);
1715:     for (adaptIter = 0; ftime < finalTime; adaptIter++) {
1716:       PetscLogDouble bytes;

1718:       PetscMemoryGetCurrentUsage(&bytes);
1719:       PetscInfo(ts, "AMR time step loop %" PetscInt_FMT ": memory used %g\n", adaptIter, bytes);
1720:       PetscFVSetLimiter(fvm, noneLimiter);
1721:       adaptToleranceFVM(fvm, ts, X, refineTag, coarsenTag, user, &tsNew, &solNew);
1722:       PetscFVSetLimiter(fvm, limiter);
1723:       if (tsNew) {
1724:         PetscInfo(ts, "AMR used\n");
1725:         DMDestroy(&dm);
1726:         VecDestroy(&X);
1727:         TSDestroy(&ts);
1728:         ts = tsNew;
1729:         X  = solNew;
1730:         TSSetFromOptions(ts);
1731:         VecGetDM(X, &dm);
1732:         PetscObjectReference((PetscObject)dm);
1733:         DMConvert(dm, DMPLEX, &plex);
1734:         DMPlexGetGeometryFVM(dm, NULL, NULL, &minRadius);
1735:         DMDestroy(&plex);
1736:         MPI_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
1738:         dt = cfl * minRadius / mod->maxspeed;
1739:         TSSetStepNumber(ts, nsteps);
1740:         TSSetTime(ts, ftime);
1741:         TSSetTimeStep(ts, dt);
1742:       } else {
1743:         PetscInfo(ts, "AMR not used\n");
1744:       }
1745:       user->monitorStepOffset = nsteps;
1746:       TSSetMaxSteps(ts, nsteps + adaptInterval);
1747:       TSSolve(ts, X);
1748:       TSGetSolveTime(ts, &ftime);
1749:       TSGetStepNumber(ts, &nsteps);
1750:     }
1751:   }
1752:   TSGetConvergedReason(ts, &reason);
1753:   PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps);
1754:   TSDestroy(&ts);

1756:   VecTaggerDestroy(&refineTag);
1757:   VecTaggerDestroy(&coarsenTag);
1758:   PetscFunctionListDestroy(&PhysicsList);
1759:   PetscFunctionListDestroy(&PhysicsRiemannList_SW);
1760:   FunctionalLinkDestroy(&user->model->functionalRegistry);
1761:   PetscFree(user->model->functionalMonitored);
1762:   PetscFree(user->model->functionalCall);
1763:   PetscFree(user->model->physics->data);
1764:   PetscFree(user->model->physics);
1765:   PetscFree(user->model);
1766:   PetscFree(user);
1767:   VecDestroy(&X);
1768:   PetscLimiterDestroy(&limiter);
1769:   PetscLimiterDestroy(&noneLimiter);
1770:   PetscFVDestroy(&fvm);
1771:   DMDestroy(&dm);
1772:   PetscFinalize();
1773:   return 0;
1774: }

1776: /* Godunov fluxs */
1777: PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
1778: {
1779:   /* System generated locals */
1780:   PetscScalar ret_val;

1782:   if (PetscRealPart(*test) > 0.) goto L10;
1783:   ret_val = *b;
1784:   return ret_val;
1785: L10:
1786:   ret_val = *a;
1787:   return ret_val;
1788: } /* cvmgp_ */

1790: PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
1791: {
1792:   /* System generated locals */
1793:   PetscScalar ret_val;

1795:   if (PetscRealPart(*test) < 0.) goto L10;
1796:   ret_val = *b;
1797:   return ret_val;
1798: L10:
1799:   ret_val = *a;
1800:   return ret_val;
1801: } /* cvmgm_ */

1803: int riem1mdt(PetscScalar *gaml, PetscScalar *gamr, PetscScalar *rl, PetscScalar *pl, PetscScalar *uxl, PetscScalar *rr, PetscScalar *pr, PetscScalar *uxr, PetscScalar *rstarl, PetscScalar *rstarr, PetscScalar *pstar, PetscScalar *ustar)
1804: {
1805:   /* Initialized data */

1807:   static PetscScalar smallp = 1e-8;

1809:   /* System generated locals */
1810:   int         i__1;
1811:   PetscScalar d__1, d__2;

1813:   /* Local variables */
1814:   static int         i0;
1815:   static PetscScalar cl, cr, wl, zl, wr, zr, pst, durl, skpr1, skpr2;
1816:   static int         iwave;
1817:   static PetscScalar gascl4, gascr4, cstarl, dpstar, cstarr;
1818:   /* static PetscScalar csqrl, csqrr, gascl1, gascl2, gascl3, gascr1, gascr2, gascr3; */
1819:   static int         iterno;
1820:   static PetscScalar ustarl, ustarr, rarepr1, rarepr2;

1822:   /* gascl1 = *gaml - 1.; */
1823:   /* gascl2 = (*gaml + 1.) * .5; */
1824:   /* gascl3 = gascl2 / *gaml; */
1825:   gascl4 = 1. / (*gaml - 1.);

1827:   /* gascr1 = *gamr - 1.; */
1828:   /* gascr2 = (*gamr + 1.) * .5; */
1829:   /* gascr3 = gascr2 / *gamr; */
1830:   gascr4 = 1. / (*gamr - 1.);
1831:   iterno = 10;
1832:   /*        find pstar: */
1833:   cl = PetscSqrtScalar(*gaml * *pl / *rl);
1834:   cr = PetscSqrtScalar(*gamr * *pr / *rr);
1835:   wl = *rl * cl;
1836:   wr = *rr * cr;
1837:   /* csqrl = wl * wl; */
1838:   /* csqrr = wr * wr; */
1839:   *pstar  = (wl * *pr + wr * *pl) / (wl + wr);
1840:   *pstar  = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1841:   pst     = *pl / *pr;
1842:   skpr1   = cr * (pst - 1.) * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1843:   d__1    = (*gamr - 1.) / (*gamr * 2.);
1844:   rarepr2 = gascr4 * 2. * cr * (1. - PetscPowScalar(pst, d__1));
1845:   pst     = *pr / *pl;
1846:   skpr2   = cl * (pst - 1.) * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1847:   d__1    = (*gaml - 1.) / (*gaml * 2.);
1848:   rarepr1 = gascl4 * 2. * cl * (1. - PetscPowScalar(pst, d__1));
1849:   durl    = *uxr - *uxl;
1850:   if (PetscRealPart(*pr) < PetscRealPart(*pl)) {
1851:     if (PetscRealPart(durl) >= PetscRealPart(rarepr1)) {
1852:       iwave = 100;
1853:     } else if (PetscRealPart(durl) <= PetscRealPart(-skpr1)) {
1854:       iwave = 300;
1855:     } else {
1856:       iwave = 400;
1857:     }
1858:   } else {
1859:     if (PetscRealPart(durl) >= PetscRealPart(rarepr2)) {
1860:       iwave = 100;
1861:     } else if (PetscRealPart(durl) <= PetscRealPart(-skpr2)) {
1862:       iwave = 300;
1863:     } else {
1864:       iwave = 200;
1865:     }
1866:   }
1867:   if (iwave == 100) {
1868:     /*     1-wave: rarefaction wave, 3-wave: rarefaction wave */
1869:     /*     case (100) */
1870:     i__1 = iterno;
1871:     for (i0 = 1; i0 <= i__1; ++i0) {
1872:       d__1    = *pstar / *pl;
1873:       d__2    = 1. / *gaml;
1874:       *rstarl = *rl * PetscPowScalar(d__1, d__2);
1875:       cstarl  = PetscSqrtScalar(*gaml * *pstar / *rstarl);
1876:       ustarl  = *uxl - gascl4 * 2. * (cstarl - cl);
1877:       zl      = *rstarl * cstarl;
1878:       d__1    = *pstar / *pr;
1879:       d__2    = 1. / *gamr;
1880:       *rstarr = *rr * PetscPowScalar(d__1, d__2);
1881:       cstarr  = PetscSqrtScalar(*gamr * *pstar / *rstarr);
1882:       ustarr  = *uxr + gascr4 * 2. * (cstarr - cr);
1883:       zr      = *rstarr * cstarr;
1884:       dpstar  = zl * zr * (ustarr - ustarl) / (zl + zr);
1885:       *pstar -= dpstar;
1886:       *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1887:       if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
1888: #if 0
1889:         break;
1890: #endif
1891:       }
1892:     }
1893:     /*     1-wave: shock wave, 3-wave: rarefaction wave */
1894:   } else if (iwave == 200) {
1895:     /*     case (200) */
1896:     i__1 = iterno;
1897:     for (i0 = 1; i0 <= i__1; ++i0) {
1898:       pst     = *pstar / *pl;
1899:       ustarl  = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1900:       zl      = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
1901:       d__1    = *pstar / *pr;
1902:       d__2    = 1. / *gamr;
1903:       *rstarr = *rr * PetscPowScalar(d__1, d__2);
1904:       cstarr  = PetscSqrtScalar(*gamr * *pstar / *rstarr);
1905:       zr      = *rstarr * cstarr;
1906:       ustarr  = *uxr + gascr4 * 2. * (cstarr - cr);
1907:       dpstar  = zl * zr * (ustarr - ustarl) / (zl + zr);
1908:       *pstar -= dpstar;
1909:       *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1910:       if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
1911: #if 0
1912:         break;
1913: #endif
1914:       }
1915:     }
1916:     /*     1-wave: shock wave, 3-wave: shock */
1917:   } else if (iwave == 300) {
1918:     /*     case (300) */
1919:     i__1 = iterno;
1920:     for (i0 = 1; i0 <= i__1; ++i0) {
1921:       pst    = *pstar / *pl;
1922:       ustarl = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1923:       zl     = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
1924:       pst    = *pstar / *pr;
1925:       ustarr = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1926:       zr     = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
1927:       dpstar = zl * zr * (ustarr - ustarl) / (zl + zr);
1928:       *pstar -= dpstar;
1929:       *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1930:       if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
1931: #if 0
1932:         break;
1933: #endif
1934:       }
1935:     }
1936:     /*     1-wave: rarefaction wave, 3-wave: shock */
1937:   } else if (iwave == 400) {
1938:     /*     case (400) */
1939:     i__1 = iterno;
1940:     for (i0 = 1; i0 <= i__1; ++i0) {
1941:       d__1    = *pstar / *pl;
1942:       d__2    = 1. / *gaml;
1943:       *rstarl = *rl * PetscPowScalar(d__1, d__2);
1944:       cstarl  = PetscSqrtScalar(*gaml * *pstar / *rstarl);
1945:       ustarl  = *uxl - gascl4 * 2. * (cstarl - cl);
1946:       zl      = *rstarl * cstarl;
1947:       pst     = *pstar / *pr;
1948:       ustarr  = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1949:       zr      = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
1950:       dpstar  = zl * zr * (ustarr - ustarl) / (zl + zr);
1951:       *pstar -= dpstar;
1952:       *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1953:       if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
1954: #if 0
1955:               break;
1956: #endif
1957:       }
1958:     }
1959:   }

1961:   *ustar = (zl * ustarr + zr * ustarl) / (zl + zr);
1962:   if (PetscRealPart(*pstar) > PetscRealPart(*pl)) {
1963:     pst     = *pstar / *pl;
1964:     *rstarl = ((*gaml + 1.) * pst + *gaml - 1.) / ((*gaml - 1.) * pst + *gaml + 1.) * *rl;
1965:   }
1966:   if (PetscRealPart(*pstar) > PetscRealPart(*pr)) {
1967:     pst     = *pstar / *pr;
1968:     *rstarr = ((*gamr + 1.) * pst + *gamr - 1.) / ((*gamr - 1.) * pst + *gamr + 1.) * *rr;
1969:   }
1970:   return iwave;
1971: }

1973: PetscScalar sign(PetscScalar x)
1974: {
1975:   if (PetscRealPart(x) > 0) return 1.0;
1976:   if (PetscRealPart(x) < 0) return -1.0;
1977:   return 0.0;
1978: }
1979: /*        Riemann Solver */
1980: /* -------------------------------------------------------------------- */
1981: int riemannsolver(PetscScalar *xcen, PetscScalar *xp, PetscScalar *dtt, PetscScalar *rl, PetscScalar *uxl, PetscScalar *pl, PetscScalar *utl, PetscScalar *ubl, PetscScalar *gaml, PetscScalar *rho1l, PetscScalar *rr, PetscScalar *uxr, PetscScalar *pr, PetscScalar *utr, PetscScalar *ubr, PetscScalar *gamr, PetscScalar *rho1r, PetscScalar *rx, PetscScalar *uxm, PetscScalar *px, PetscScalar *utx, PetscScalar *ubx, PetscScalar *gam, PetscScalar *rho1)
1982: {
1983:   /* System generated locals */
1984:   PetscScalar d__1, d__2;

1986:   /* Local variables */
1987:   static PetscScalar s, c0, p0, r0, u0, w0, x0, x2, ri, cx, sgn0, wsp0, gasc1, gasc2, gasc3, gasc4;
1988:   static PetscScalar cstar, pstar, rstar, ustar, xstar, wspst, ushock, streng, rstarl, rstarr, rstars;
1989:   int                iwave;

1991:   if (*rl == *rr && *pr == *pl && *uxl == *uxr && *gaml == *gamr) {
1992:     *rx  = *rl;
1993:     *px  = *pl;
1994:     *uxm = *uxl;
1995:     *gam = *gaml;
1996:     x2   = *xcen + *uxm * *dtt;

1998:     if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
1999:       *utx  = *utr;
2000:       *ubx  = *ubr;
2001:       *rho1 = *rho1r;
2002:     } else {
2003:       *utx  = *utl;
2004:       *ubx  = *ubl;
2005:       *rho1 = *rho1l;
2006:     }
2007:     return 0;
2008:   }
2009:   iwave = riem1mdt(gaml, gamr, rl, pl, uxl, rr, pr, uxr, &rstarl, &rstarr, &pstar, &ustar);

2011:   x2   = *xcen + ustar * *dtt;
2012:   d__1 = *xp - x2;
2013:   sgn0 = sign(d__1);
2014:   /*            x is in 3-wave if sgn0 = 1 */
2015:   /*            x is in 1-wave if sgn0 = -1 */
2016:   r0     = cvmgm_(rl, rr, &sgn0);
2017:   p0     = cvmgm_(pl, pr, &sgn0);
2018:   u0     = cvmgm_(uxl, uxr, &sgn0);
2019:   *gam   = cvmgm_(gaml, gamr, &sgn0);
2020:   gasc1  = *gam - 1.;
2021:   gasc2  = (*gam + 1.) * .5;
2022:   gasc3  = gasc2 / *gam;
2023:   gasc4  = 1. / (*gam - 1.);
2024:   c0     = PetscSqrtScalar(*gam * p0 / r0);
2025:   streng = pstar - p0;
2026:   w0     = *gam * r0 * p0 * (gasc3 * streng / p0 + 1.);
2027:   rstars = r0 / (1. - r0 * streng / w0);
2028:   d__1   = p0 / pstar;
2029:   d__2   = -1. / *gam;
2030:   rstarr = r0 * PetscPowScalar(d__1, d__2);
2031:   rstar  = cvmgm_(&rstarr, &rstars, &streng);
2032:   w0     = PetscSqrtScalar(w0);
2033:   cstar  = PetscSqrtScalar(*gam * pstar / rstar);
2034:   wsp0   = u0 + sgn0 * c0;
2035:   wspst  = ustar + sgn0 * cstar;
2036:   ushock = ustar + sgn0 * w0 / rstar;
2037:   wspst  = cvmgp_(&ushock, &wspst, &streng);
2038:   wsp0   = cvmgp_(&ushock, &wsp0, &streng);
2039:   x0     = *xcen + wsp0 * *dtt;
2040:   xstar  = *xcen + wspst * *dtt;
2041:   /*           using gas formula to evaluate rarefaction wave */
2042:   /*            ri : reiman invariant */
2043:   ri   = u0 - sgn0 * 2. * gasc4 * c0;
2044:   cx   = sgn0 * .5 * gasc1 / gasc2 * ((*xp - *xcen) / *dtt - ri);
2045:   *uxm = ri + sgn0 * 2. * gasc4 * cx;
2046:   s    = p0 / PetscPowScalar(r0, *gam);
2047:   d__1 = cx * cx / (*gam * s);
2048:   *rx  = PetscPowScalar(d__1, gasc4);
2049:   *px  = cx * cx * *rx / *gam;
2050:   d__1 = sgn0 * (x0 - *xp);
2051:   *rx  = cvmgp_(rx, &r0, &d__1);
2052:   d__1 = sgn0 * (x0 - *xp);
2053:   *px  = cvmgp_(px, &p0, &d__1);
2054:   d__1 = sgn0 * (x0 - *xp);
2055:   *uxm = cvmgp_(uxm, &u0, &d__1);
2056:   d__1 = sgn0 * (xstar - *xp);
2057:   *rx  = cvmgm_(rx, &rstar, &d__1);
2058:   d__1 = sgn0 * (xstar - *xp);
2059:   *px  = cvmgm_(px, &pstar, &d__1);
2060:   d__1 = sgn0 * (xstar - *xp);
2061:   *uxm = cvmgm_(uxm, &ustar, &d__1);
2062:   if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
2063:     *utx  = *utr;
2064:     *ubx  = *ubr;
2065:     *rho1 = *rho1r;
2066:   } else {
2067:     *utx  = *utl;
2068:     *ubx  = *ubl;
2069:     *rho1 = *rho1l;
2070:   }
2071:   return iwave;
2072: }
2073: int godunovflux(const PetscScalar *ul, const PetscScalar *ur, PetscScalar *flux, const PetscReal *nn, const int *ndim, const PetscReal *gamma)
2074: {
2075:   /* System generated locals */
2076:   int         i__1, iwave;
2077:   PetscScalar d__1, d__2, d__3;

2079:   /* Local variables */
2080:   static int         k;
2081:   static PetscScalar bn[3], fn, ft, tg[3], pl, rl, pm, pr, rr, xp, ubl, ubm, ubr, dtt, unm, tmp, utl, utm, uxl, utr, uxr, gaml, gamm, gamr, xcen, rhom, rho1l, rho1m, rho1r;

2083:   /* Function Body */
2084:   xcen = 0.;
2085:   xp   = 0.;
2086:   i__1 = *ndim;
2087:   for (k = 1; k <= i__1; ++k) {
2088:     tg[k - 1] = 0.;
2089:     bn[k - 1] = 0.;
2090:   }
2091:   dtt = 1.;
2092:   if (*ndim == 3) {
2093:     if (nn[0] == 0. && nn[1] == 0.) {
2094:       tg[0] = 1.;
2095:     } else {
2096:       tg[0] = -nn[1];
2097:       tg[1] = nn[0];
2098:     }
2099:     /*           tmp=dsqrt(tg(1)**2+tg(2)**2) */
2100:     /*           tg=tg/tmp */
2101:     bn[0] = -nn[2] * tg[1];
2102:     bn[1] = nn[2] * tg[0];
2103:     bn[2] = nn[0] * tg[1] - nn[1] * tg[0];
2104:     /* Computing 2nd power */
2105:     d__1 = bn[0];
2106:     /* Computing 2nd power */
2107:     d__2 = bn[1];
2108:     /* Computing 2nd power */
2109:     d__3 = bn[2];
2110:     tmp  = PetscSqrtScalar(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
2111:     i__1 = *ndim;
2112:     for (k = 1; k <= i__1; ++k) bn[k - 1] /= tmp;
2113:   } else if (*ndim == 2) {
2114:     tg[0] = -nn[1];
2115:     tg[1] = nn[0];
2116:     /*           tmp=dsqrt(tg(1)**2+tg(2)**2) */
2117:     /*           tg=tg/tmp */
2118:     bn[0] = 0.;
2119:     bn[1] = 0.;
2120:     bn[2] = 1.;
2121:   }
2122:   rl   = ul[0];
2123:   rr   = ur[0];
2124:   uxl  = 0.;
2125:   uxr  = 0.;
2126:   utl  = 0.;
2127:   utr  = 0.;
2128:   ubl  = 0.;
2129:   ubr  = 0.;
2130:   i__1 = *ndim;
2131:   for (k = 1; k <= i__1; ++k) {
2132:     uxl += ul[k] * nn[k - 1];
2133:     uxr += ur[k] * nn[k - 1];
2134:     utl += ul[k] * tg[k - 1];
2135:     utr += ur[k] * tg[k - 1];
2136:     ubl += ul[k] * bn[k - 1];
2137:     ubr += ur[k] * bn[k - 1];
2138:   }
2139:   uxl /= rl;
2140:   uxr /= rr;
2141:   utl /= rl;
2142:   utr /= rr;
2143:   ubl /= rl;
2144:   ubr /= rr;

2146:   gaml = *gamma;
2147:   gamr = *gamma;
2148:   /* Computing 2nd power */
2149:   d__1 = uxl;
2150:   /* Computing 2nd power */
2151:   d__2 = utl;
2152:   /* Computing 2nd power */
2153:   d__3 = ubl;
2154:   pl   = (*gamma - 1.) * (ul[*ndim + 1] - rl * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
2155:   /* Computing 2nd power */
2156:   d__1 = uxr;
2157:   /* Computing 2nd power */
2158:   d__2 = utr;
2159:   /* Computing 2nd power */
2160:   d__3  = ubr;
2161:   pr    = (*gamma - 1.) * (ur[*ndim + 1] - rr * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
2162:   rho1l = rl;
2163:   rho1r = rr;

2165:   iwave = riemannsolver(&xcen, &xp, &dtt, &rl, &uxl, &pl, &utl, &ubl, &gaml, &rho1l, &rr, &uxr, &pr, &utr, &ubr, &gamr, &rho1r, &rhom, &unm, &pm, &utm, &ubm, &gamm, &rho1m);

2167:   flux[0] = rhom * unm;
2168:   fn      = rhom * unm * unm + pm;
2169:   ft      = rhom * unm * utm;
2170:   /*           flux(2)=fn*nn(1)+ft*nn(2) */
2171:   /*           flux(3)=fn*tg(1)+ft*tg(2) */
2172:   flux[1] = fn * nn[0] + ft * tg[0];
2173:   flux[2] = fn * nn[1] + ft * tg[1];
2174:   /*           flux(2)=rhom*unm*(unm)+pm */
2175:   /*           flux(3)=rhom*(unm)*utm */
2176:   if (*ndim == 3) flux[3] = rhom * unm * ubm;
2177:   flux[*ndim + 1] = (rhom * .5 * (unm * unm + utm * utm + ubm * ubm) + gamm / (gamm - 1.) * pm) * unm;
2178:   return iwave;
2179: } /* godunovflux_ */

2181: /* Subroutine to set up the initial conditions for the */
2182: /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
2183: /* ----------------------------------------------------------------------- */
2184: int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
2185: {
2186:   int j, k;
2187:   /*      Wc=matmul(lv,Ueq) 3 vars */
2188:   for (k = 0; k < 3; ++k) {
2189:     wc[k] = 0.;
2190:     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
2191:   }
2192:   return 0;
2193: }
2194: /* ----------------------------------------------------------------------- */
2195: int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
2196: {
2197:   int k, j;
2198:   /*      V=matmul(rv,WC) 3 vars */
2199:   for (k = 0; k < 3; ++k) {
2200:     v[k] = 0.;
2201:     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
2202:   }
2203:   return 0;
2204: }
2205: /* ---------------------------------------------------------------------- */
2206: int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
2207: {
2208:   int       j, k;
2209:   PetscReal rho, csnd, p0;
2210:   /* PetscScalar u; */

2212:   for (k = 0; k < 3; ++k)
2213:     for (j = 0; j < 3; ++j) {
2214:       lv[k][j] = 0.;
2215:       rv[k][j] = 0.;
2216:     }
2217:   rho = ueq[0];
2218:   /* u = ueq[1]; */
2219:   p0       = ueq[2];
2220:   csnd     = PetscSqrtReal(gamma * p0 / rho);
2221:   lv[0][1] = rho * .5;
2222:   lv[0][2] = -.5 / csnd;
2223:   lv[1][0] = csnd;
2224:   lv[1][2] = -1. / csnd;
2225:   lv[2][1] = rho * .5;
2226:   lv[2][2] = .5 / csnd;
2227:   rv[0][0] = -1. / csnd;
2228:   rv[1][0] = 1. / rho;
2229:   rv[2][0] = -csnd;
2230:   rv[0][1] = 1. / csnd;
2231:   rv[0][2] = 1. / csnd;
2232:   rv[1][2] = 1. / rho;
2233:   rv[2][2] = csnd;
2234:   return 0;
2235: }

2237: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
2238: {
2239:   PetscReal p0, u0, wcp[3], wc[3];
2240:   PetscReal lv[3][3];
2241:   PetscReal vp[3];
2242:   PetscReal rv[3][3];
2243:   PetscReal eps, ueq[3], rho0, twopi;

2245:   /* Function Body */
2246:   twopi  = 2. * PETSC_PI;
2247:   eps    = 1e-4;    /* perturbation */
2248:   rho0   = 1e3;     /* density of water */
2249:   p0     = 101325.; /* init pressure of 1 atm (?) */
2250:   u0     = 0.;
2251:   ueq[0] = rho0;
2252:   ueq[1] = u0;
2253:   ueq[2] = p0;
2254:   /* Project initial state to characteristic variables */
2255:   eigenvectors(rv, lv, ueq, gamma);
2256:   projecteqstate(wc, ueq, lv);
2257:   wcp[0] = wc[0];
2258:   wcp[1] = wc[1];
2259:   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
2260:   projecttoprim(vp, wcp, rv);
2261:   ux->r     = vp[0];         /* density */
2262:   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
2263:   ux->ru[1] = 0.;
2264: #if defined DIM > 2
2265:   if (dim > 2) ux->ru[2] = 0.;
2266: #endif
2267:   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
2268:   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
2269:   return 0;
2270: }

2272: /*TEST

2274:   testset:
2275:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0

2277:     test:
2278:       suffix: adv_2d_tri_0
2279:       requires: triangle
2280:       TODO: how did this ever get in main when there is no support for this
2281:       args: -ufv_vtk_interval 0 -simplex -dm_refine 3 -dm_plex_faces 1,1 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2283:     test:
2284:       suffix: adv_2d_tri_1
2285:       requires: triangle
2286:       TODO: how did this ever get in main when there is no support for this
2287:       args: -ufv_vtk_interval 0 -simplex -dm_refine 5 -dm_plex_faces 1,1 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1

2289:     test:
2290:       suffix: tut_1
2291:       requires: exodusii
2292:       nsize: 1
2293:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

2295:     test:
2296:       suffix: tut_2
2297:       requires: exodusii
2298:       nsize: 1
2299:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw

2301:     test:
2302:       suffix: tut_3
2303:       requires: exodusii
2304:       nsize: 4
2305:       args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -monitor Error -advect_sol_type bump -petscfv_type leastsquares -petsclimiter_type sin

2307:     test:
2308:       suffix: tut_4
2309:       requires: exodusii
2310:       nsize: 4
2311:       args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -physics sw -monitor Height,Energy -petscfv_type leastsquares -petsclimiter_type minmod

2313:   testset:
2314:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

2316:     # 2D Advection 0-10
2317:     test:
2318:       suffix: 0
2319:       requires: exodusii
2320:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

2322:     test:
2323:       suffix: 1
2324:       requires: exodusii
2325:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo

2327:     test:
2328:       suffix: 2
2329:       requires: exodusii
2330:       nsize: 2
2331:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

2333:     test:
2334:       suffix: 3
2335:       requires: exodusii
2336:       nsize: 2
2337:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo

2339:     test:
2340:       suffix: 4
2341:       requires: exodusii
2342:       nsize: 8
2343:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo

2345:     test:
2346:       suffix: 5
2347:       requires: exodusii
2348:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1

2350:     test:
2351:       suffix: 7
2352:       requires: exodusii
2353:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1

2355:     test:
2356:       suffix: 8
2357:       requires: exodusii
2358:       nsize: 2
2359:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1

2361:     test:
2362:       suffix: 9
2363:       requires: exodusii
2364:       nsize: 8
2365:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1

2367:     test:
2368:       suffix: 10
2369:       requires: exodusii
2370:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo

2372:   # 2D Shallow water
2373:   testset:
2374:     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0

2376:     test:
2377:       suffix: sw_0
2378:       requires: exodusii
2379:       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
2380:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
2381:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2382:             -monitor height,energy

2384:     test:
2385:       suffix: sw_1
2386:       nsize: 2
2387:       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
2388:             -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_faces 24,12 -dm_plex_box_lower 0,1 -dm_plex_box_upper 1,3 -dm_distribute_overlap 1 \
2389:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2390:             -monitor height,energy

2392:     test:
2393:       suffix: sw_hll
2394:       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
2395:             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
2396:             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2397:             -monitor height,energy

2399:   testset:
2400:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

2402:     # 2D Advection: p4est
2403:     test:
2404:       suffix: p4est_advec_2d
2405:       requires: p4est
2406:       args: -ufv_vtk_interval 0 -dm_type p4est -dm_forest_minimum_refinement 1 -dm_forest_initial_refinement 2 -dm_p4est_refine_pattern hash   -dm_forest_maximum_refinement 5

2408:     # Advection in a box
2409:     test:
2410:       suffix: adv_2d_quad_0
2411:       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2413:     test:
2414:       suffix: adv_2d_quad_1
2415:       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
2416:       timeoutfactor: 3

2418:     test:
2419:       suffix: adv_2d_quad_p4est_0
2420:       requires: p4est
2421:       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2423:     test:
2424:       suffix: adv_2d_quad_p4est_1
2425:       requires: p4est
2426:       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow   3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
2427:       timeoutfactor: 3

2429:     test:
2430:       suffix: adv_2d_quad_p4est_adapt_0
2431:       requires: p4est !__float128 #broken for quad precision
2432:       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow   3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 -ufv_use_amr -refine_vec_tagger_box 0.005,inf -coarsen_vec_tagger_box   0,1.e-5 -petscfv_type leastsquares -ts_max_time 0.01
2433:       timeoutfactor: 3

2435:     test:
2436:       suffix: adv_0
2437:       requires: exodusii
2438:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -bc_inflow 100,101,200 -bc_outflow 201

2440:     test:
2441:       suffix: shock_0
2442:       requires: p4est !single !complex
2443:       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
2444:       -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 6 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
2445:       -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
2446:       -bc_wall 1,2,3,4 -physics euler -eu_type iv_shock -ufv_cfl 10 -eu_alpha 60. -eu_gamma 1.4 -eu_amach 2.02 -eu_rho2 3. \
2447:       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
2448:       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2449:       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
2450:       timeoutfactor: 3

2452:     # Test GLVis visualization of PetscFV fields
2453:     test:
2454:       suffix: glvis_adv_2d_tet
2455:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
2456:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
2457:             -ts_monitor_solution glvis: -ts_max_steps 0

2459:     test:
2460:       suffix: glvis_adv_2d_quad
2461:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
2462:             -dm_refine 5 -dm_plex_separate_marker \
2463:             -ts_monitor_solution glvis: -ts_max_steps 0

2465: TEST*/