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.

 37: The example also shows how to handle AMR in a time-dependent TS solver.
 38: F*/
 39: #include <petscdmplex.h>
 40: #include <petscdmforest.h>
 41: #include <petscds.h>
 42: #include <petscts.h>

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

 46: static PetscFunctionList PhysicsList, PhysicsRiemannList_SW;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

198:   PetscFunctionBeginUser;
199:   xG[0] = advect->inflowState;
200:   PetscFunctionReturn(PETSC_SUCCESS);
201: }

203: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
204: {
205:   PetscFunctionBeginUser;
206:   xG[0] = xI[0];
207:   PetscFunctionReturn(PETSC_SUCCESS);
208: }

210: 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)
211: {
212:   Physics_Advect *advect = (Physics_Advect *)phys->data;
213:   PetscReal       wind[DIM], wn;

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

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

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

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

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

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

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

305:   PetscFunctionBeginUser;
306:   PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
307:   f[advect->functional.Solution] = PetscRealPart(y[0]);
308:   f[advect->functional.Error]    = PetscAbsScalar(y[0] - yexact[0]);
309:   PetscFunctionReturn(PETSC_SUCCESS);
310: }

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

317:   PetscFunctionBeginUser;
318:   /* Register "canned" boundary conditions and defaults for where to apply. */
319:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
320:   PetscCall(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));
321:   PetscCall(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));
322:   PetscFunctionReturn(PETSC_SUCCESS);
323: }

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

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

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

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

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

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

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

421: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
422: {
423:   PetscFunctionBeginUser;
424:   xG[0] = xI[0];
425:   xG[1] = -xI[1];
426:   xG[2] = -xI[2];
427:   PetscFunctionReturn(PETSC_SUCCESS);
428: }

430: 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)
431: {
432:   Physics_SW *sw = (Physics_SW *)phys->data;
433:   PetscReal   aL, aR;
434:   PetscReal   nn[DIM];
435: #if !defined(PETSC_USE_COMPLEX)
436:   const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
437: #else
438:   SWNodeUnion   uLreal, uRreal;
439:   const SWNode *uL = &uLreal.swnode;
440:   const SWNode *uR = &uRreal.swnode;
441: #endif
442:   SWNodeUnion fL, fR;
443:   PetscInt    i;
444:   PetscReal   zero = 0.;

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

480: 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)
481: {
482:   Physics_SW *sw = (Physics_SW *)phys->data;
483:   PetscReal   cL, cR, speed;
484:   PetscReal   nn[DIM];
485: #if !defined(PETSC_USE_COMPLEX)
486:   const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
487: #else
488:   SWNodeUnion   uLreal, uRreal;
489:   const SWNode *uL = &uLreal.swnode;
490:   const SWNode *uR = &uRreal.swnode;
491: #endif
492:   SWNodeUnion fL, fR;
493:   PetscInt    i;
494:   PetscReal   zero = 0.;

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

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

521:   PetscFunctionBeginUser;
522:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
523:   dx[0] = x[0] - 1.5;
524:   dx[1] = x[1] - 1.0;
525:   r     = Norm2Real(dx);
526:   sigma = 0.5;
527:   u[0]  = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma)));
528:   u[1]  = 0.0;
529:   u[2]  = 0.0;
530:   PetscFunctionReturn(PETSC_SUCCESS);
531: }

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

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

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

555:   PetscFunctionBeginUser;
556:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
557:   PetscCall(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));
558:   PetscFunctionReturn(PETSC_SUCCESS);
559: }

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

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

572:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov));
573:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL));

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

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

592:   PetscFunctionReturn(PETSC_SUCCESS);
593: }

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

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

643: /* initial condition */
644: int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
645: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
646: {
647:   PetscInt       i;
648:   Physics        phys = (Physics)ctx;
649:   Physics_Euler *eu   = (Physics_Euler *)phys->data;
650:   EulerNode     *uu   = (EulerNode *)u;
651:   PetscReal      p0, gamma, c;
652:   PetscFunctionBeginUser;
653:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);

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

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

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

703:   PetscFunctionReturn(PETSC_SUCCESS);
704: }

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

710:   PetscFunctionBeginUser;
711:   ru2  = DotDIMReal(x->ru, x->ru);
712:   (*p) = (x->E - 0.5 * ru2 / x->r) * (gamma - 1.0); /* (E - rho V^2/2)(gamma-1) = e rho (gamma-1) */
713:   PetscFunctionReturn(PETSC_SUCCESS);
714: }

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

720:   PetscFunctionBeginUser;
721:   PetscCall(Pressure_PG(*gamma, x, &p));
722:   PetscCheck(p >= 0., PETSC_COMM_WORLD, PETSC_ERR_SUP, "negative pressure time %g -- NEED TO FIX!!!!!!", (double)p);
723:   /* pars[EULER_PAR_GAMMA] = heat capacity ratio */
724:   (*c) = PetscSqrtReal(*gamma * p / x->r);
725:   PetscFunctionReturn(PETSC_SUCCESS);
726: }

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

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

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

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

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

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

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

835:   PetscFunctionBeginUser;
836:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
837:   if (eu->type == EULER_LINEAR_WAVE) {
838:     const PetscInt wallids[] = {100, 101};
839:     PetscCall(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));
840:   } else {
841:     const PetscInt wallids[] = {100, 101, 200, 300};
842:     PetscCall(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));
843:   }
844:   PetscFunctionReturn(PETSC_SUCCESS);
845: }

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

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

909:   PetscFunctionReturn(PETSC_SUCCESS);
910: }

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

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

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

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

957:     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
958:     p[0] = rank;
959:   }
960:   PetscCall(VecRestoreArray(*partition, &part));
961:   PetscFunctionReturn(PETSC_SUCCESS);
962: }

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

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

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

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

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

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

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

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

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

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

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

1108: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
1109: {
1110:   FunctionalLink l, next;

1112:   PetscFunctionBeginUser;
1113:   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
1114:   l     = *link;
1115:   *link = NULL;
1116:   for (; l; l = next) {
1117:     next = l->next;
1118:     PetscCall(PetscFree(l->name));
1119:     PetscCall(PetscFree(l));
1120:   }
1121:   PetscFunctionReturn(PETSC_SUCCESS);
1122: }

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

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

1140:   PetscFunctionBeginUser;
1141:   func[0] = SolutionFunctional;
1142:   ctx[0]  = (void *)mod;
1143:   PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X));
1144:   PetscFunctionReturn(PETSC_SUCCESS);
1145: }

1147: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1148: {
1149:   PetscFunctionBeginUser;
1150:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1151:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
1152:   PetscCall(PetscViewerFileSetName(*viewer, filename));
1153:   PetscFunctionReturn(PETSC_SUCCESS);
1154: }

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

1164:   PetscFunctionBeginUser;
1165:   PetscCall(PetscObjectSetName((PetscObject)X, "u"));
1166:   PetscCall(VecGetDM(X, &dm));
1167:   PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));

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

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

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

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

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

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

1287: typedef struct {
1288:   PetscFV      fvm;
1289:   VecTagger    refineTag;
1290:   VecTagger    coarsenTag;
1291:   DM           adaptedDM;
1292:   User         user;
1293:   PetscReal    cfl;
1294:   PetscLimiter limiter;
1295:   PetscLimiter noneLimiter;
1296: } TransferCtx;

1298: static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
1299: {
1300:   TransferCtx       *tctx       = (TransferCtx *)ctx;
1301:   PetscFV            fvm        = tctx->fvm;
1302:   VecTagger          refineTag  = tctx->refineTag;
1303:   VecTagger          coarsenTag = tctx->coarsenTag;
1304:   User               user       = tctx->user;
1305:   DM                 dm, gradDM, plex, cellDM, adaptedDM = NULL;
1306:   Vec                cellGeom, faceGeom;
1307:   PetscBool          isForest, computeGradient;
1308:   Vec                grad, locGrad, locX, errVec;
1309:   PetscInt           cStart, cEnd, c, dim, nRefine, nCoarsen;
1310:   PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1311:   PetscScalar       *errArray;
1312:   const PetscScalar *pointVals;
1313:   const PetscScalar *pointGrads;
1314:   const PetscScalar *pointGeom;
1315:   DMLabel            adaptLabel = NULL;
1316:   IS                 refineIS, coarsenIS;

1318:   PetscFunctionBeginUser;
1319:   *resize = PETSC_FALSE;
1320:   PetscCall(VecGetDM(sol, &dm));
1321:   PetscCall(DMGetDimension(dm, &dim));
1322:   PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
1323:   PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
1324:   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
1325:   PetscCall(DMIsForest(dm, &isForest));
1326:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1327:   PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
1328:   PetscCall(DMCreateLocalVector(plex, &locX));
1329:   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
1330:   PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
1331:   PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
1332:   PetscCall(DMCreateGlobalVector(gradDM, &grad));
1333:   PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
1334:   PetscCall(DMCreateLocalVector(gradDM, &locGrad));
1335:   PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
1336:   PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
1337:   PetscCall(VecDestroy(&grad));
1338:   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
1339:   PetscCall(VecGetArrayRead(locGrad, &pointGrads));
1340:   PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
1341:   PetscCall(VecGetArrayRead(locX, &pointVals));
1342:   PetscCall(VecGetDM(cellGeom, &cellDM));
1343:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
1344:   PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
1345:   PetscCall(VecSetUp(errVec));
1346:   PetscCall(VecGetArray(errVec, &errArray));
1347:   for (c = cStart; c < cEnd; c++) {
1348:     PetscReal        errInd = 0.;
1349:     PetscScalar     *pointGrad;
1350:     PetscScalar     *pointVal;
1351:     PetscFVCellGeom *cg;

1353:     PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
1354:     PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
1355:     PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));

1357:     PetscCall((user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1358:     errArray[c - cStart] = errInd;
1359:     minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
1360:     minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
1361:   }
1362:   PetscCall(VecRestoreArray(errVec, &errArray));
1363:   PetscCall(VecRestoreArrayRead(locX, &pointVals));
1364:   PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
1365:   PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
1366:   PetscCall(VecDestroy(&locGrad));
1367:   PetscCall(VecDestroy(&locX));
1368:   PetscCall(DMDestroy(&plex));

1370:   PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
1371:   PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
1372:   PetscCall(ISGetSize(refineIS, &nRefine));
1373:   PetscCall(ISGetSize(coarsenIS, &nCoarsen));
1374:   if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
1375:   if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
1376:   PetscCall(ISDestroy(&coarsenIS));
1377:   PetscCall(ISDestroy(&refineIS));
1378:   PetscCall(VecDestroy(&errVec));

1380:   PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
1381:   PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1382:   minMaxInd[1] = -minMaxInd[1];
1383:   PetscCall(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1384:   PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1385:   if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
1386:     PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1387:   }
1388:   PetscCall(DMLabelDestroy(&adaptLabel));
1389:   if (adaptedDM) {
1390:     tctx->adaptedDM = adaptedDM;
1391:     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nstep, nRefine, nCoarsen));
1392:     *resize = PETSC_TRUE;
1393:   } else {
1394:     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1395:   }
1396:   PetscFunctionReturn(PETSC_SUCCESS);
1397: }

1399: static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
1400: {
1401:   TransferCtx *tctx = (TransferCtx *)ctx;
1402:   DM           dm;
1403:   PetscReal    time;

1405:   PetscFunctionBeginUser;
1406:   PetscCall(TSGetDM(ts, &dm));
1407:   PetscCall(TSGetTime(ts, &time));
1408:   PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
1409:   for (PetscInt i = 0; i < nv; i++) {
1410:     const char *name;

1412:     PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
1413:     PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
1414:     PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
1415:     PetscCall(PetscObjectSetName((PetscObject)vecsout[i], name));
1416:   }
1417:   PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */

1419:   Model     mod  = tctx->user->model;
1420:   Physics   phys = mod->physics;
1421:   PetscReal minRadius;

1423:   PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius));
1424:   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1425:   PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed");

1427:   PetscReal dt = tctx->cfl * minRadius / mod->maxspeed;
1428:   PetscCall(TSSetTimeStep(ts, dt));

1430:   PetscCall(TSSetDM(ts, tctx->adaptedDM));
1431:   PetscCall(DMDestroy(&tctx->adaptedDM));

1433:   PetscFunctionReturn(PETSC_SUCCESS);
1434: }

1436: int main(int argc, char **argv)
1437: {
1438:   MPI_Comm          comm;
1439:   PetscDS           prob;
1440:   PetscFV           fvm;
1441:   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1442:   User              user;
1443:   Model             mod;
1444:   Physics           phys;
1445:   DM                dm, plex;
1446:   PetscReal         ftime, cfl, dt, minRadius;
1447:   PetscInt          dim, nsteps;
1448:   TS                ts;
1449:   TSConvergedReason reason;
1450:   Vec               X;
1451:   PetscViewer       viewer;
1452:   PetscBool         vtkCellGeom, useAMR;
1453:   PetscInt          adaptInterval;
1454:   char              physname[256] = "advect";
1455:   VecTagger         refineTag = NULL, coarsenTag = NULL;
1456:   TransferCtx       tctx;

1458:   PetscFunctionBeginUser;
1459:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1460:   comm = PETSC_COMM_WORLD;

1462:   PetscCall(PetscNew(&user));
1463:   PetscCall(PetscNew(&user->model));
1464:   PetscCall(PetscNew(&user->model->physics));
1465:   mod           = user->model;
1466:   phys          = mod->physics;
1467:   mod->comm     = comm;
1468:   useAMR        = PETSC_FALSE;
1469:   adaptInterval = 1;

1471:   /* Register physical models to be available on the command line */
1472:   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
1473:   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
1474:   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));

1476:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1477:   {
1478:     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1479:     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1480:     user->vtkInterval = 1;
1481:     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1482:     user->vtkmon = PETSC_TRUE;
1483:     PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1484:     vtkCellGeom = PETSC_FALSE;
1485:     PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
1486:     PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
1487:     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
1488:     PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
1489:     PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1490:   }
1491:   PetscOptionsEnd();

1493:   if (useAMR) {
1494:     VecTaggerBox refineBox, coarsenBox;

1496:     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1497:     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;

1499:     PetscCall(VecTaggerCreate(comm, &refineTag));
1500:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
1501:     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
1502:     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
1503:     PetscCall(VecTaggerSetFromOptions(refineTag));
1504:     PetscCall(VecTaggerSetUp(refineTag));
1505:     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));

1507:     PetscCall(VecTaggerCreate(comm, &coarsenTag));
1508:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
1509:     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
1510:     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
1511:     PetscCall(VecTaggerSetFromOptions(coarsenTag));
1512:     PetscCall(VecTaggerSetUp(coarsenTag));
1513:     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1514:   }

1516:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1517:   {
1518:     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
1519:     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
1520:     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
1521:     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
1522:     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1523:     /* Count number of fields and dofs */
1524:     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1525:     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1526:     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1527:   }
1528:   PetscOptionsEnd();

1530:   /* Create mesh */
1531:   {
1532:     PetscInt i;

1534:     PetscCall(DMCreate(comm, &dm));
1535:     PetscCall(DMSetType(dm, DMPLEX));
1536:     PetscCall(DMSetFromOptions(dm));
1537:     for (i = 0; i < DIM; i++) {
1538:       mod->bounds[2 * i]     = 0.;
1539:       mod->bounds[2 * i + 1] = 1.;
1540:     };
1541:     dim = DIM;
1542:     { /* a null name means just do a hex box */
1543:       PetscInt  cells[3] = {1, 1, 1}, n = 3;
1544:       PetscBool flg2, skew              = PETSC_FALSE;
1545:       PetscInt  nret2 = 2 * DIM;
1546:       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
1547:       PetscCall(PetscOptionsRealArray("-grid_bounds", "bounds of the mesh in each direction (i.e., x_min,x_max,y_min,y_max", "", mod->bounds, &nret2, &flg2));
1548:       PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
1549:       PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1550:       PetscOptionsEnd();
1551:       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1552:       if (flg2) {
1553:         PetscInt     dimEmbed, i;
1554:         PetscInt     nCoords;
1555:         PetscScalar *coords;
1556:         Vec          coordinates;

1558:         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1559:         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1560:         PetscCall(VecGetLocalSize(coordinates, &nCoords));
1561:         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
1562:         PetscCall(VecGetArray(coordinates, &coords));
1563:         for (i = 0; i < nCoords; i += dimEmbed) {
1564:           PetscInt j;

1566:           PetscScalar *coord = &coords[i];
1567:           for (j = 0; j < dimEmbed; j++) {
1568:             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1569:             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1570:               if (cells[0] == 2 && i == 8) {
1571:                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
1572:               } else if (cells[0] == 3) {
1573:                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1574:                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1575:                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1576:               }
1577:             }
1578:           }
1579:         }
1580:         PetscCall(VecRestoreArray(coordinates, &coords));
1581:         PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1582:       }
1583:     }
1584:   }
1585:   PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
1586:   PetscCall(DMGetDimension(dm, &dim));

1588:   /* set up BCs, functions, tags */
1589:   PetscCall(DMCreateLabel(dm, "Face Sets"));
1590:   mod->errorIndicator = ErrorIndicator_Simple;

1592:   {
1593:     DM gdm;

1595:     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
1596:     PetscCall(DMDestroy(&dm));
1597:     dm = gdm;
1598:     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1599:   }

1601:   PetscCall(PetscFVCreate(comm, &fvm));
1602:   PetscCall(PetscFVSetFromOptions(fvm));
1603:   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
1604:   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
1605:   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1606:   {
1607:     PetscInt f, dof;
1608:     for (f = 0, dof = 0; f < phys->nfields; f++) {
1609:       PetscInt newDof = phys->field_desc[f].dof;

1611:       if (newDof == 1) {
1612:         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
1613:       } else {
1614:         PetscInt j;

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

1619:           PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
1620:           PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1621:         }
1622:       }
1623:       dof += newDof;
1624:     }
1625:   }
1626:   /* FV is now structured with one field having all physics as components */
1627:   PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
1628:   PetscCall(DMCreateDS(dm));
1629:   PetscCall(DMGetDS(dm, &prob));
1630:   PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
1631:   PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
1632:   PetscCall((*mod->setupbc)(dm, prob, phys));
1633:   PetscCall(PetscDSSetFromOptions(prob));
1634:   {
1635:     char      convType[256];
1636:     PetscBool flg;

1638:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
1639:     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1640:     PetscOptionsEnd();
1641:     if (flg) {
1642:       DM dmConv;

1644:       PetscCall(DMConvert(dm, convType, &dmConv));
1645:       if (dmConv) {
1646:         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
1647:         PetscCall(DMDestroy(&dm));
1648:         dm = dmConv;
1649:         PetscCall(DMSetFromOptions(dm));
1650:       }
1651:     }
1652:   }

1654:   PetscCall(initializeTS(dm, user, &ts));

1656:   PetscCall(DMCreateGlobalVector(dm, &X));
1657:   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1658:   PetscCall(SetInitialCondition(dm, X, user));
1659:   if (useAMR) {
1660:     PetscInt adaptIter;

1662:     /* use no limiting when reconstructing gradients for adaptivity */
1663:     PetscCall(PetscFVGetLimiter(fvm, &limiter));
1664:     PetscCall(PetscObjectReference((PetscObject)limiter));
1665:     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
1666:     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));

1668:     /* Refinement context */
1669:     tctx.fvm         = fvm;
1670:     tctx.refineTag   = refineTag;
1671:     tctx.coarsenTag  = coarsenTag;
1672:     tctx.adaptedDM   = NULL;
1673:     tctx.user        = user;
1674:     tctx.noneLimiter = noneLimiter;
1675:     tctx.limiter     = limiter;
1676:     tctx.cfl         = cfl;

1678:     /* Do some initial refinement steps */
1679:     for (adaptIter = 0;; ++adaptIter) {
1680:       PetscLogDouble bytes;
1681:       PetscBool      resize;

1683:       PetscCall(PetscMemoryGetCurrentUsage(&bytes));
1684:       PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
1685:       PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
1686:       PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));

1688:       PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
1689:       if (!resize) break;
1690:       PetscCall(DMDestroy(&dm));
1691:       PetscCall(VecDestroy(&X));
1692:       dm             = tctx.adaptedDM;
1693:       tctx.adaptedDM = NULL;
1694:       PetscCall(TSSetDM(ts, dm));
1695:       PetscCall(DMCreateGlobalVector(dm, &X));
1696:       PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1697:       PetscCall(SetInitialCondition(dm, X, user));
1698:     }
1699:   }

1701:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1702:   if (vtkCellGeom) {
1703:     DM  dmCell;
1704:     Vec cellgeom, partition;

1706:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1707:     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
1708:     PetscCall(VecView(cellgeom, viewer));
1709:     PetscCall(PetscViewerDestroy(&viewer));
1710:     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
1711:     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
1712:     PetscCall(VecView(partition, viewer));
1713:     PetscCall(PetscViewerDestroy(&viewer));
1714:     PetscCall(VecDestroy(&partition));
1715:     PetscCall(DMDestroy(&dmCell));
1716:   }
1717:   /* collect max maxspeed from all processes -- todo */
1718:   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
1719:   PetscCall(DMDestroy(&plex));
1720:   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1721:   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1722:   dt = cfl * minRadius / mod->maxspeed;
1723:   PetscCall(TSSetTimeStep(ts, dt));
1724:   PetscCall(TSSetFromOptions(ts));

1726:   /* When using adaptive mesh refinement
1727:      specify callbacks to refine the solution
1728:      and interpolate data from old to new mesh */
1729:   if (useAMR) { PetscCall(TSSetResize(ts, adaptToleranceFVMSetUp, Transfer, &tctx)); }
1730:   PetscCall(TSSetSolution(ts, X));
1731:   PetscCall(VecDestroy(&X));
1732:   PetscCall(TSSolve(ts, NULL));
1733:   PetscCall(TSGetSolveTime(ts, &ftime));
1734:   PetscCall(TSGetStepNumber(ts, &nsteps));
1735:   PetscCall(TSGetConvergedReason(ts, &reason));
1736:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
1737:   PetscCall(TSDestroy(&ts));

1739:   PetscCall(VecTaggerDestroy(&refineTag));
1740:   PetscCall(VecTaggerDestroy(&coarsenTag));
1741:   PetscCall(PetscFunctionListDestroy(&PhysicsList));
1742:   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1743:   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
1744:   PetscCall(PetscFree(user->model->functionalMonitored));
1745:   PetscCall(PetscFree(user->model->functionalCall));
1746:   PetscCall(PetscFree(user->model->physics->data));
1747:   PetscCall(PetscFree(user->model->physics));
1748:   PetscCall(PetscFree(user->model));
1749:   PetscCall(PetscFree(user));
1750:   PetscCall(PetscLimiterDestroy(&limiter));
1751:   PetscCall(PetscLimiterDestroy(&noneLimiter));
1752:   PetscCall(PetscFVDestroy(&fvm));
1753:   PetscCall(DMDestroy(&dm));
1754:   PetscCall(PetscFinalize());
1755:   return 0;
1756: }

1758: /* Godunov fluxs */
1759: PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
1760: {
1761:   /* System generated locals */
1762:   PetscScalar ret_val;

1764:   if (PetscRealPart(*test) > 0.) goto L10;
1765:   ret_val = *b;
1766:   return ret_val;
1767: L10:
1768:   ret_val = *a;
1769:   return ret_val;
1770: } /* cvmgp_ */

1772: PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
1773: {
1774:   /* System generated locals */
1775:   PetscScalar ret_val;

1777:   if (PetscRealPart(*test) < 0.) goto L10;
1778:   ret_val = *b;
1779:   return ret_val;
1780: L10:
1781:   ret_val = *a;
1782:   return ret_val;
1783: } /* cvmgm_ */

1785: 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)
1786: {
1787:   /* Initialized data */

1789:   static PetscScalar smallp = 1e-8;

1791:   /* System generated locals */
1792:   int         i__1;
1793:   PetscScalar d__1, d__2;

1795:   /* Local variables */
1796:   static int         i0;
1797:   static PetscScalar cl, cr, wl, zl, wr, zr, pst, durl, skpr1, skpr2;
1798:   static int         iwave;
1799:   static PetscScalar gascl4, gascr4, cstarl, dpstar, cstarr;
1800:   /* static PetscScalar csqrl, csqrr, gascl1, gascl2, gascl3, gascr1, gascr2, gascr3; */
1801:   static int         iterno;
1802:   static PetscScalar ustarl, ustarr, rarepr1, rarepr2;

1804:   /* gascl1 = *gaml - 1.; */
1805:   /* gascl2 = (*gaml + 1.) * .5; */
1806:   /* gascl3 = gascl2 / *gaml; */
1807:   gascl4 = 1. / (*gaml - 1.);

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

1943:   *ustar = (zl * ustarr + zr * ustarl) / (zl + zr);
1944:   if (PetscRealPart(*pstar) > PetscRealPart(*pl)) {
1945:     pst     = *pstar / *pl;
1946:     *rstarl = ((*gaml + 1.) * pst + *gaml - 1.) / ((*gaml - 1.) * pst + *gaml + 1.) * *rl;
1947:   }
1948:   if (PetscRealPart(*pstar) > PetscRealPart(*pr)) {
1949:     pst     = *pstar / *pr;
1950:     *rstarr = ((*gamr + 1.) * pst + *gamr - 1.) / ((*gamr - 1.) * pst + *gamr + 1.) * *rr;
1951:   }
1952:   return iwave;
1953: }

1955: PetscScalar sign(PetscScalar x)
1956: {
1957:   if (PetscRealPart(x) > 0) return 1.0;
1958:   if (PetscRealPart(x) < 0) return -1.0;
1959:   return 0.0;
1960: }
1961: /*        Riemann Solver */
1962: /* -------------------------------------------------------------------- */
1963: 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)
1964: {
1965:   /* System generated locals */
1966:   PetscScalar d__1, d__2;

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

1973:   if (*rl == *rr && *pr == *pl && *uxl == *uxr && *gaml == *gamr) {
1974:     *rx  = *rl;
1975:     *px  = *pl;
1976:     *uxm = *uxl;
1977:     *gam = *gaml;
1978:     x2   = *xcen + *uxm * *dtt;

1980:     if (PetscRealPart(*xp) >= PetscRealPart(x2)) {
1981:       *utx  = *utr;
1982:       *ubx  = *ubr;
1983:       *rho1 = *rho1r;
1984:     } else {
1985:       *utx  = *utl;
1986:       *ubx  = *ubl;
1987:       *rho1 = *rho1l;
1988:     }
1989:     return 0;
1990:   }
1991:   iwave = riem1mdt(gaml, gamr, rl, pl, uxl, rr, pr, uxr, &rstarl, &rstarr, &pstar, &ustar);

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

2061:   /* Local variables */
2062:   static int         k;
2063:   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;

2065:   /* Function Body */
2066:   xcen = 0.;
2067:   xp   = 0.;
2068:   i__1 = *ndim;
2069:   for (k = 1; k <= i__1; ++k) {
2070:     tg[k - 1] = 0.;
2071:     bn[k - 1] = 0.;
2072:   }
2073:   dtt = 1.;
2074:   if (*ndim == 3) {
2075:     if (nn[0] == 0. && nn[1] == 0.) {
2076:       tg[0] = 1.;
2077:     } else {
2078:       tg[0] = -nn[1];
2079:       tg[1] = nn[0];
2080:     }
2081:     /*           tmp=dsqrt(tg(1)**2+tg(2)**2) */
2082:     /*           tg=tg/tmp */
2083:     bn[0] = -nn[2] * tg[1];
2084:     bn[1] = nn[2] * tg[0];
2085:     bn[2] = nn[0] * tg[1] - nn[1] * tg[0];
2086:     /* Computing 2nd power */
2087:     d__1 = bn[0];
2088:     /* Computing 2nd power */
2089:     d__2 = bn[1];
2090:     /* Computing 2nd power */
2091:     d__3 = bn[2];
2092:     tmp  = PetscSqrtScalar(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
2093:     i__1 = *ndim;
2094:     for (k = 1; k <= i__1; ++k) bn[k - 1] /= tmp;
2095:   } else if (*ndim == 2) {
2096:     tg[0] = -nn[1];
2097:     tg[1] = nn[0];
2098:     /*           tmp=dsqrt(tg(1)**2+tg(2)**2) */
2099:     /*           tg=tg/tmp */
2100:     bn[0] = 0.;
2101:     bn[1] = 0.;
2102:     bn[2] = 1.;
2103:   }
2104:   rl   = ul[0];
2105:   rr   = ur[0];
2106:   uxl  = 0.;
2107:   uxr  = 0.;
2108:   utl  = 0.;
2109:   utr  = 0.;
2110:   ubl  = 0.;
2111:   ubr  = 0.;
2112:   i__1 = *ndim;
2113:   for (k = 1; k <= i__1; ++k) {
2114:     uxl += ul[k] * nn[k - 1];
2115:     uxr += ur[k] * nn[k - 1];
2116:     utl += ul[k] * tg[k - 1];
2117:     utr += ur[k] * tg[k - 1];
2118:     ubl += ul[k] * bn[k - 1];
2119:     ubr += ur[k] * bn[k - 1];
2120:   }
2121:   uxl /= rl;
2122:   uxr /= rr;
2123:   utl /= rl;
2124:   utr /= rr;
2125:   ubl /= rl;
2126:   ubr /= rr;

2128:   gaml = *gamma;
2129:   gamr = *gamma;
2130:   /* Computing 2nd power */
2131:   d__1 = uxl;
2132:   /* Computing 2nd power */
2133:   d__2 = utl;
2134:   /* Computing 2nd power */
2135:   d__3 = ubl;
2136:   pl   = (*gamma - 1.) * (ul[*ndim + 1] - rl * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
2137:   /* Computing 2nd power */
2138:   d__1 = uxr;
2139:   /* Computing 2nd power */
2140:   d__2 = utr;
2141:   /* Computing 2nd power */
2142:   d__3  = ubr;
2143:   pr    = (*gamma - 1.) * (ur[*ndim + 1] - rr * .5 * (d__1 * d__1 + d__2 * d__2 + d__3 * d__3));
2144:   rho1l = rl;
2145:   rho1r = rr;

2147:   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);

2149:   flux[0] = rhom * unm;
2150:   fn      = rhom * unm * unm + pm;
2151:   ft      = rhom * unm * utm;
2152:   /*           flux(2)=fn*nn(1)+ft*nn(2) */
2153:   /*           flux(3)=fn*tg(1)+ft*tg(2) */
2154:   flux[1] = fn * nn[0] + ft * tg[0];
2155:   flux[2] = fn * nn[1] + ft * tg[1];
2156:   /*           flux(2)=rhom*unm*(unm)+pm */
2157:   /*           flux(3)=rhom*(unm)*utm */
2158:   if (*ndim == 3) flux[3] = rhom * unm * ubm;
2159:   flux[*ndim + 1] = (rhom * .5 * (unm * unm + utm * utm + ubm * ubm) + gamm / (gamm - 1.) * pm) * unm;
2160:   return iwave;
2161: } /* godunovflux_ */

2163: /* Subroutine to set up the initial conditions for the */
2164: /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
2165: /* ----------------------------------------------------------------------- */
2166: int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
2167: {
2168:   int j, k;
2169:   /*      Wc=matmul(lv,Ueq) 3 vars */
2170:   for (k = 0; k < 3; ++k) {
2171:     wc[k] = 0.;
2172:     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
2173:   }
2174:   return 0;
2175: }
2176: /* ----------------------------------------------------------------------- */
2177: int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
2178: {
2179:   int k, j;
2180:   /*      V=matmul(rv,WC) 3 vars */
2181:   for (k = 0; k < 3; ++k) {
2182:     v[k] = 0.;
2183:     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
2184:   }
2185:   return 0;
2186: }
2187: /* ---------------------------------------------------------------------- */
2188: int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
2189: {
2190:   int       j, k;
2191:   PetscReal rho, csnd, p0;
2192:   /* PetscScalar u; */

2194:   for (k = 0; k < 3; ++k)
2195:     for (j = 0; j < 3; ++j) {
2196:       lv[k][j] = 0.;
2197:       rv[k][j] = 0.;
2198:     }
2199:   rho = ueq[0];
2200:   /* u = ueq[1]; */
2201:   p0       = ueq[2];
2202:   csnd     = PetscSqrtReal(gamma * p0 / rho);
2203:   lv[0][1] = rho * .5;
2204:   lv[0][2] = -.5 / csnd;
2205:   lv[1][0] = csnd;
2206:   lv[1][2] = -1. / csnd;
2207:   lv[2][1] = rho * .5;
2208:   lv[2][2] = .5 / csnd;
2209:   rv[0][0] = -1. / csnd;
2210:   rv[1][0] = 1. / rho;
2211:   rv[2][0] = -csnd;
2212:   rv[0][1] = 1. / csnd;
2213:   rv[0][2] = 1. / csnd;
2214:   rv[1][2] = 1. / rho;
2215:   rv[2][2] = csnd;
2216:   return 0;
2217: }

2219: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
2220: {
2221:   PetscReal p0, u0, wcp[3], wc[3];
2222:   PetscReal lv[3][3];
2223:   PetscReal vp[3];
2224:   PetscReal rv[3][3];
2225:   PetscReal eps, ueq[3], rho0, twopi;

2227:   /* Function Body */
2228:   twopi  = 2. * PETSC_PI;
2229:   eps    = 1e-4;    /* perturbation */
2230:   rho0   = 1e3;     /* density of water */
2231:   p0     = 101325.; /* init pressure of 1 atm (?) */
2232:   u0     = 0.;
2233:   ueq[0] = rho0;
2234:   ueq[1] = u0;
2235:   ueq[2] = p0;
2236:   /* Project initial state to characteristic variables */
2237:   eigenvectors(rv, lv, ueq, gamma);
2238:   projecteqstate(wc, ueq, lv);
2239:   wcp[0] = wc[0];
2240:   wcp[1] = wc[1];
2241:   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
2242:   projecttoprim(vp, wcp, rv);
2243:   ux->r     = vp[0];         /* density */
2244:   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
2245:   ux->ru[1] = 0.;
2246: #if defined DIM > 2
2247:   if (dim > 2) ux->ru[2] = 0.;
2248: #endif
2249:   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
2250:   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
2251:   return 0;
2252: }

2254: /*TEST

2256:   testset:
2257:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0

2259:     test:
2260:       suffix: adv_2d_tri_0
2261:       requires: triangle
2262:       TODO: how did this ever get in main when there is no support for this
2263:       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

2265:     test:
2266:       suffix: adv_2d_tri_1
2267:       requires: triangle
2268:       TODO: how did this ever get in main when there is no support for this
2269:       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

2271:     test:
2272:       suffix: tut_1
2273:       requires: exodusii
2274:       nsize: 1
2275:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

2277:     test:
2278:       suffix: tut_2
2279:       requires: exodusii
2280:       nsize: 1
2281:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw

2283:     test:
2284:       suffix: tut_3
2285:       requires: exodusii
2286:       nsize: 4
2287:       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

2289:     test:
2290:       suffix: tut_4
2291:       requires: exodusii
2292:       nsize: 4
2293:       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

2295:   testset:
2296:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

2298:     # 2D Advection 0-10
2299:     test:
2300:       suffix: 0
2301:       requires: exodusii
2302:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

2304:     test:
2305:       suffix: 1
2306:       requires: exodusii
2307:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo

2309:     test:
2310:       suffix: 2
2311:       requires: exodusii
2312:       nsize: 2
2313:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

2315:     test:
2316:       suffix: 3
2317:       requires: exodusii
2318:       nsize: 2
2319:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo

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

2327:     test:
2328:       suffix: 5
2329:       requires: exodusii
2330:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1

2332:     test:
2333:       suffix: 7
2334:       requires: exodusii
2335:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1

2337:     test:
2338:       suffix: 8
2339:       requires: exodusii
2340:       nsize: 2
2341:       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

2343:     test:
2344:       suffix: 9
2345:       requires: exodusii
2346:       nsize: 8
2347:       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

2349:     test:
2350:       suffix: 10
2351:       requires: exodusii
2352:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo

2354:   # 2D Shallow water
2355:   testset:
2356:     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0

2358:     test:
2359:       suffix: sw_0
2360:       requires: exodusii
2361:       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
2362:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
2363:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2364:             -monitor height,energy

2366:     test:
2367:       suffix: sw_1
2368:       nsize: 2
2369:       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
2370:             -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 \
2371:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2372:             -monitor height,energy

2374:     test:
2375:       suffix: sw_hll
2376:       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
2377:             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
2378:             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2379:             -monitor height,energy

2381:   testset:
2382:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

2384:     # 2D Advection: p4est
2385:     test:
2386:       suffix: p4est_advec_2d
2387:       requires: p4est
2388:       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

2390:     # Advection in a box
2391:     test:
2392:       suffix: adv_2d_quad_0
2393:       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2395:     test:
2396:       suffix: adv_2d_quad_1
2397:       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
2398:       timeoutfactor: 3

2400:     test:
2401:       suffix: adv_2d_quad_p4est_0
2402:       requires: p4est
2403:       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2405:     test:
2406:       suffix: adv_2d_quad_p4est_1
2407:       requires: p4est
2408:       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
2409:       timeoutfactor: 3

2411:     test:
2412:       suffix: adv_2d_quad_p4est_adapt_0
2413:       requires: p4est !__float128 #broken for quad precision
2414:       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
2415:       timeoutfactor: 3

2417:     test:
2418:       suffix: adv_0
2419:       requires: exodusii
2420:       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

2422:     test:
2423:       suffix: shock_0
2424:       requires: p4est !single !complex
2425:       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
2426:       -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 6 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
2427:       -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
2428:       -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. \
2429:       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
2430:       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2431:       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
2432:       timeoutfactor: 3

2434:     # Test GLVis visualization of PetscFV fields
2435:     test:
2436:       suffix: glvis_adv_2d_tet
2437:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
2438:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
2439:             -ts_monitor_solution glvis: -ts_max_steps 0

2441:     test:
2442:       suffix: glvis_adv_2d_quad
2443:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
2444:             -dm_refine 5 -dm_plex_separate_marker \
2445:             -ts_monitor_solution glvis: -ts_max_steps 0

2447: TEST*/