Actual source code: ex11.c

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

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

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

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

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

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

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

 44: static PetscFunctionList PhysicsList, PhysicsRiemannList_SW;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

590:   PetscFunctionReturn(PETSC_SUCCESS);
591: }

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

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

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

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

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

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

701:   PetscFunctionReturn(PETSC_SUCCESS);
702: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

907:   PetscFunctionReturn(PETSC_SUCCESS);
908: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1337:     PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
1338:     PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
1339:     PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));

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

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

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

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

1414:   PetscFunctionBeginUser;
1415:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1416:   comm = PETSC_COMM_WORLD;

1418:   PetscCall(PetscNew(&user));
1419:   PetscCall(PetscNew(&user->model));
1420:   PetscCall(PetscNew(&user->model->physics));
1421:   mod           = user->model;
1422:   phys          = mod->physics;
1423:   mod->comm     = comm;
1424:   useAMR        = PETSC_FALSE;
1425:   adaptInterval = 1;

1427:   /* Register physical models to be available on the command line */
1428:   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
1429:   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
1430:   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));

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

1449:   if (useAMR) {
1450:     VecTaggerBox refineBox, coarsenBox;

1452:     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1453:     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;

1455:     PetscCall(VecTaggerCreate(comm, &refineTag));
1456:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
1457:     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
1458:     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
1459:     PetscCall(VecTaggerSetFromOptions(refineTag));
1460:     PetscCall(VecTaggerSetUp(refineTag));
1461:     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));

1463:     PetscCall(VecTaggerCreate(comm, &coarsenTag));
1464:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
1465:     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
1466:     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
1467:     PetscCall(VecTaggerSetFromOptions(coarsenTag));
1468:     PetscCall(VecTaggerSetUp(coarsenTag));
1469:     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1470:   }

1472:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1473:   {
1474:     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
1475:     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
1476:     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
1477:     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
1478:     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1479:     /* Count number of fields and dofs */
1480:     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1481:     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1482:     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1483:   }
1484:   PetscOptionsEnd();

1486:   /* Create mesh */
1487:   {
1488:     PetscInt i;

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

1514:         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1515:         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1516:         PetscCall(VecGetLocalSize(coordinates, &nCoords));
1517:         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
1518:         PetscCall(VecGetArray(coordinates, &coords));
1519:         for (i = 0; i < nCoords; i += dimEmbed) {
1520:           PetscInt j;

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

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

1548:   {
1549:     DM gdm;

1551:     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
1552:     PetscCall(DMDestroy(&dm));
1553:     dm = gdm;
1554:     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1555:   }

1557:   PetscCall(PetscFVCreate(comm, &fvm));
1558:   PetscCall(PetscFVSetFromOptions(fvm));
1559:   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
1560:   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
1561:   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1562:   {
1563:     PetscInt f, dof;
1564:     for (f = 0, dof = 0; f < phys->nfields; f++) {
1565:       PetscInt newDof = phys->field_desc[f].dof;

1567:       if (newDof == 1) {
1568:         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
1569:       } else {
1570:         PetscInt j;

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

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

1594:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
1595:     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1596:     PetscOptionsEnd();
1597:     if (flg) {
1598:       DM dmConv;

1600:       PetscCall(DMConvert(dm, convType, &dmConv));
1601:       if (dmConv) {
1602:         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
1603:         PetscCall(DMDestroy(&dm));
1604:         dm = dmConv;
1605:         PetscCall(DMSetFromOptions(dm));
1606:       }
1607:     }
1608:   }

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

1612:   PetscCall(DMCreateGlobalVector(dm, &X));
1613:   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1614:   PetscCall(SetInitialCondition(dm, X, user));
1615:   if (useAMR) {
1616:     PetscInt adaptIter;

1618:     /* use no limiting when reconstructing gradients for adaptivity */
1619:     PetscCall(PetscFVGetLimiter(fvm, &limiter));
1620:     PetscCall(PetscObjectReference((PetscObject)limiter));
1621:     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
1622:     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));

1624:     PetscCall(PetscFVSetLimiter(fvm, noneLimiter));
1625:     for (adaptIter = 0;; ++adaptIter) {
1626:       PetscLogDouble bytes;
1627:       TS             tsNew = NULL;

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

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

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

1681:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1682:   if (vtkCellGeom) {
1683:     DM  dmCell;
1684:     Vec cellgeom, partition;

1686:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1687:     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
1688:     PetscCall(VecView(cellgeom, viewer));
1689:     PetscCall(PetscViewerDestroy(&viewer));
1690:     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
1691:     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
1692:     PetscCall(VecView(partition, viewer));
1693:     PetscCall(PetscViewerDestroy(&viewer));
1694:     PetscCall(VecDestroy(&partition));
1695:     PetscCall(DMDestroy(&dmCell));
1696:   }
1697:   /* collect max maxspeed from all processes -- todo */
1698:   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
1699:   PetscCall(DMDestroy(&plex));
1700:   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1701:   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1702:   dt = cfl * minRadius / mod->maxspeed;
1703:   PetscCall(TSSetTimeStep(ts, dt));
1704:   PetscCall(TSSetFromOptions(ts));
1705:   if (!useAMR) {
1706:     PetscCall(TSSolve(ts, X));
1707:     PetscCall(TSGetSolveTime(ts, &ftime));
1708:     PetscCall(TSGetStepNumber(ts, &nsteps));
1709:   } else {
1710:     PetscReal finalTime;
1711:     PetscInt  adaptIter;
1712:     TS        tsNew  = NULL;
1713:     Vec       solNew = NULL;

1715:     PetscCall(TSGetMaxTime(ts, &finalTime));
1716:     PetscCall(TSSetMaxSteps(ts, adaptInterval));
1717:     PetscCall(TSSolve(ts, X));
1718:     PetscCall(TSGetSolveTime(ts, &ftime));
1719:     PetscCall(TSGetStepNumber(ts, &nsteps));
1720:     for (adaptIter = 0; ftime < finalTime; adaptIter++) {
1721:       PetscLogDouble bytes;

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

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

1781: /* Godunov fluxs */
1782: PetscScalar cvmgp_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
1783: {
1784:   /* System generated locals */
1785:   PetscScalar ret_val;

1787:   if (PetscRealPart(*test) > 0.) goto L10;
1788:   ret_val = *b;
1789:   return ret_val;
1790: L10:
1791:   ret_val = *a;
1792:   return ret_val;
1793: } /* cvmgp_ */

1795: PetscScalar cvmgm_(PetscScalar *a, PetscScalar *b, PetscScalar *test)
1796: {
1797:   /* System generated locals */
1798:   PetscScalar ret_val;

1800:   if (PetscRealPart(*test) < 0.) goto L10;
1801:   ret_val = *b;
1802:   return ret_val;
1803: L10:
1804:   ret_val = *a;
1805:   return ret_val;
1806: } /* cvmgm_ */

1808: 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)
1809: {
1810:   /* Initialized data */

1812:   static PetscScalar smallp = 1e-8;

1814:   /* System generated locals */
1815:   int         i__1;
1816:   PetscScalar d__1, d__2;

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

1827:   /* gascl1 = *gaml - 1.; */
1828:   /* gascl2 = (*gaml + 1.) * .5; */
1829:   /* gascl3 = gascl2 / *gaml; */
1830:   gascl4 = 1. / (*gaml - 1.);

1832:   /* gascr1 = *gamr - 1.; */
1833:   /* gascr2 = (*gamr + 1.) * .5; */
1834:   /* gascr3 = gascr2 / *gamr; */
1835:   gascr4 = 1. / (*gamr - 1.);
1836:   iterno = 10;
1837:   /*        find pstar: */
1838:   cl = PetscSqrtScalar(*gaml * *pl / *rl);
1839:   cr = PetscSqrtScalar(*gamr * *pr / *rr);
1840:   wl = *rl * cl;
1841:   wr = *rr * cr;
1842:   /* csqrl = wl * wl; */
1843:   /* csqrr = wr * wr; */
1844:   *pstar  = (wl * *pr + wr * *pl) / (wl + wr);
1845:   *pstar  = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1846:   pst     = *pl / *pr;
1847:   skpr1   = cr * (pst - 1.) * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1848:   d__1    = (*gamr - 1.) / (*gamr * 2.);
1849:   rarepr2 = gascr4 * 2. * cr * (1. - PetscPowScalar(pst, d__1));
1850:   pst     = *pr / *pl;
1851:   skpr2   = cl * (pst - 1.) * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1852:   d__1    = (*gaml - 1.) / (*gaml * 2.);
1853:   rarepr1 = gascl4 * 2. * cl * (1. - PetscPowScalar(pst, d__1));
1854:   durl    = *uxr - *uxl;
1855:   if (PetscRealPart(*pr) < PetscRealPart(*pl)) {
1856:     if (PetscRealPart(durl) >= PetscRealPart(rarepr1)) {
1857:       iwave = 100;
1858:     } else if (PetscRealPart(durl) <= PetscRealPart(-skpr1)) {
1859:       iwave = 300;
1860:     } else {
1861:       iwave = 400;
1862:     }
1863:   } else {
1864:     if (PetscRealPart(durl) >= PetscRealPart(rarepr2)) {
1865:       iwave = 100;
1866:     } else if (PetscRealPart(durl) <= PetscRealPart(-skpr2)) {
1867:       iwave = 300;
1868:     } else {
1869:       iwave = 200;
1870:     }
1871:   }
1872:   if (iwave == 100) {
1873:     /*     1-wave: rarefaction wave, 3-wave: rarefaction wave */
1874:     /*     case (100) */
1875:     i__1 = iterno;
1876:     for (i0 = 1; i0 <= i__1; ++i0) {
1877:       d__1    = *pstar / *pl;
1878:       d__2    = 1. / *gaml;
1879:       *rstarl = *rl * PetscPowScalar(d__1, d__2);
1880:       cstarl  = PetscSqrtScalar(*gaml * *pstar / *rstarl);
1881:       ustarl  = *uxl - gascl4 * 2. * (cstarl - cl);
1882:       zl      = *rstarl * cstarl;
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:       ustarr  = *uxr + gascr4 * 2. * (cstarr - cr);
1888:       zr      = *rstarr * cstarr;
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: rarefaction wave */
1899:   } else if (iwave == 200) {
1900:     /*     case (200) */
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:       d__1    = *pstar / *pr;
1907:       d__2    = 1. / *gamr;
1908:       *rstarr = *rr * PetscPowScalar(d__1, d__2);
1909:       cstarr  = PetscSqrtScalar(*gamr * *pstar / *rstarr);
1910:       zr      = *rstarr * cstarr;
1911:       ustarr  = *uxr + gascr4 * 2. * (cstarr - cr);
1912:       dpstar  = zl * zr * (ustarr - ustarl) / (zl + zr);
1913:       *pstar -= dpstar;
1914:       *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1915:       if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
1916: #if 0
1917:         break;
1918: #endif
1919:       }
1920:     }
1921:     /*     1-wave: shock wave, 3-wave: shock */
1922:   } else if (iwave == 300) {
1923:     /*     case (300) */
1924:     i__1 = iterno;
1925:     for (i0 = 1; i0 <= i__1; ++i0) {
1926:       pst    = *pstar / *pl;
1927:       ustarl = *uxl - (pst - 1.) * cl * PetscSqrtScalar(2. / (*gaml * (*gaml - 1. + (*gaml + 1.) * pst)));
1928:       zl     = *pl / cl * PetscSqrtScalar(*gaml * 2. * (*gaml - 1. + (*gaml + 1.) * pst)) * (*gaml - 1. + (*gaml + 1.) * pst) / (*gaml * 3. - 1. + (*gaml + 1.) * pst);
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:     /*     1-wave: rarefaction wave, 3-wave: shock */
1942:   } else if (iwave == 400) {
1943:     /*     case (400) */
1944:     i__1 = iterno;
1945:     for (i0 = 1; i0 <= i__1; ++i0) {
1946:       d__1    = *pstar / *pl;
1947:       d__2    = 1. / *gaml;
1948:       *rstarl = *rl * PetscPowScalar(d__1, d__2);
1949:       cstarl  = PetscSqrtScalar(*gaml * *pstar / *rstarl);
1950:       ustarl  = *uxl - gascl4 * 2. * (cstarl - cl);
1951:       zl      = *rstarl * cstarl;
1952:       pst     = *pstar / *pr;
1953:       ustarr  = *uxr + (pst - 1.) * cr * PetscSqrtScalar(2. / (*gamr * (*gamr - 1. + (*gamr + 1.) * pst)));
1954:       zr      = *pr / cr * PetscSqrtScalar(*gamr * 2. * (*gamr - 1. + (*gamr + 1.) * pst)) * (*gamr - 1. + (*gamr + 1.) * pst) / (*gamr * 3. - 1. + (*gamr + 1.) * pst);
1955:       dpstar  = zl * zr * (ustarr - ustarl) / (zl + zr);
1956:       *pstar -= dpstar;
1957:       *pstar = PetscMax(PetscRealPart(*pstar), PetscRealPart(smallp));
1958:       if (PetscAbsScalar(dpstar) / PetscRealPart(*pstar) <= 1e-8) {
1959: #if 0
1960:               break;
1961: #endif
1962:       }
1963:     }
1964:   }

1966:   *ustar = (zl * ustarr + zr * ustarl) / (zl + zr);
1967:   if (PetscRealPart(*pstar) > PetscRealPart(*pl)) {
1968:     pst     = *pstar / *pl;
1969:     *rstarl = ((*gaml + 1.) * pst + *gaml - 1.) / ((*gaml - 1.) * pst + *gaml + 1.) * *rl;
1970:   }
1971:   if (PetscRealPart(*pstar) > PetscRealPart(*pr)) {
1972:     pst     = *pstar / *pr;
1973:     *rstarr = ((*gamr + 1.) * pst + *gamr - 1.) / ((*gamr - 1.) * pst + *gamr + 1.) * *rr;
1974:   }
1975:   return iwave;
1976: }

1978: PetscScalar sign(PetscScalar x)
1979: {
1980:   if (PetscRealPart(x) > 0) return 1.0;
1981:   if (PetscRealPart(x) < 0) return -1.0;
1982:   return 0.0;
1983: }
1984: /*        Riemann Solver */
1985: /* -------------------------------------------------------------------- */
1986: 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)
1987: {
1988:   /* System generated locals */
1989:   PetscScalar d__1, d__2;

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

1996:   if (*rl == *rr && *pr == *pl && *uxl == *uxr && *gaml == *gamr) {
1997:     *rx  = *rl;
1998:     *px  = *pl;
1999:     *uxm = *uxl;
2000:     *gam = *gaml;
2001:     x2   = *xcen + *uxm * *dtt;

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

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

2084:   /* Local variables */
2085:   static int         k;
2086:   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;

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

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

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

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

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

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

2242: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
2243: {
2244:   PetscReal p0, u0, wcp[3], wc[3];
2245:   PetscReal lv[3][3];
2246:   PetscReal vp[3];
2247:   PetscReal rv[3][3];
2248:   PetscReal eps, ueq[3], rho0, twopi;

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

2277: /*TEST

2279:   testset:
2280:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0

2282:     test:
2283:       suffix: adv_2d_tri_0
2284:       requires: triangle
2285:       TODO: how did this ever get in main when there is no support for this
2286:       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

2288:     test:
2289:       suffix: adv_2d_tri_1
2290:       requires: triangle
2291:       TODO: how did this ever get in main when there is no support for this
2292:       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

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

2300:     test:
2301:       suffix: tut_2
2302:       requires: exodusii
2303:       nsize: 1
2304:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw

2306:     test:
2307:       suffix: tut_3
2308:       requires: exodusii
2309:       nsize: 4
2310:       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

2312:     test:
2313:       suffix: tut_4
2314:       requires: exodusii
2315:       nsize: 4
2316:       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

2318:   testset:
2319:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

2321:     # 2D Advection 0-10
2322:     test:
2323:       suffix: 0
2324:       requires: exodusii
2325:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

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

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

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

2344:     test:
2345:       suffix: 4
2346:       requires: exodusii
2347:       nsize: 4
2348:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple

2350:     test:
2351:       suffix: 5
2352:       requires: exodusii
2353:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1

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

2360:     test:
2361:       suffix: 8
2362:       requires: exodusii
2363:       nsize: 2
2364:       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

2366:     test:
2367:       suffix: 9
2368:       requires: exodusii
2369:       nsize: 8
2370:       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

2372:     test:
2373:       suffix: 10
2374:       requires: exodusii
2375:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo

2377:   # 2D Shallow water
2378:   testset:
2379:     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0

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

2389:     test:
2390:       suffix: sw_1
2391:       nsize: 2
2392:       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
2393:             -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 \
2394:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2395:             -monitor height,energy

2397:     test:
2398:       suffix: sw_hll
2399:       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
2400:             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
2401:             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2402:             -monitor height,energy

2404:   testset:
2405:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

2407:     # 2D Advection: p4est
2408:     test:
2409:       suffix: p4est_advec_2d
2410:       requires: p4est
2411:       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

2413:     # Advection in a box
2414:     test:
2415:       suffix: adv_2d_quad_0
2416:       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2418:     test:
2419:       suffix: adv_2d_quad_1
2420:       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
2421:       timeoutfactor: 3

2423:     test:
2424:       suffix: adv_2d_quad_p4est_0
2425:       requires: p4est
2426:       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

2428:     test:
2429:       suffix: adv_2d_quad_p4est_1
2430:       requires: p4est
2431:       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
2432:       timeoutfactor: 3

2434:     test:
2435:       suffix: adv_2d_quad_p4est_adapt_0
2436:       requires: p4est !__float128 #broken for quad precision
2437:       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
2438:       timeoutfactor: 3

2440:     test:
2441:       suffix: adv_0
2442:       requires: exodusii
2443:       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

2445:     test:
2446:       suffix: shock_0
2447:       requires: p4est !single !complex
2448:       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
2449:       -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 6 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
2450:       -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
2451:       -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. \
2452:       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
2453:       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
2454:       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
2455:       timeoutfactor: 3

2457:     # Test GLVis visualization of PetscFV fields
2458:     test:
2459:       suffix: glvis_adv_2d_tet
2460:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
2461:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
2462:             -ts_monitor_solution glvis: -ts_max_steps 0

2464:     test:
2465:       suffix: glvis_adv_2d_quad
2466:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
2467:             -dm_refine 5 -dm_plex_separate_marker \
2468:             -ts_monitor_solution glvis: -ts_max_steps 0

2470: TEST*/