Actual source code: ex11.c

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

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

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

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

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

 35: The mesh is read in from an ExodusII file, usually generated by Cubit.

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

 44: #include "ex11.h"

 46: static PetscFunctionList PhysicsList, PhysicsRiemannList_SW, PhysicsRiemannList_Euler;

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

 58: typedef struct _n_FunctionalLink *FunctionalLink;
 59: struct _n_FunctionalLink {
 60:   char              *name;
 61:   FunctionalFunction func;
 62:   void              *ctx;
 63:   PetscInt           offset;
 64:   FunctionalLink     next;
 65: };

 67: struct _n_Model {
 68:   MPI_Comm         comm; /* Does not do collective communication, but some error conditions can be collective */
 69:   Physics          physics;
 70:   FunctionalLink   functionalRegistry;
 71:   PetscInt         maxComputed;
 72:   PetscInt         numMonitored;
 73:   FunctionalLink  *functionalMonitored;
 74:   PetscInt         numCall;
 75:   FunctionalLink  *functionalCall;
 76:   SolutionFunction solution;
 77:   SetUpBCFunction  setupbc;
 78:   void            *solutionctx;
 79:   PetscReal        maxspeed; /* estimate of global maximum speed (for CFL calculation) */
 80:   PetscReal        bounds[2 * DIM];
 81:   PetscErrorCode (*errorIndicator)(PetscInt, PetscReal, PetscInt, const PetscScalar[], const PetscScalar[], PetscReal *, void *);
 82:   void *errorCtx;
 83:   PetscErrorCode (*setupCEED)(DM, Physics);
 84: };

 86: struct _n_User {
 87:   PetscInt  vtkInterval;                        /* For monitor */
 88:   char      outputBasename[PETSC_MAX_PATH_LEN]; /* Basename for output files */
 89:   PetscInt  monitorStepOffset;
 90:   Model     model;
 91:   PetscBool vtkmon;
 92: };

 94: #ifdef PETSC_HAVE_LIBCEED
 95: // Free a plain data context that was allocated using PETSc; returning libCEED error codes
 96: static int FreeContextPetsc(void *data)
 97: {
 98:   if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed");
 99:   return CEED_ERROR_SUCCESS;
100: }
101: #endif

103: /******************* Advect ********************/
104: typedef enum {
105:   ADVECT_SOL_TILTED,
106:   ADVECT_SOL_BUMP,
107:   ADVECT_SOL_BUMP_CAVITY
108: } AdvectSolType;
109: static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "BUMP_CAVITY", "AdvectSolType", "ADVECT_SOL_", 0};
110: typedef enum {
111:   ADVECT_SOL_BUMP_CONE,
112:   ADVECT_SOL_BUMP_COS
113: } AdvectSolBumpType;
114: static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};

116: typedef struct {
117:   PetscReal wind[DIM];
118: } Physics_Advect_Tilted;
119: typedef struct {
120:   PetscReal         center[DIM];
121:   PetscReal         radius;
122:   AdvectSolBumpType type;
123: } Physics_Advect_Bump;

125: typedef struct {
126:   PetscReal     inflowState;
127:   AdvectSolType soltype;
128:   union
129:   {
130:     Physics_Advect_Tilted tilted;
131:     Physics_Advect_Bump   bump;
132:   } sol;
133:   struct {
134:     PetscInt Solution;
135:     PetscInt Error;
136:   } functional;
137: } Physics_Advect;

139: static const struct FieldDescription PhysicsFields_Advect[] = {
140:   {"U",  1},
141:   {NULL, 0}
142: };

144: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
145: {
146:   Physics         phys   = (Physics)ctx;
147:   Physics_Advect *advect = (Physics_Advect *)phys->data;

149:   PetscFunctionBeginUser;
150:   xG[0] = advect->inflowState;
151:   PetscFunctionReturn(PETSC_SUCCESS);
152: }

154: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
155: {
156:   PetscFunctionBeginUser;
157:   xG[0] = xI[0];
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: 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)
162: {
163:   Physics_Advect *advect = (Physics_Advect *)phys->data;
164:   PetscReal       wind[DIM], wn;

166:   switch (advect->soltype) {
167:   case ADVECT_SOL_TILTED: {
168:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
169:     wind[0]                       = tilted->wind[0];
170:     wind[1]                       = tilted->wind[1];
171:   } break;
172:   case ADVECT_SOL_BUMP:
173:     wind[0] = -qp[1];
174:     wind[1] = qp[0];
175:     break;
176:   case ADVECT_SOL_BUMP_CAVITY: {
177:     PetscInt  i;
178:     PetscReal comp2[3] = {0., 0., 0.}, rad2;

180:     rad2 = 0.;
181:     for (i = 0; i < dim; i++) {
182:       comp2[i] = qp[i] * qp[i];
183:       rad2 += comp2[i];
184:     }

186:     wind[0] = -qp[1];
187:     wind[1] = qp[0];
188:     if (rad2 > 1.) {
189:       PetscInt  maxI     = 0;
190:       PetscReal maxComp2 = comp2[0];

192:       for (i = 1; i < dim; i++) {
193:         if (comp2[i] > maxComp2) {
194:           maxI     = i;
195:           maxComp2 = comp2[i];
196:         }
197:       }
198:       wind[maxI] = 0.;
199:     }
200:   } break;
201:   default: {
202:     PetscInt i;
203:     for (i = 0; i < DIM; ++i) wind[i] = 0.0;
204:   }
205:     /* default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
206:   }
207:   wn      = Dot2Real(wind, n);
208:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
209: }

211: static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
212: {
213:   Physics         phys   = (Physics)ctx;
214:   Physics_Advect *advect = (Physics_Advect *)phys->data;

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

250: static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx)
251: {
252:   Physics         phys      = (Physics)ctx;
253:   Physics_Advect *advect    = (Physics_Advect *)phys->data;
254:   PetscScalar     yexact[1] = {0.0};

256:   PetscFunctionBeginUser;
257:   PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
258:   f[advect->functional.Solution] = PetscRealPart(y[0]);
259:   f[advect->functional.Error]    = PetscAbsScalar(y[0] - yexact[0]);
260:   PetscFunctionReturn(PETSC_SUCCESS);
261: }

263: static PetscErrorCode SetUpBC_Advect(DM dm, PetscDS prob, Physics phys)
264: {
265:   const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
266:   DMLabel        label;

268:   PetscFunctionBeginUser;
269:   /* Register "canned" boundary conditions and defaults for where to apply. */
270:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
271:   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));
272:   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));
273:   PetscFunctionReturn(PETSC_SUCCESS);
274: }

276: static PetscErrorCode PhysicsCreate_Advect(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
277: {
278:   Physics_Advect *advect;

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

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

327: /******************* Shallow Water ********************/
328: static const struct FieldDescription PhysicsFields_SW[] = {
329:   {"Height",   1  },
330:   {"Momentum", DIM},
331:   {NULL,       0  }
332: };

334: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
335: {
336:   PetscFunctionBeginUser;
337:   xG[0] = xI[0];
338:   xG[1] = -xI[1];
339:   xG[2] = -xI[2];
340:   PetscFunctionReturn(PETSC_SUCCESS);
341: }

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

347:   PetscFunctionBeginUser;
348:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
349:   dx[0] = x[0] - 1.5;
350:   dx[1] = x[1] - 1.0;
351:   r     = Norm2Real(dx);
352:   sigma = 0.5;
353:   u[0]  = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma)));
354:   u[1]  = 0.0;
355:   u[2]  = 0.0;
356:   PetscFunctionReturn(PETSC_SUCCESS);
357: }

359: static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
360: {
361:   Physics       phys = (Physics)ctx;
362:   Physics_SW   *sw   = (Physics_SW *)phys->data;
363:   const SWNode *x    = (const SWNode *)xx;
364:   PetscReal     u[2];
365:   PetscReal     h;

367:   PetscFunctionBeginUser;
368:   h = x->h;
369:   Scale2Real(1. / x->h, x->uh, u);
370:   f[sw->functional.Height] = h;
371:   f[sw->functional.Speed]  = Norm2Real(u) + PetscSqrtReal(sw->gravity * h);
372:   f[sw->functional.Energy] = 0.5 * (Dot2Real(x->uh, u) + sw->gravity * PetscSqr(h));
373:   PetscFunctionReturn(PETSC_SUCCESS);
374: }

376: static PetscErrorCode SetUpBC_SW(DM dm, PetscDS prob, Physics phys)
377: {
378:   const PetscInt wallids[] = {100, 101, 200, 300};
379:   DMLabel        label;

381:   PetscFunctionBeginUser;
382:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
383:   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));
384:   PetscFunctionReturn(PETSC_SUCCESS);
385: }

387: #ifdef PETSC_HAVE_LIBCEED
388: static PetscErrorCode CreateQFunctionContext_SW(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
389: {
390:   Physics_SW *in = (Physics_SW *)phys->data;
391:   Physics_SW *sw;

393:   PetscFunctionBeginUser;
394:   PetscCall(PetscCalloc1(1, &sw));

396:   sw->gravity = in->gravity;

398:   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
399:   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*sw), sw));
400:   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
401:   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gravity", offsetof(Physics_SW, gravity), 1, "Acceleration due to gravity"));
402:   PetscFunctionReturn(PETSC_SUCCESS);
403: }
404: #endif

406: static PetscErrorCode SetupCEED_SW(DM dm, Physics physics)
407: {
408: #ifdef PETSC_HAVE_LIBCEED
409:   Ceed                 ceed;
410:   CeedQFunctionContext qfCtx;
411: #endif

413:   PetscFunctionBegin;
414: #ifdef PETSC_HAVE_LIBCEED
415:   PetscCall(DMGetCeed(dm, &ceed));
416:   PetscCall(CreateQFunctionContext_SW(physics, ceed, &qfCtx));
417:   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_SW_Rusanov_CEED, PhysicsRiemann_SW_Rusanov_CEED_loc, qfCtx));
418:   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
419: #endif
420:   PetscFunctionReturn(PETSC_SUCCESS);
421: }

423: static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
424: {
425:   Physics_SW *sw;
426:   char        sw_riemann[64] = "rusanov";

428:   PetscFunctionBeginUser;
429:   phys->field_desc = PhysicsFields_SW;
430:   PetscCall(PetscNew(&sw));
431:   phys->data     = sw;
432:   mod->setupbc   = SetUpBC_SW;
433:   mod->setupCEED = SetupCEED_SW;

435:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov));
436:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL));
437: #ifdef PETSC_HAVE_LIBCEED
438:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov_ceed", PhysicsRiemann_SW_Rusanov_CEED));
439: #endif

441:   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
442:   {
443:     void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
444:     sw->gravity = 1.0;
445:     PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
446:     PetscCall(PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL));
447:     PetscCall(PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW));
448:     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_SW;
449:   }
450:   PetscOptionsHeadEnd();
451:   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */

453:   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
454:   PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
455:   PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
456:   PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
457:   PetscFunctionReturn(PETSC_SUCCESS);
458: }

460: /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
461: /* An initial-value and self-similar solutions of the compressible Euler equations */
462: /* Ravi Samtaney and D. I. Pullin */
463: /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
464: static const struct FieldDescription PhysicsFields_Euler[] = {
465:   {"Density",  1  },
466:   {"Momentum", DIM},
467:   {"Energy",   1  },
468:   {NULL,       0  }
469: };

471: /* initial condition */
472: int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
473: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
474: {
475:   PetscInt       i;
476:   Physics        phys = (Physics)ctx;
477:   Physics_Euler *eu   = (Physics_Euler *)phys->data;
478:   EulerNode     *uu   = (EulerNode *)u;
479:   PetscReal      p0, gamma, c = 0.;

481:   PetscFunctionBeginUser;
482:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);

484:   for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
485:   /* set E and rho */
486:   gamma = eu->gamma;

488:   if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) {
489:     /******************* Euler Density Shock ********************/
490:     /* On initial-value and self-similar solutions of the compressible Euler equations */
491:     /* Ravi Samtaney and D. I. Pullin */
492:     /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
493:     /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity,  */
494:     p0 = 1.;
495:     if (x[0] < 0.0 + x[1] * eu->itana) {
496:       if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */
497:         PetscReal amach, rho, press, gas1, p1;
498:         amach     = eu->amach;
499:         rho       = 1.;
500:         press     = p0;
501:         p1        = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0));
502:         gas1      = (gamma - 1.0) / (gamma + 1.0);
503:         uu->r     = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0);
504:         uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach);
505:         uu->E     = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0];
506:       } else {      /* left of discontinuity (0) */
507:         uu->r = 1.; /* rho = 1 */
508:         uu->E = p0 / (gamma - 1.0);
509:       }
510:     } else { /* right of discontinuity (2) */
511:       uu->r = eu->rhoR;
512:       uu->E = p0 / (gamma - 1.0);
513:     }
514:   } else if (eu->type == EULER_SHOCK_TUBE) {
515:     /* 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. */
516:     if (x[0] < 0.0) {
517:       uu->r = 8.;
518:       uu->E = 10. / (gamma - 1.);
519:     } else {
520:       uu->r = 1.;
521:       uu->E = 1. / (gamma - 1.);
522:     }
523:   } else if (eu->type == EULER_LINEAR_WAVE) {
524:     initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
525:   } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type);

527:   /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */
528:   PetscCall(SpeedOfSound_PG(gamma, uu, &c));
529:   c = (uu->ru[0] / uu->r) + c;
530:   if (c > phys->maxspeed) phys->maxspeed = c;
531:   PetscFunctionReturn(PETSC_SUCCESS);
532: }

534: /* PetscReal* => EulerNode* conversion */
535: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
536: {
537:   PetscInt         i;
538:   const EulerNode *xI   = (const EulerNode *)a_xI;
539:   EulerNode       *xG   = (EulerNode *)a_xG;
540:   Physics          phys = (Physics)ctx;
541:   Physics_Euler   *eu   = (Physics_Euler *)phys->data;

543:   PetscFunctionBeginUser;
544:   xG->r = xI->r;                                     /* ghost cell density - same */
545:   xG->E = xI->E;                                     /* ghost cell energy - same */
546:   if (n[1] != 0.) {                                  /* top and bottom */
547:     xG->ru[0] = xI->ru[0];                           /* copy tang to wall */
548:     xG->ru[1] = -xI->ru[1];                          /* reflect perp to t/b wall */
549:   } else {                                           /* sides */
550:     for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
551:   }
552:   if (eu->type == EULER_LINEAR_WAVE) { /* debug */
553: #if 0
554:     PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]);
555: #endif
556:   }
557:   PetscFunctionReturn(PETSC_SUCCESS);
558: }

560: static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
561: {
562:   Physics          phys = (Physics)ctx;
563:   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
564:   const EulerNode *x    = (const EulerNode *)xx;
565:   PetscReal        p;

567:   PetscFunctionBeginUser;
568:   f[eu->monitor.Density]  = x->r;
569:   f[eu->monitor.Momentum] = NormDIM(x->ru);
570:   f[eu->monitor.Energy]   = x->E;
571:   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
572:   PetscCall(Pressure_PG(eu->gamma, x, &p));
573:   f[eu->monitor.Pressure] = p;
574:   PetscFunctionReturn(PETSC_SUCCESS);
575: }

577: static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys)
578: {
579:   Physics_Euler *eu = (Physics_Euler *)phys->data;
580:   DMLabel        label;

582:   PetscFunctionBeginUser;
583:   PetscCall(DMGetLabel(dm, "Face Sets", &label));
584:   if (eu->type == EULER_LINEAR_WAVE) {
585:     const PetscInt wallids[] = {100, 101};
586:     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));
587:   } else {
588:     const PetscInt wallids[] = {100, 101, 200, 300};
589:     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));
590:   }
591:   PetscFunctionReturn(PETSC_SUCCESS);
592: }

594: #ifdef PETSC_HAVE_LIBCEED
595: static PetscErrorCode CreateQFunctionContext_Euler(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
596: {
597:   Physics_Euler *in = (Physics_Euler *)phys->data;
598:   Physics_Euler *eu;

600:   PetscFunctionBeginUser;
601:   PetscCall(PetscCalloc1(1, &eu));

603:   eu->gamma = in->gamma;

605:   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
606:   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*eu), eu));
607:   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
608:   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gamma", offsetof(Physics_Euler, gamma), 1, "Heat capacity ratio"));
609:   PetscFunctionReturn(PETSC_SUCCESS);
610: }
611: #endif

613: static PetscErrorCode SetupCEED_Euler(DM dm, Physics physics)
614: {
615: #ifdef PETSC_HAVE_LIBCEED
616:   Ceed                 ceed;
617:   CeedQFunctionContext qfCtx;
618: #endif

620:   PetscFunctionBegin;
621: #ifdef PETSC_HAVE_LIBCEED
622:   PetscCall(DMGetCeed(dm, &ceed));
623:   PetscCall(CreateQFunctionContext_Euler(physics, ceed, &qfCtx));
624:   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_Euler_Godunov_CEED, PhysicsRiemann_Euler_Godunov_CEED_loc, qfCtx));
625:   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
626: #endif
627:   PetscFunctionReturn(PETSC_SUCCESS);
628: }

630: static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
631: {
632:   Physics_Euler *eu;

634:   PetscFunctionBeginUser;
635:   phys->field_desc = PhysicsFields_Euler;
636:   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Euler_Godunov;
637:   PetscCall(PetscNew(&eu));
638:   phys->data     = eu;
639:   mod->setupbc   = SetUpBC_Euler;
640:   mod->setupCEED = SetupCEED_Euler;

642:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov", PhysicsRiemann_Euler_Godunov));
643: #ifdef PETSC_HAVE_LIBCEED
644:   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov_ceed", PhysicsRiemann_Euler_Godunov_CEED));
645: #endif

647:   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
648:   {
649:     void (*PhysicsRiemann_Euler)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
650:     PetscReal alpha;
651:     char      type[64]       = "linear_wave";
652:     char      eu_riemann[64] = "godunov";
653:     PetscBool is;
654:     eu->gamma = 1.4;
655:     eu->amach = 2.02;
656:     eu->rhoR  = 3.0;
657:     eu->itana = 0.57735026918963; /* angle of Euler self similar (SS) shock */
658:     PetscCall(PetscOptionsFList("-eu_riemann", "Riemann solver", "", PhysicsRiemannList_Euler, eu_riemann, eu_riemann, sizeof eu_riemann, NULL));
659:     PetscCall(PetscFunctionListFind(PhysicsRiemannList_Euler, eu_riemann, &PhysicsRiemann_Euler));
660:     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Euler;
661:     PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->gamma, &eu->gamma, NULL));
662:     PetscCall(PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->amach, &eu->amach, NULL));
663:     PetscCall(PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->rhoR, &eu->rhoR, NULL));
664:     alpha = 60.;
665:     PetscCall(PetscOptionsRangeReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL, 0.0, 90.0));
666:     eu->itana = 1. / PetscTanReal(alpha * PETSC_PI / 180.0);
667:     PetscCall(PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL));
668:     PetscCall(PetscStrcmp(type, "linear_wave", &is));
669:     if (is) {
670:       /* Remember this should be periodic */
671:       eu->type = EULER_LINEAR_WAVE;
672:       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave"));
673:     } else {
674:       PetscCheck(DIM == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "DIM must be 2 unless linear wave test %s", type);
675:       PetscCall(PetscStrcmp(type, "iv_shock", &is));
676:       if (is) {
677:         eu->type = EULER_IV_SHOCK;
678:         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock"));
679:       } else {
680:         PetscCall(PetscStrcmp(type, "ss_shock", &is));
681:         if (is) {
682:           eu->type = EULER_SS_SHOCK;
683:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock"));
684:         } else {
685:           PetscCall(PetscStrcmp(type, "shock_tube", &is));
686:           if (is) eu->type = EULER_SHOCK_TUBE;
687:           else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type);
688:           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube"));
689:         }
690:       }
691:     }
692:   }
693:   PetscOptionsHeadEnd();
694:   phys->maxspeed = 0.; /* will get set in solution */
695:   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
696:   PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
697:   PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
698:   PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
699:   PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
700:   PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
701:   PetscFunctionReturn(PETSC_SUCCESS);
702: }

704: static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx)
705: {
706:   PetscReal err = 0.;
707:   PetscInt  i, j;

709:   PetscFunctionBeginUser;
710:   for (i = 0; i < numComps; i++) {
711:     for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
712:   }
713:   *error = volume * err;
714:   PetscFunctionReturn(PETSC_SUCCESS);
715: }

717: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
718: {
719:   PetscSF      sfPoint;
720:   PetscSection coordSection;
721:   Vec          coordinates;
722:   PetscSection sectionCell;
723:   PetscScalar *part;
724:   PetscInt     cStart, cEnd, c;
725:   PetscMPIInt  rank;

727:   PetscFunctionBeginUser;
728:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
729:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
730:   PetscCall(DMClone(dm, dmCell));
731:   PetscCall(DMGetPointSF(dm, &sfPoint));
732:   PetscCall(DMSetPointSF(*dmCell, sfPoint));
733:   PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
734:   PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
735:   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
736:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
737:   PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
738:   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
739:   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
740:   PetscCall(PetscSectionSetUp(sectionCell));
741:   PetscCall(DMSetLocalSection(*dmCell, sectionCell));
742:   PetscCall(PetscSectionDestroy(&sectionCell));
743:   PetscCall(DMCreateLocalVector(*dmCell, partition));
744:   PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
745:   PetscCall(VecGetArray(*partition, &part));
746:   for (c = cStart; c < cEnd; ++c) {
747:     PetscScalar *p;

749:     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
750:     p[0] = rank;
751:   }
752:   PetscCall(VecRestoreArray(*partition, &part));
753:   PetscFunctionReturn(PETSC_SUCCESS);
754: }

756: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
757: {
758:   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
759:   PetscSection       coordSection;
760:   Vec                coordinates, facegeom, cellgeom;
761:   PetscSection       sectionMass;
762:   PetscScalar       *m;
763:   const PetscScalar *fgeom, *cgeom, *coords;
764:   PetscInt           vStart, vEnd, v;

766:   PetscFunctionBeginUser;
767:   PetscCall(DMConvert(dm, DMPLEX, &plex));
768:   PetscCall(DMGetCoordinateSection(dm, &coordSection));
769:   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
770:   PetscCall(DMClone(dm, &dmMass));
771:   PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
772:   PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
773:   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
774:   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
775:   PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
776:   for (v = vStart; v < vEnd; ++v) {
777:     PetscInt numFaces;

779:     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
780:     PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
781:   }
782:   PetscCall(PetscSectionSetUp(sectionMass));
783:   PetscCall(DMSetLocalSection(dmMass, sectionMass));
784:   PetscCall(PetscSectionDestroy(&sectionMass));
785:   PetscCall(DMGetLocalVector(dmMass, massMatrix));
786:   PetscCall(VecGetArray(*massMatrix, &m));
787:   PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
788:   PetscCall(VecGetDM(facegeom, &dmFace));
789:   PetscCall(VecGetArrayRead(facegeom, &fgeom));
790:   PetscCall(VecGetDM(cellgeom, &dmCell));
791:   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
792:   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
793:   PetscCall(VecGetArrayRead(coordinates, &coords));
794:   for (v = vStart; v < vEnd; ++v) {
795:     const PetscInt  *faces;
796:     PetscFVFaceGeom *fgA, *fgB, *cg;
797:     PetscScalar     *vertex;
798:     PetscInt         numFaces, sides[2], f, g;

800:     PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
801:     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
802:     PetscCall(DMPlexGetSupport(dmMass, v, &faces));
803:     for (f = 0; f < numFaces; ++f) {
804:       sides[0] = faces[f];
805:       PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
806:       for (g = 0; g < numFaces; ++g) {
807:         const PetscInt *cells = NULL;
808:         PetscReal       area  = 0.0;
809:         PetscInt        numCells;

811:         sides[1] = faces[g];
812:         PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
813:         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
814:         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
815:         PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
816:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
817:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
818:         m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
819:         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
820:       }
821:     }
822:   }
823:   PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
824:   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
825:   PetscCall(VecRestoreArrayRead(coordinates, &coords));
826:   PetscCall(VecRestoreArray(*massMatrix, &m));
827:   PetscCall(DMDestroy(&dmMass));
828:   PetscCall(DMDestroy(&plex));
829:   PetscFunctionReturn(PETSC_SUCCESS);
830: }

832: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
833: static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
834: {
835:   PetscFunctionBeginUser;
836:   mod->solution    = func;
837:   mod->solutionctx = ctx;
838:   PetscFunctionReturn(PETSC_SUCCESS);
839: }

841: static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
842: {
843:   FunctionalLink link, *ptr;
844:   PetscInt       lastoffset = -1;

846:   PetscFunctionBeginUser;
847:   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
848:   PetscCall(PetscNew(&link));
849:   PetscCall(PetscStrallocpy(name, &link->name));
850:   link->offset = lastoffset + 1;
851:   link->func   = func;
852:   link->ctx    = ctx;
853:   link->next   = NULL;
854:   *ptr         = link;
855:   *offset      = link->offset;
856:   PetscFunctionReturn(PETSC_SUCCESS);
857: }

859: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject)
860: {
861:   PetscInt       i, j;
862:   FunctionalLink link;
863:   char          *names[256];

865:   PetscFunctionBeginUser;
866:   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
867:   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
868:   /* Create list of functionals that will be computed somehow */
869:   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
870:   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
871:   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
872:   mod->numCall = 0;
873:   for (i = 0; i < mod->numMonitored; i++) {
874:     for (link = mod->functionalRegistry; link; link = link->next) {
875:       PetscBool match;
876:       PetscCall(PetscStrcasecmp(names[i], link->name, &match));
877:       if (match) break;
878:     }
879:     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
880:     mod->functionalMonitored[i] = link;
881:     for (j = 0; j < i; j++) {
882:       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
883:     }
884:     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
885:   next_name:
886:     PetscCall(PetscFree(names[i]));
887:   }

889:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
890:   mod->maxComputed = -1;
891:   for (link = mod->functionalRegistry; link; link = link->next) {
892:     for (i = 0; i < mod->numCall; i++) {
893:       FunctionalLink call = mod->functionalCall[i];
894:       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
895:     }
896:   }
897:   PetscFunctionReturn(PETSC_SUCCESS);
898: }

900: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
901: {
902:   FunctionalLink l, next;

904:   PetscFunctionBeginUser;
905:   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
906:   l     = *link;
907:   *link = NULL;
908:   for (; l; l = next) {
909:     next = l->next;
910:     PetscCall(PetscFree(l->name));
911:     PetscCall(PetscFree(l));
912:   }
913:   PetscFunctionReturn(PETSC_SUCCESS);
914: }

916: /* put the solution callback into a functional callback */
917: static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
918: {
919:   Model mod;

921:   PetscFunctionBeginUser;
922:   mod = (Model)modctx;
923:   PetscCall((*mod->solution)(mod, time, x, u, mod->solutionctx));
924:   PetscFunctionReturn(PETSC_SUCCESS);
925: }

927: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
928: {
929:   PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
930:   void *ctx[1];
931:   Model mod = user->model;

933:   PetscFunctionBeginUser;
934:   func[0] = SolutionFunctional;
935:   ctx[0]  = (void *)mod;
936:   PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X));
937:   PetscFunctionReturn(PETSC_SUCCESS);
938: }

940: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
941: {
942:   PetscFunctionBeginUser;
943:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
944:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
945:   PetscCall(PetscViewerFileSetName(*viewer, filename));
946:   PetscFunctionReturn(PETSC_SUCCESS);
947: }

949: static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
950: {
951:   User        user = (User)ctx;
952:   DM          dm, plex;
953:   PetscViewer viewer;
954:   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
955:   PetscReal   xnorm;
956:   PetscBool   rollback;

958:   PetscFunctionBeginUser;
959:   PetscCall(TSGetStepRollBack(ts, &rollback));
960:   if (rollback) PetscFunctionReturn(PETSC_SUCCESS);
961:   PetscCall(PetscObjectSetName((PetscObject)X, "u"));
962:   PetscCall(VecGetDM(X, &dm));
963:   PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));

965:   if (stepnum >= 0) stepnum += user->monitorStepOffset;
966:   if (stepnum >= 0) { /* No summary for final time */
967:     Model              mod = user->model;
968:     Vec                cellgeom;
969:     PetscInt           c, cStart, cEnd, fcount, i;
970:     size_t             ftableused, ftablealloc;
971:     const PetscScalar *cgeom, *x;
972:     DM                 dmCell;
973:     DMLabel            vtkLabel;
974:     PetscReal         *fmin, *fmax, *fintegral, *ftmp;

976:     PetscCall(DMConvert(dm, DMPLEX, &plex));
977:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
978:     fcount = mod->maxComputed + 1;
979:     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
980:     for (i = 0; i < fcount; i++) {
981:       fmin[i]      = PETSC_MAX_REAL;
982:       fmax[i]      = PETSC_MIN_REAL;
983:       fintegral[i] = 0;
984:     }
985:     PetscCall(VecGetDM(cellgeom, &dmCell));
986:     PetscCall(DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd));
987:     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
988:     PetscCall(VecGetArrayRead(X, &x));
989:     PetscCall(DMGetLabel(dm, "vtk", &vtkLabel));
990:     for (c = cStart; c < cEnd; ++c) {
991:       PetscFVCellGeom   *cg;
992:       const PetscScalar *cx     = NULL;
993:       PetscInt           vtkVal = 0;

995:       /* not that these two routines as currently implemented work for any dm with a
996:        * localSection/globalSection */
997:       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
998:       PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
999:       if (vtkLabel) PetscCall(DMLabelGetValue(vtkLabel, c, &vtkVal));
1000:       if (!vtkVal || !cx) continue; /* ghost, or not a global cell */
1001:       for (i = 0; i < mod->numCall; i++) {
1002:         FunctionalLink flink = mod->functionalCall[i];
1003:         PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
1004:       }
1005:       for (i = 0; i < fcount; i++) {
1006:         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1007:         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1008:         fintegral[i] += cg->volume * ftmp[i];
1009:       }
1010:     }
1011:     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
1012:     PetscCall(VecRestoreArrayRead(X, &x));
1013:     PetscCall(DMDestroy(&plex));
1014:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmin, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1015:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmax, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1016:     PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fintegral, (PetscMPIInt)fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));

1018:     ftablealloc = fcount * 100;
1019:     ftableused  = 0;
1020:     PetscCall(PetscMalloc1(ftablealloc, &ftable));
1021:     for (i = 0; i < mod->numMonitored; i++) {
1022:       size_t         countused;
1023:       char           buffer[256], *p;
1024:       FunctionalLink flink = mod->functionalMonitored[i];
1025:       PetscInt       id    = flink->offset;
1026:       if (i % 3) {
1027:         PetscCall(PetscArraycpy(buffer, "  ", 2));
1028:         p = buffer + 2;
1029:       } else if (i) {
1030:         char newline[] = "\n";
1031:         PetscCall(PetscMemcpy(buffer, newline, sizeof(newline) - 1));
1032:         p = buffer + sizeof(newline) - 1;
1033:       } else {
1034:         p = buffer;
1035:       }
1036:       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]));
1037:       countused--;
1038:       countused += p - buffer;
1039:       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1040:         char *ftablenew;
1041:         ftablealloc = 2 * ftablealloc + countused;
1042:         PetscCall(PetscMalloc(ftablealloc, &ftablenew));
1043:         PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
1044:         PetscCall(PetscFree(ftable));
1045:         ftable = ftablenew;
1046:       }
1047:       PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
1048:       ftableused += countused;
1049:       ftable[ftableused] = 0;
1050:     }
1051:     PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));

1053:     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
1054:     PetscCall(PetscFree(ftable));
1055:   }
1056:   if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
1057:   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1058:     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
1059:       PetscCall(TSGetStepNumber(ts, &stepnum));
1060:     }
1061:     PetscCall(PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum));
1062:     PetscCall(OutputVTK(dm, filename, &viewer));
1063:     PetscCall(VecView(X, viewer));
1064:     PetscCall(PetscViewerDestroy(&viewer));
1065:   }
1066:   PetscFunctionReturn(PETSC_SUCCESS);
1067: }

1069: static PetscErrorCode initializeTS(DM dm, User user, TS *ts)
1070: {
1071:   PetscBool useCeed;

1073:   PetscFunctionBeginUser;
1074:   PetscCall(TSCreate(PetscObjectComm((PetscObject)dm), ts));
1075:   PetscCall(TSSetType(*ts, TSSSP));
1076:   PetscCall(TSSetDM(*ts, dm));
1077:   if (user->vtkmon) PetscCall(TSMonitorSet(*ts, MonitorVTK, user, NULL));
1078:   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1079:   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user));
1080:   if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user));
1081:   else PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
1082:   PetscCall(TSSetMaxTime(*ts, 2.0));
1083:   PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
1084:   PetscFunctionReturn(PETSC_SUCCESS);
1085: }

1087: typedef struct {
1088:   PetscFV      fvm;
1089:   VecTagger    refineTag;
1090:   VecTagger    coarsenTag;
1091:   DM           adaptedDM;
1092:   User         user;
1093:   PetscReal    cfl;
1094:   PetscLimiter limiter;
1095:   PetscLimiter noneLimiter;
1096: } TransferCtx;

1098: static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
1099: {
1100:   TransferCtx       *tctx       = (TransferCtx *)ctx;
1101:   PetscFV            fvm        = tctx->fvm;
1102:   VecTagger          refineTag  = tctx->refineTag;
1103:   VecTagger          coarsenTag = tctx->coarsenTag;
1104:   User               user       = tctx->user;
1105:   DM                 dm, gradDM, plex, cellDM, adaptedDM = NULL;
1106:   Vec                cellGeom, faceGeom;
1107:   PetscBool          computeGradient;
1108:   Vec                grad, locGrad, locX, errVec;
1109:   PetscInt           cStart, cEnd, c, dim, nRefine, nCoarsen;
1110:   PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1111:   PetscScalar       *errArray;
1112:   const PetscScalar *pointVals;
1113:   const PetscScalar *pointGrads;
1114:   const PetscScalar *pointGeom;
1115:   DMLabel            adaptLabel = NULL;
1116:   IS                 refineIS, coarsenIS;

1118:   PetscFunctionBeginUser;
1119:   *resize = PETSC_FALSE;
1120:   PetscCall(VecGetDM(sol, &dm));
1121:   PetscCall(DMGetDimension(dm, &dim));
1122:   PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
1123:   PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
1124:   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
1125:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1126:   PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
1127:   PetscCall(DMCreateLocalVector(plex, &locX));
1128:   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
1129:   PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
1130:   PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
1131:   PetscCall(DMCreateGlobalVector(gradDM, &grad));
1132:   PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
1133:   PetscCall(DMCreateLocalVector(gradDM, &locGrad));
1134:   PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
1135:   PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
1136:   PetscCall(VecDestroy(&grad));
1137:   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
1138:   PetscCall(VecGetArrayRead(locGrad, &pointGrads));
1139:   PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
1140:   PetscCall(VecGetArrayRead(locX, &pointVals));
1141:   PetscCall(VecGetDM(cellGeom, &cellDM));
1142:   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
1143:   PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
1144:   PetscCall(VecSetUp(errVec));
1145:   PetscCall(VecGetArray(errVec, &errArray));
1146:   for (c = cStart; c < cEnd; c++) {
1147:     PetscReal        errInd = 0.;
1148:     PetscScalar     *pointGrad;
1149:     PetscScalar     *pointVal;
1150:     PetscFVCellGeom *cg;

1152:     PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
1153:     PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
1154:     PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));

1156:     PetscCall((*user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1157:     errArray[c - cStart] = errInd;
1158:     minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
1159:     minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
1160:   }
1161:   PetscCall(VecRestoreArray(errVec, &errArray));
1162:   PetscCall(VecRestoreArrayRead(locX, &pointVals));
1163:   PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
1164:   PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
1165:   PetscCall(VecDestroy(&locGrad));
1166:   PetscCall(VecDestroy(&locX));
1167:   PetscCall(DMDestroy(&plex));

1169:   PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
1170:   PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
1171:   PetscCall(ISGetSize(refineIS, &nRefine));
1172:   PetscCall(ISGetSize(coarsenIS, &nCoarsen));
1173:   if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
1174:   if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
1175:   PetscCall(ISDestroy(&coarsenIS));
1176:   PetscCall(ISDestroy(&refineIS));
1177:   PetscCall(VecDestroy(&errVec));

1179:   PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
1180:   PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1181:   minMaxInd[1] = -minMaxInd[1];
1182:   PetscCallMPI(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1183:   PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1184:   if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
1185:     PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1186:   }
1187:   PetscCall(DMLabelDestroy(&adaptLabel));
1188:   if (adaptedDM) {
1189:     tctx->adaptedDM = adaptedDM;
1190:     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nstep, nRefine, nCoarsen));
1191:     *resize = PETSC_TRUE;
1192:   } else {
1193:     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1194:   }
1195:   PetscFunctionReturn(PETSC_SUCCESS);
1196: }

1198: static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
1199: {
1200:   TransferCtx *tctx = (TransferCtx *)ctx;
1201:   DM           dm;
1202:   PetscReal    time;

1204:   PetscFunctionBeginUser;
1205:   PetscCall(TSGetDM(ts, &dm));
1206:   PetscCall(TSGetTime(ts, &time));
1207:   PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
1208:   for (PetscInt i = 0; i < nv; i++) {
1209:     PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
1210:     PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
1211:   }
1212:   PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */

1214:   Model     mod  = tctx->user->model;
1215:   Physics   phys = mod->physics;
1216:   PetscReal minRadius;

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

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

1225:   PetscCall(TSSetDM(ts, tctx->adaptedDM));
1226:   PetscCall(DMDestroy(&tctx->adaptedDM));
1227:   PetscFunctionReturn(PETSC_SUCCESS);
1228: }

1230: int main(int argc, char **argv)
1231: {
1232:   MPI_Comm          comm;
1233:   PetscDS           prob;
1234:   PetscFV           fvm;
1235:   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1236:   User              user;
1237:   Model             mod;
1238:   Physics           phys;
1239:   DM                dm, plex;
1240:   PetscReal         ftime, cfl, dt, minRadius;
1241:   PetscInt          dim, nsteps;
1242:   TS                ts;
1243:   TSConvergedReason reason;
1244:   Vec               X;
1245:   PetscViewer       viewer;
1246:   PetscBool         vtkCellGeom, useAMR;
1247:   PetscInt          adaptInterval;
1248:   char              physname[256] = "advect";
1249:   VecTagger         refineTag = NULL, coarsenTag = NULL;
1250:   TransferCtx       tctx;

1252:   PetscFunctionBeginUser;
1253:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1254:   comm = PETSC_COMM_WORLD;

1256:   PetscCall(PetscNew(&user));
1257:   PetscCall(PetscNew(&user->model));
1258:   PetscCall(PetscNew(&user->model->physics));
1259:   mod           = user->model;
1260:   phys          = mod->physics;
1261:   mod->comm     = comm;
1262:   useAMR        = PETSC_FALSE;
1263:   adaptInterval = 1;

1265:   /* Register physical models to be available on the command line */
1266:   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
1267:   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
1268:   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));

1270:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1271:   {
1272:     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1273:     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1274:     user->vtkInterval = 1;
1275:     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1276:     user->vtkmon = PETSC_TRUE;
1277:     PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1278:     vtkCellGeom = PETSC_FALSE;
1279:     PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
1280:     PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
1281:     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
1282:     PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
1283:     PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1284:   }
1285:   PetscOptionsEnd();

1287:   if (useAMR) {
1288:     VecTaggerBox refineBox, coarsenBox;

1290:     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1291:     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;

1293:     PetscCall(VecTaggerCreate(comm, &refineTag));
1294:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
1295:     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
1296:     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
1297:     PetscCall(VecTaggerSetFromOptions(refineTag));
1298:     PetscCall(VecTaggerSetUp(refineTag));
1299:     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));

1301:     PetscCall(VecTaggerCreate(comm, &coarsenTag));
1302:     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
1303:     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
1304:     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
1305:     PetscCall(VecTaggerSetFromOptions(coarsenTag));
1306:     PetscCall(VecTaggerSetUp(coarsenTag));
1307:     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1308:   }

1310:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1311:   {
1312:     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
1313:     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
1314:     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
1315:     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
1316:     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1317:     /* Count number of fields and dofs */
1318:     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1319:     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1320:     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1321:   }
1322:   PetscOptionsEnd();

1324:   /* Create mesh */
1325:   {
1326:     PetscInt i;

1328:     PetscCall(DMCreate(comm, &dm));
1329:     PetscCall(DMSetType(dm, DMPLEX));
1330:     PetscCall(DMSetFromOptions(dm));
1331:     for (i = 0; i < DIM; i++) {
1332:       mod->bounds[2 * i]     = 0.;
1333:       mod->bounds[2 * i + 1] = 1.;
1334:     }
1335:     dim = DIM;
1336:     { /* a null name means just do a hex box */
1337:       PetscInt  cells[3] = {1, 1, 1}, n = 3;
1338:       PetscBool flg2, skew              = PETSC_FALSE;
1339:       PetscInt  nret2 = 2 * DIM;
1340:       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
1341:       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));
1342:       PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
1343:       PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1344:       PetscOptionsEnd();
1345:       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1346:       if (flg2) {
1347:         PetscInt     dimEmbed, i;
1348:         PetscInt     nCoords;
1349:         PetscScalar *coords;
1350:         Vec          coordinates;

1352:         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1353:         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1354:         PetscCall(VecGetLocalSize(coordinates, &nCoords));
1355:         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
1356:         PetscCall(VecGetArray(coordinates, &coords));
1357:         for (i = 0; i < nCoords; i += dimEmbed) {
1358:           PetscInt j;

1360:           PetscScalar *coord = &coords[i];
1361:           for (j = 0; j < dimEmbed; j++) {
1362:             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1363:             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1364:               if (cells[0] == 2 && i == 8) {
1365:                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
1366:               } else if (cells[0] == 3) {
1367:                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1368:                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1369:                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1370:               }
1371:             }
1372:           }
1373:         }
1374:         PetscCall(VecRestoreArray(coordinates, &coords));
1375:         PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1376:       }
1377:     }
1378:   }
1379:   PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
1380:   PetscCall(DMGetDimension(dm, &dim));

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

1386:   {
1387:     DM gdm;

1389:     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
1390:     PetscCall(DMDestroy(&dm));
1391:     dm = gdm;
1392:     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1393:   }

1395:   PetscCall(PetscFVCreate(comm, &fvm));
1396:   PetscCall(PetscFVSetFromOptions(fvm));
1397:   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
1398:   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
1399:   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1400:   {
1401:     PetscInt f, dof;
1402:     for (f = 0, dof = 0; f < phys->nfields; f++) {
1403:       PetscInt newDof = phys->field_desc[f].dof;

1405:       if (newDof == 1) {
1406:         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
1407:       } else {
1408:         PetscInt j;

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

1413:           PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
1414:           PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1415:         }
1416:       }
1417:       dof += newDof;
1418:     }
1419:   }
1420:   /* FV is now structured with one field having all physics as components */
1421:   PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
1422:   PetscCall(DMCreateDS(dm));
1423:   PetscCall(DMGetDS(dm, &prob));
1424:   PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
1425:   PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
1426:   PetscCall((*mod->setupbc)(dm, prob, phys));
1427:   PetscCall(PetscDSSetFromOptions(prob));
1428:   {
1429:     char      convType[256];
1430:     PetscBool flg;

1432:     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
1433:     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1434:     PetscOptionsEnd();
1435:     if (flg) {
1436:       DM dmConv;

1438:       PetscCall(DMConvert(dm, convType, &dmConv));
1439:       if (dmConv) {
1440:         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
1441:         PetscCall(DMDestroy(&dm));
1442:         dm = dmConv;
1443:         PetscCall(DMSetFromOptions(dm));
1444:       }
1445:     }
1446:   }
1447: #ifdef PETSC_HAVE_LIBCEED
1448:   {
1449:     PetscBool useCeed;
1450:     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1451:     if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1452:   }
1453: #endif

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

1457:   PetscCall(DMCreateGlobalVector(dm, &X));
1458:   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1459:   PetscCall(SetInitialCondition(dm, X, user));
1460:   if (useAMR) {
1461:     PetscInt adaptIter;

1463:     /* use no limiting when reconstructing gradients for adaptivity */
1464:     PetscCall(PetscFVGetLimiter(fvm, &limiter));
1465:     PetscCall(PetscObjectReference((PetscObject)limiter));
1466:     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
1467:     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));

1469:     /* Refinement context */
1470:     tctx.fvm         = fvm;
1471:     tctx.refineTag   = refineTag;
1472:     tctx.coarsenTag  = coarsenTag;
1473:     tctx.adaptedDM   = NULL;
1474:     tctx.user        = user;
1475:     tctx.noneLimiter = noneLimiter;
1476:     tctx.limiter     = limiter;
1477:     tctx.cfl         = cfl;

1479:     /* Do some initial refinement steps */
1480:     for (adaptIter = 0;; ++adaptIter) {
1481:       PetscLogDouble bytes;
1482:       PetscBool      resize;

1484:       PetscCall(PetscMemoryGetCurrentUsage(&bytes));
1485:       PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
1486:       PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
1487:       PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));

1489:       PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
1490:       if (!resize) break;
1491:       PetscCall(DMDestroy(&dm));
1492:       PetscCall(VecDestroy(&X));
1493:       dm             = tctx.adaptedDM;
1494:       tctx.adaptedDM = NULL;
1495:       PetscCall(TSSetDM(ts, dm));
1496:       PetscCall(DMCreateGlobalVector(dm, &X));
1497:       PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1498:       PetscCall(SetInitialCondition(dm, X, user));
1499:     }
1500:   }

1502:   PetscCall(DMConvert(dm, DMPLEX, &plex));
1503:   if (vtkCellGeom) {
1504:     DM  dmCell;
1505:     Vec cellgeom, partition;

1507:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1508:     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
1509:     PetscCall(VecView(cellgeom, viewer));
1510:     PetscCall(PetscViewerDestroy(&viewer));
1511:     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
1512:     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
1513:     PetscCall(VecView(partition, viewer));
1514:     PetscCall(PetscViewerDestroy(&viewer));
1515:     PetscCall(VecDestroy(&partition));
1516:     PetscCall(DMDestroy(&dmCell));
1517:   }
1518:   /* collect max maxspeed from all processes -- todo */
1519:   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
1520:   PetscCall(DMDestroy(&plex));
1521:   PetscCallMPI(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1522:   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1523:   dt = cfl * minRadius / mod->maxspeed;
1524:   PetscCall(TSSetTimeStep(ts, dt));
1525:   PetscCall(TSSetFromOptions(ts));

1527:   /* When using adaptive mesh refinement
1528:      specify callbacks to refine the solution
1529:      and interpolate data from old to new mesh
1530:      When mesh adaption is requested, the step will be rolled back
1531:   */
1532:   if (useAMR) PetscCall(TSSetResize(ts, PETSC_TRUE, adaptToleranceFVMSetUp, Transfer, &tctx));
1533:   PetscCall(TSSetSolution(ts, X));
1534:   PetscCall(VecDestroy(&X));
1535:   PetscCall(TSSolve(ts, NULL));
1536:   PetscCall(TSGetSolveTime(ts, &ftime));
1537:   PetscCall(TSGetStepNumber(ts, &nsteps));
1538:   PetscCall(TSGetConvergedReason(ts, &reason));
1539:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
1540:   PetscCall(TSDestroy(&ts));

1542:   PetscCall(VecTaggerDestroy(&refineTag));
1543:   PetscCall(VecTaggerDestroy(&coarsenTag));
1544:   PetscCall(PetscFunctionListDestroy(&PhysicsList));
1545:   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1546:   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
1547:   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
1548:   PetscCall(PetscFree(user->model->functionalMonitored));
1549:   PetscCall(PetscFree(user->model->functionalCall));
1550:   PetscCall(PetscFree(user->model->physics->data));
1551:   PetscCall(PetscFree(user->model->physics));
1552:   PetscCall(PetscFree(user->model));
1553:   PetscCall(PetscFree(user));
1554:   PetscCall(PetscLimiterDestroy(&limiter));
1555:   PetscCall(PetscLimiterDestroy(&noneLimiter));
1556:   PetscCall(PetscFVDestroy(&fvm));
1557:   PetscCall(DMDestroy(&dm));
1558:   PetscCall(PetscFinalize());
1559:   return 0;
1560: }

1562: /* Subroutine to set up the initial conditions for the */
1563: /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1564: /* ----------------------------------------------------------------------- */
1565: int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1566: {
1567:   int j, k;
1568:   /*      Wc=matmul(lv,Ueq) 3 vars */
1569:   for (k = 0; k < 3; ++k) {
1570:     wc[k] = 0.;
1571:     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1572:   }
1573:   return 0;
1574: }
1575: /* ----------------------------------------------------------------------- */
1576: int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1577: {
1578:   int k, j;
1579:   /*      V=matmul(rv,WC) 3 vars */
1580:   for (k = 0; k < 3; ++k) {
1581:     v[k] = 0.;
1582:     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1583:   }
1584:   return 0;
1585: }
1586: /* ---------------------------------------------------------------------- */
1587: int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1588: {
1589:   int       j, k;
1590:   PetscReal rho, csnd, p0;
1591:   /* PetscScalar u; */

1593:   for (k = 0; k < 3; ++k)
1594:     for (j = 0; j < 3; ++j) {
1595:       lv[k][j] = 0.;
1596:       rv[k][j] = 0.;
1597:     }
1598:   rho = ueq[0];
1599:   /* u = ueq[1]; */
1600:   p0       = ueq[2];
1601:   csnd     = PetscSqrtReal(gamma * p0 / rho);
1602:   lv[0][1] = rho * .5;
1603:   lv[0][2] = -.5 / csnd;
1604:   lv[1][0] = csnd;
1605:   lv[1][2] = -1. / csnd;
1606:   lv[2][1] = rho * .5;
1607:   lv[2][2] = .5 / csnd;
1608:   rv[0][0] = -1. / csnd;
1609:   rv[1][0] = 1. / rho;
1610:   rv[2][0] = -csnd;
1611:   rv[0][1] = 1. / csnd;
1612:   rv[0][2] = 1. / csnd;
1613:   rv[1][2] = 1. / rho;
1614:   rv[2][2] = csnd;
1615:   return 0;
1616: }

1618: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1619: {
1620:   PetscReal p0, u0, wcp[3], wc[3];
1621:   PetscReal lv[3][3];
1622:   PetscReal vp[3];
1623:   PetscReal rv[3][3];
1624:   PetscReal eps, ueq[3], rho0, twopi;

1626:   /* Function Body */
1627:   twopi  = 2. * PETSC_PI;
1628:   eps    = 1e-4;    /* perturbation */
1629:   rho0   = 1e3;     /* density of water */
1630:   p0     = 101325.; /* init pressure of 1 atm (?) */
1631:   u0     = 0.;
1632:   ueq[0] = rho0;
1633:   ueq[1] = u0;
1634:   ueq[2] = p0;
1635:   /* Project initial state to characteristic variables */
1636:   eigenvectors(rv, lv, ueq, gamma);
1637:   projecteqstate(wc, ueq, lv);
1638:   wcp[0] = wc[0];
1639:   wcp[1] = wc[1];
1640:   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1641:   projecttoprim(vp, wcp, rv);
1642:   ux->r     = vp[0];         /* density */
1643:   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1644:   ux->ru[1] = 0.;
1645: #if defined DIM > 2
1646:   if (dim > 2) ux->ru[2] = 0.;
1647: #endif
1648:   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1649:   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1650:   return 0;
1651: }

1653: /*TEST

1655:   testset:
1656:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0

1658:     test:
1659:       suffix: adv_2d_tri_0
1660:       requires: triangle
1661:       TODO: how did this ever get in main when there is no support for this
1662:       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

1664:     test:
1665:       suffix: adv_2d_tri_1
1666:       requires: triangle
1667:       TODO: how did this ever get in main when there is no support for this
1668:       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

1670:     test:
1671:       suffix: tut_1
1672:       requires: exodusii
1673:       nsize: 1
1674:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -

1676:     test:
1677:       suffix: tut_2
1678:       requires: exodusii
1679:       nsize: 1
1680:       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw

1682:     test:
1683:       suffix: tut_3
1684:       requires: exodusii
1685:       nsize: 4
1686:       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

1688:     test:
1689:       suffix: tut_4
1690:       requires: exodusii !single
1691:       nsize: 4
1692:       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

1694:   testset:
1695:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

1697:     # 2D Advection 0-10
1698:     test:
1699:       suffix: 0
1700:       requires: exodusii
1701:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

1703:     test:
1704:       suffix: 1
1705:       requires: exodusii
1706:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo

1708:     test:
1709:       suffix: 2
1710:       requires: exodusii
1711:       nsize: 2
1712:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo

1714:     test:
1715:       suffix: 3
1716:       requires: exodusii
1717:       nsize: 2
1718:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo

1720:     test:
1721:       suffix: 4
1722:       requires: exodusii
1723:       nsize: 4
1724:       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple

1726:     test:
1727:       suffix: 5
1728:       requires: exodusii
1729:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1

1731:     test:
1732:       suffix: 7
1733:       requires: exodusii
1734:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1

1736:     test:
1737:       suffix: 8
1738:       requires: exodusii
1739:       nsize: 2
1740:       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

1742:     test:
1743:       suffix: 9
1744:       requires: exodusii
1745:       nsize: 8
1746:       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

1748:     test:
1749:       suffix: 10
1750:       requires: exodusii
1751:       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo

1753:   # 2D Shallow water
1754:   testset:
1755:     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0

1757:     test:
1758:       suffix: sw_0
1759:       requires: exodusii
1760:       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1761:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1762:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1763:             -monitor height,energy

1765:     test:
1766:       suffix: sw_ceed
1767:       requires: exodusii libceed
1768:       args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1769:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1770:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1771:             -monitor height,energy

1773:     test:
1774:       suffix: sw_ceed_small
1775:       requires: exodusii libceed
1776:       args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \
1777:             -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_lower 0,1 -dm_plex_box_upper 6.28,3 -dm_plex_box_faces 8,2 \
1778:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1779:             -monitor height,energy

1781:     test:
1782:       suffix: sw_1
1783:       nsize: 2
1784:       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1785:             -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 \
1786:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1787:             -monitor height,energy

1789:     test:
1790:       suffix: sw_hll
1791:       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1792:             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1793:             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1794:             -monitor height,energy

1796:   # 2D Euler
1797:   testset:
1798:     args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1799:           -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy

1801:     test:
1802:       suffix: euler_0
1803:       requires: exodusii !complex !single
1804:       args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1805:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1806:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10

1808:     test:
1809:       suffix: euler_ceed
1810:       requires: exodusii libceed
1811:       args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1812:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1813:             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10

1815:   testset:
1816:     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1

1818:     # 2D Advection: p4est
1819:     test:
1820:       suffix: p4est_advec_2d
1821:       requires: p4est
1822:       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

1824:     # Advection in a box
1825:     test:
1826:       suffix: adv_2d_quad_0
1827:       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

1829:     test:
1830:       suffix: adv_2d_quad_1
1831:       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
1832:       timeoutfactor: 3

1834:     test:
1835:       suffix: adv_2d_quad_p4est_0
1836:       requires: p4est
1837:       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3

1839:     test:
1840:       suffix: adv_2d_quad_p4est_1
1841:       requires: p4est
1842:       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
1843:       timeoutfactor: 3

1845:     test: # broken for quad precision
1846:       suffix: adv_2d_quad_p4est_adapt_0
1847:       requires: p4est !__float128
1848:       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
1849:       timeoutfactor: 3

1851:     test:
1852:       suffix: adv_0
1853:       requires: exodusii
1854:       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

1856:     # Run with -dm_forest_maximum_refinement 6 -ts_max_time 0.5 instead to get the full movie
1857:     test:
1858:       suffix: shock_0
1859:       requires: p4est !single !complex
1860:       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
1861:       -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 2 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
1862:       -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
1863:       -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. \
1864:       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
1865:       -ts_max_steps 3 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1866:       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy

1868:     # Test GLVis visualization of PetscFV fields
1869:     test:
1870:       suffix: glvis_adv_2d_tri
1871:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1872:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1873:             -ts_monitor_solution glvis: -ts_max_steps 0

1875:     test:
1876:       suffix: glvis_adv_2d_quad
1877:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
1878:             -dm_refine 5 -dm_plex_separate_marker \
1879:             -ts_monitor_solution glvis: -ts_max_steps 0

1881:     # Test CGNS file writing for PetscFV fields
1882:     test:
1883:       suffix: cgns_adv_2d_tri
1884:       requires: cgns
1885:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1886:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1887:             -ts_monitor_solution cgns:sol.cgns -ts_max_steps 0

1889:     test:
1890:       suffix: cgns_adv_2d_quad
1891:       requires: cgns
1892:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
1893:             -dm_refine 5 -dm_plex_separate_marker \
1894:             -ts_monitor_solution cgns:sol.cgns -ts_max_steps 0

1896:     # Test VTK file writing for PetscFV fields with -ts_monitor_solution_vtk
1897:     test:
1898:       suffix: vtk_adv_2d_tri
1899:       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1900:             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1901:             -ts_monitor_solution_vtk 'bar-%03d.vtu'  -ts_monitor_solution_vtk_interval  2 -ts_max_steps 5

1903: TEST*/