Actual source code: ex11_sa.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 <petscts.h>
 38: #include <petscfv.h>
 39: #include <petscdmplex.h>
 40: #include <petscsf.h>
 41: #include <petscblaslapack.h>

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

 45: static PetscFunctionList PhysicsList;

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

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

 54: /* 'User' implements a discretization of a continuous model. */
 55: typedef struct _n_User *User;

 57: typedef PetscErrorCode (*RiemannFunction)(const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscScalar *, void *);
 58: typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
 59: typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
 60: typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
 61: static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
 62: static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
 63: static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);

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

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

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

 88: struct _n_Model {
 89:   MPI_Comm         comm; /* Does not do collective communicaton, but some error conditions can be collective */
 90:   Physics          physics;
 91:   FunctionalLink   functionalRegistry;
 92:   PetscInt         maxComputed;
 93:   PetscInt         numMonitored;
 94:   FunctionalLink  *functionalMonitored;
 95:   PetscInt         numCall;
 96:   FunctionalLink  *functionalCall;
 97:   SolutionFunction solution;
 98:   void            *solutionctx;
 99:   PetscReal        maxspeed; /* estimate of global maximum speed (for CFL calculation) */
100: };

102: struct _n_User {
103:   PetscInt numSplitFaces;
104:   PetscInt vtkInterval; /* For monitor */
105:   Model    model;
106: };

108: static inline PetscScalar DotDIM(const PetscScalar *x, const PetscScalar *y)
109: {
110:   PetscInt    i;
111:   PetscScalar prod = 0.0;

113:   for (i = 0; i < DIM; i++) prod += x[i] * y[i];
114:   return prod;
115: }
116: static inline PetscReal NormDIM(const PetscScalar *x)
117: {
118:   return PetscSqrtReal(PetscAbsScalar(DotDIM(x, x)));
119: }
120: static inline void axDIM(const PetscScalar a, PetscScalar *x)
121: {
122:   PetscInt i;
123:   for (i = 0; i < DIM; i++) x[i] *= a;
124: }
125: static inline void waxDIM(const PetscScalar a, const PetscScalar *x, PetscScalar *w)
126: {
127:   PetscInt i;
128:   for (i = 0; i < DIM; i++) w[i] = x[i] * a;
129: }
130: static inline void NormalSplitDIM(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
131: { /* Split x into normal and tangential components */
132:   PetscInt    i;
133:   PetscScalar c;
134:   c = DotDIM(x, n) / DotDIM(n, n);
135:   for (i = 0; i < DIM; i++) {
136:     xn[i] = c * n[i];
137:     xt[i] = x[i] - xn[i];
138:   }
139: }

141: static inline PetscScalar Dot2(const PetscScalar *x, const PetscScalar *y)
142: {
143:   return x[0] * y[0] + x[1] * y[1];
144: }
145: static inline PetscReal Norm2(const PetscScalar *x)
146: {
147:   return PetscSqrtReal(PetscAbsScalar(Dot2(x, x)));
148: }
149: static inline void Normalize2(PetscScalar *x)
150: {
151:   PetscReal a = 1. / Norm2(x);
152:   x[0] *= a;
153:   x[1] *= a;
154: }
155: static inline void Waxpy2(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
156: {
157:   w[0] = a * x[0] + y[0];
158:   w[1] = a * x[1] + y[1];
159: }
160: static inline void Scale2(PetscScalar a, const PetscScalar *x, PetscScalar *y)
161: {
162:   y[0] = a * x[0];
163:   y[1] = a * x[1];
164: }

166: static inline void WaxpyD(PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
167: {
168:   PetscInt d;
169:   for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
170: }
171: static inline PetscScalar DotD(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
172: {
173:   PetscScalar sum = 0.0;
174:   PetscInt    d;
175:   for (d = 0; d < dim; ++d) sum += x[d] * y[d];
176:   return sum;
177: }
178: static inline PetscReal NormD(PetscInt dim, const PetscScalar *x)
179: {
180:   return PetscSqrtReal(PetscAbsScalar(DotD(dim, x, x)));
181: }

183: static inline void NormalSplit(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
184: { /* Split x into normal and tangential components */
185:   Scale2(Dot2(x, n) / Dot2(n, n), n, xn);
186:   Waxpy2(-1, xn, x, xt);
187: }

189: /******************* Advect ********************/
190: typedef enum {
191:   ADVECT_SOL_TILTED,
192:   ADVECT_SOL_BUMP
193: } AdvectSolType;
194: static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "AdvectSolType", "ADVECT_SOL_", 0};
195: typedef enum {
196:   ADVECT_SOL_BUMP_CONE,
197:   ADVECT_SOL_BUMP_COS
198: } AdvectSolBumpType;
199: static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};

201: typedef struct {
202:   PetscReal wind[DIM];
203: } Physics_Advect_Tilted;
204: typedef struct {
205:   PetscReal         center[DIM];
206:   PetscReal         radius;
207:   AdvectSolBumpType type;
208: } Physics_Advect_Bump;

210: typedef struct {
211:   PetscReal     inflowState;
212:   AdvectSolType soltype;
213:   union
214:   {
215:     Physics_Advect_Tilted tilted;
216:     Physics_Advect_Bump   bump;
217:   } sol;
218:   struct {
219:     PetscInt Error;
220:   } functional;
221: } Physics_Advect;

223: static const struct FieldDescription PhysicsFields_Advect[] = {
224:   {"U",  1},
225:   {NULL, 0}
226: };

228: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
229: {
230:   Physics         phys   = (Physics)ctx;
231:   Physics_Advect *advect = (Physics_Advect *)phys->data;

234:   xG[0] = advect->inflowState;
235:   return 0;
236: }

238: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
239: {
241:   xG[0] = xI[0];
242:   return 0;
243: }

245: static PetscErrorCode PhysicsRiemann_Advect(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
246: {
247:   Physics_Advect *advect = (Physics_Advect *)phys->data;
248:   PetscReal       wind[DIM], wn;

251:   switch (advect->soltype) {
252:   case ADVECT_SOL_TILTED: {
253:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
254:     wind[0]                       = tilted->wind[0];
255:     wind[1]                       = tilted->wind[1];
256:   } break;
257:   case ADVECT_SOL_BUMP:
258:     wind[0] = -qp[1];
259:     wind[1] = qp[0];
260:     break;
261:   default:
262:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for solution type %s", AdvectSolBumpTypes[advect->soltype]);
263:   }
264:   wn      = Dot2(wind, n);
265:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
266:   return 0;
267: }

269: static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
270: {
271:   Physics         phys   = (Physics)ctx;
272:   Physics_Advect *advect = (Physics_Advect *)phys->data;

275:   switch (advect->soltype) {
276:   case ADVECT_SOL_TILTED: {
277:     PetscReal              x0[DIM];
278:     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
279:     Waxpy2(-time, tilted->wind, x, x0);
280:     if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
281:     else u[0] = advect->inflowState;
282:   } break;
283:   case ADVECT_SOL_BUMP: {
284:     Physics_Advect_Bump *bump = &advect->sol.bump;
285:     PetscReal            x0[DIM], v[DIM], r, cost, sint;
286:     cost  = PetscCosReal(time);
287:     sint  = PetscSinReal(time);
288:     x0[0] = cost * x[0] + sint * x[1];
289:     x0[1] = -sint * x[0] + cost * x[1];
290:     Waxpy2(-1, bump->center, x0, v);
291:     r = Norm2(v);
292:     switch (bump->type) {
293:     case ADVECT_SOL_BUMP_CONE:
294:       u[0] = PetscMax(1 - r / bump->radius, 0);
295:       break;
296:     case ADVECT_SOL_BUMP_COS:
297:       u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
298:       break;
299:     }
300:   } break;
301:   default:
302:     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
303:   }
304:   return 0;
305: }

307: static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscScalar *x, const PetscScalar *y, PetscReal *f, void *ctx)
308: {
309:   Physics         phys   = (Physics)ctx;
310:   Physics_Advect *advect = (Physics_Advect *)phys->data;
311:   PetscScalar     yexact[1];

314:   PhysicsSolution_Advect(mod, time, x, yexact, phys);
315:   f[advect->functional.Error] = PetscAbsScalar(y[0] - yexact[0]);
316:   return 0;
317: }

319: static PetscErrorCode PhysicsCreate_Advect(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
320: {
321:   Physics_Advect *advect;

324:   phys->field_desc = PhysicsFields_Advect;
325:   phys->riemann    = (RiemannFunction)PhysicsRiemann_Advect;
326:   PetscNew(&advect);
327:   phys->data = advect;
328:   PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
329:   {
330:     PetscInt two = 2, dof = 1;
331:     advect->soltype = ADVECT_SOL_TILTED;
332:     PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL);
333:     switch (advect->soltype) {
334:     case ADVECT_SOL_TILTED: {
335:       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
336:       two                           = 2;
337:       tilted->wind[0]               = 0.0;
338:       tilted->wind[1]               = 1.0;
339:       PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL);
340:       advect->inflowState = -2.0;
341:       PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL);
342:       phys->maxspeed = Norm2(tilted->wind);
343:     } break;
344:     case ADVECT_SOL_BUMP: {
345:       Physics_Advect_Bump *bump = &advect->sol.bump;
346:       two                       = 2;
347:       bump->center[0]           = 2.;
348:       bump->center[1]           = 0.;
349:       PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL);
350:       bump->radius = 0.9;
351:       PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL);
352:       bump->type = ADVECT_SOL_BUMP_CONE;
353:       PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL);
354:       phys->maxspeed = 3.; /* radius of mesh, kludge */
355:     } break;
356:     }
357:   }
358:   PetscOptionsHeadEnd();
359:   {
360:     const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
361:     DMLabel        label;

363:     DMGetLabel(dm, "Face Sets", &label);
364:     /* Register "canned" boundary conditions and defaults for where to apply. */
365:     PetscDSAddBoundary(prob, PETSC_TRUE, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Inflow, NULL, phys, NULL);
366:     PetscDSAddBoundary(prob, PETSC_TRUE, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Outflow, NULL, phys, NULL);
367:     /* Initial/transient solution with default boundary conditions */
368:     ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys);
369:     /* Register "canned" functionals */
370:     ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys);
371:   }
372:   return 0;
373: }

375: /******************* Shallow Water ********************/
376: typedef struct {
377:   PetscReal gravity;
378:   PetscReal boundaryHeight;
379:   struct {
380:     PetscInt Height;
381:     PetscInt Speed;
382:     PetscInt Energy;
383:   } functional;
384: } Physics_SW;
385: typedef struct {
386:   PetscScalar vals[0];
387:   PetscScalar h;
388:   PetscScalar uh[DIM];
389: } SWNode;

391: static const struct FieldDescription PhysicsFields_SW[] = {
392:   {"Height",   1  },
393:   {"Momentum", DIM},
394:   {NULL,       0  }
395: };

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

409:   Scale2(1. / x->h, x->uh, u);
410:   uhn  = Dot2(x->uh, n);
411:   f->h = uhn;
412:   for (i = 0; i < DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
413:   return 0;
414: }

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

425: static PetscErrorCode PhysicsRiemann_SW(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
426: {
427:   Physics_SW   *sw = (Physics_SW *)phys->data;
428:   PetscReal     cL, cR, speed, nn[DIM];
429:   const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
430:   SWNode        fL, fR;
431:   PetscInt      i;

435:   nn[0] = n[0];
436:   nn[1] = n[1];
437:   Normalize2(nn);
438:   SWFlux(phys, nn, uL, &fL);
439:   SWFlux(phys, nn, uR, &fR);
440:   cL    = PetscSqrtReal(sw->gravity * PetscRealPart(uL->h));
441:   cR    = PetscSqrtReal(sw->gravity * PetscRealPart(uR->h)); /* gravity wave speed */
442:   speed = PetscMax(PetscAbsScalar(Dot2(uL->uh, nn) / uL->h) + cL, PetscAbsScalar(Dot2(uR->uh, nn) / uR->h) + cR);
443:   for (i = 0; i < 1 + DIM; i++) flux[i] = (0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i])) * Norm2(n);
444:   return 0;
445: }

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

453:   dx[0] = x[0] - 1.5;
454:   dx[1] = x[1] - 1.0;
455:   r     = Norm2(dx);
456:   sigma = 0.5;
457:   u[0]  = 1 + 2 * PetscExpScalar(-PetscSqr(r) / (2 * PetscSqr(sigma)));
458:   u[1]  = 0.0;
459:   u[2]  = 0.0;
460:   return 0;
461: }

463: static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
464: {
465:   Physics       phys = (Physics)ctx;
466:   Physics_SW   *sw   = (Physics_SW *)phys->data;
467:   const SWNode *x    = (const SWNode *)xx;
468:   PetscScalar   u[2];
469:   PetscReal     h;

472:   h = PetscRealPart(x->h);
473:   Scale2(1. / x->h, x->uh, u);
474:   f[sw->functional.Height] = h;
475:   f[sw->functional.Speed]  = Norm2(u) + PetscSqrtReal(sw->gravity * h);
476:   f[sw->functional.Energy] = 0.5 * (Dot2(x->uh, u) + sw->gravity * PetscSqr(h));
477:   return 0;
478: }

480: static PetscErrorCode PhysicsCreate_SW(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
481: {
482:   Physics_SW *sw;

485:   phys->field_desc = PhysicsFields_SW;
486:   phys->riemann    = (RiemannFunction)PhysicsRiemann_SW;
487:   PetscNew(&sw);
488:   phys->data = sw;
489:   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
490:   {
491:     sw->gravity = 1.0;
492:     PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL);
493:   }
494:   PetscOptionsHeadEnd();
495:   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */

497:   {
498:     const PetscInt wallids[] = {100, 101, 200, 300};
499:     DMLabel        label;

501:     DMGetLabel(dm, "Face Sets", &label);
502:     PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_SW_Wall, NULL, phys, NULL);
503:     ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys);
504:     ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys);
505:     ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys);
506:     ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys);
507:   }
508:   return 0;
509: }

511: /******************* Euler ********************/
512: typedef struct {
513:   PetscScalar vals[0];
514:   PetscScalar r;
515:   PetscScalar ru[DIM];
516:   PetscScalar e;
517: } EulerNode;
518: typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode *, PetscScalar *);
519: typedef struct {
520:   PetscInt        npars;
521:   PetscReal       pars[DIM];
522:   EquationOfState pressure;
523:   EquationOfState sound;
524:   struct {
525:     PetscInt Density;
526:     PetscInt Momentum;
527:     PetscInt Energy;
528:     PetscInt Pressure;
529:     PetscInt Speed;
530:   } monitor;
531: } Physics_Euler;

533: static const struct FieldDescription PhysicsFields_Euler[] = {
534:   {"Density",  1  },
535:   {"Momentum", DIM},
536:   {"Energy",   1  },
537:   {NULL,       0  }
538: };

540: static PetscErrorCode Pressure_PG(const PetscReal *pars, const EulerNode *x, PetscScalar *p)
541: {
542:   PetscScalar ru2;

545:   ru2 = DotDIM(x->ru, x->ru);
546:   ru2 /= x->r;
547:   /* kinematic dof = params[0] */
548:   (*p) = 2.0 * (x->e - 0.5 * ru2) / pars[0];
549:   return 0;
550: }

552: static PetscErrorCode SpeedOfSound_PG(const PetscReal *pars, const EulerNode *x, PetscScalar *c)
553: {
554:   PetscScalar p;

557:   /* TODO remove direct usage of Pressure_PG */
558:   Pressure_PG(pars, x, &p);
559:   /* TODO check the sign of p */
560:   /* pars[1] = heat capacity ratio */
561:   (*c) = PetscSqrtScalar(pars[1] * p / x->r);
562:   return 0;
563: }

565: /*
566:  * x = (rho,rho*(u_1),...,rho*e)^T
567:  * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
568:  *
569:  * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
570:  *
571:  * */
572: static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
573: {
574:   Physics_Euler *eu = (Physics_Euler *)phys->data;
575:   PetscScalar    u, nu, p;
576:   PetscInt       i;

579:   u = DotDIM(x->ru, x->ru);
580:   u /= (x->r * x->r);
581:   nu = DotDIM(x->ru, n);
582:   /* TODO check the sign of p */
583:   eu->pressure(eu->pars, x, &p);
584:   f->r = nu * x->r;
585:   for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p;
586:   f->e = nu * (x->e + p);
587:   return 0;
588: }

590: /* PetscReal* => EulerNode* conversion */
591: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
592: {
593:   PetscInt    i;
594:   PetscScalar xn[DIM], xt[DIM];

597:   xG[0] = xI[0];
598:   NormalSplitDIM(n, xI + 1, xn, xt);
599:   for (i = 0; i < DIM; i++) xG[i + 1] = -xn[i] + xt[i];
600:   xG[DIM + 1] = xI[DIM + 1];
601:   return 0;
602: }

604: /* PetscReal* => EulerNode* conversion */
605: static PetscErrorCode PhysicsRiemann_Euler_Rusanov(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
606: {
607:   Physics_Euler   *eu = (Physics_Euler *)phys->data;
608:   PetscScalar      cL, cR, speed;
609:   const EulerNode *uL = (const EulerNode *)xL, *uR = (const EulerNode *)xR;
610:   EulerNode        fL, fR;
611:   PetscInt         i;

615:   EulerFlux(phys, n, uL, &fL);
616:   EulerFlux(phys, n, uR, &fR);
617:   eu->sound(eu->pars, uL, &cL);
618:   eu->sound(eu->pars, uR, &cR);
619:   speed = PetscMax(cL, cR) + PetscMax(PetscAbsScalar(DotDIM(uL->ru, n) / NormDIM(n)), PetscAbsScalar(DotDIM(uR->ru, n) / NormDIM(n)));
620:   for (i = 0; i < 2 + DIM; i++) flux[i] = 0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i]);
621:   return 0;
622: }

624: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
625: {
626:   PetscInt i;

630:   u[0]       = 1.0;
631:   u[DIM + 1] = 1.0 + PetscAbsReal(x[0]);
632:   for (i = 1; i < DIM + 1; i++) u[i] = 0.0;
633:   return 0;
634: }

636: static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
637: {
638:   Physics          phys = (Physics)ctx;
639:   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
640:   const EulerNode *x    = (const EulerNode *)xx;
641:   PetscScalar      p;

644:   f[eu->monitor.Density]  = x->r;
645:   f[eu->monitor.Momentum] = NormDIM(x->ru);
646:   f[eu->monitor.Energy]   = x->e;
647:   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
648:   eu->pressure(eu->pars, x, &p);
649:   f[eu->monitor.Pressure] = p;
650:   return 0;
651: }

653: static PetscErrorCode PhysicsCreate_Euler(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
654: {
655:   Physics_Euler *eu;

658:   phys->field_desc = PhysicsFields_Euler;
659:   phys->riemann    = (RiemannFunction)PhysicsRiemann_Euler_Rusanov;
660:   PetscNew(&eu);
661:   phys->data = eu;
662:   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
663:   {
664:     eu->pars[0] = 3.0;
665:     eu->pars[1] = 1.67;
666:     PetscOptionsReal("-eu_f", "Degrees of freedom", "", eu->pars[0], &eu->pars[0], NULL);
667:     PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->pars[1], &eu->pars[1], NULL);
668:   }
669:   PetscOptionsHeadEnd();
670:   eu->pressure   = Pressure_PG;
671:   eu->sound      = SpeedOfSound_PG;
672:   phys->maxspeed = 1.0;
673:   {
674:     const PetscInt wallids[] = {100, 101, 200, 300};
675:     DMLabel        label;

677:     DMGetLabel(dm, "Face Sets", &label);
678:     PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_Euler_Wall, NULL, phys, NULL);
679:     ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys);
680:     ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys);
681:     ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys);
682:     ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys);
683:     ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys);
684:     ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys);
685:   }
686:   return 0;
687: }

689: PetscErrorCode ConstructCellBoundary(DM dm, User user)
690: {
691:   const char     *name   = "Cell Sets";
692:   const char     *bdname = "split faces";
693:   IS              regionIS, innerIS;
694:   const PetscInt *regions, *cells;
695:   PetscInt        numRegions, innerRegion, numCells, c;
696:   PetscInt        cStart, cEnd, cEndInterior, fStart, fEnd;
697:   PetscBool       hasLabel;

700:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
701:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
702:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);

704:   DMHasLabel(dm, name, &hasLabel);
705:   if (!hasLabel) return 0;
706:   DMGetLabelSize(dm, name, &numRegions);
707:   if (numRegions != 2) return 0;
708:   /* Get the inner id */
709:   DMGetLabelIdIS(dm, name, &regionIS);
710:   ISGetIndices(regionIS, &regions);
711:   innerRegion = regions[0];
712:   ISRestoreIndices(regionIS, &regions);
713:   ISDestroy(&regionIS);
714:   /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
715:   DMGetStratumIS(dm, name, innerRegion, &innerIS);
716:   ISGetLocalSize(innerIS, &numCells);
717:   ISGetIndices(innerIS, &cells);
718:   DMCreateLabel(dm, bdname);
719:   for (c = 0; c < numCells; ++c) {
720:     const PetscInt  cell = cells[c];
721:     const PetscInt *faces;
722:     PetscInt        numFaces, f;

725:     DMPlexGetConeSize(dm, cell, &numFaces);
726:     DMPlexGetCone(dm, cell, &faces);
727:     for (f = 0; f < numFaces; ++f) {
728:       const PetscInt  face = faces[f];
729:       const PetscInt *neighbors;
730:       PetscInt        nC, regionA, regionB;

733:       DMPlexGetSupportSize(dm, face, &nC);
734:       if (nC != 2) continue;
735:       DMPlexGetSupport(dm, face, &neighbors);
736:       if ((neighbors[0] >= cEndInterior) || (neighbors[1] >= cEndInterior)) continue;
739:       DMGetLabelValue(dm, name, neighbors[0], &regionA);
740:       DMGetLabelValue(dm, name, neighbors[1], &regionB);
743:       if (regionA != regionB) DMSetLabelValue(dm, bdname, faces[f], 1);
744:     }
745:   }
746:   ISRestoreIndices(innerIS, &cells);
747:   ISDestroy(&innerIS);
748:   {
749:     DMLabel label;

751:     DMGetLabel(dm, bdname, &label);
752:     DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD);
753:   }
754:   return 0;
755: }

757: /* Right now, I have just added duplicate faces, which see both cells. We can
758: - Add duplicate vertices and decouple the face cones
759: - Disconnect faces from cells across the rotation gap
760: */
761: PetscErrorCode SplitFaces(DM *dmSplit, const char labelName[], User user)
762: {
763:   DM              dm = *dmSplit, sdm;
764:   PetscSF         sfPoint, gsfPoint;
765:   PetscSection    coordSection, newCoordSection;
766:   Vec             coordinates;
767:   IS              idIS;
768:   const PetscInt *ids;
769:   PetscInt       *newpoints;
770:   PetscInt        dim, depth, maxConeSize, maxSupportSize, numLabels, numGhostCells;
771:   PetscInt        numFS, fs, pStart, pEnd, p, cEnd, cEndInterior, vStart, vEnd, v, fStart, fEnd, newf, d, l;
772:   PetscBool       hasLabel;

775:   DMHasLabel(dm, labelName, &hasLabel);
776:   if (!hasLabel) return 0;
777:   DMCreate(PetscObjectComm((PetscObject)dm), &sdm);
778:   DMSetType(sdm, DMPLEX);
779:   DMGetDimension(dm, &dim);
780:   DMSetDimension(sdm, dim);

782:   DMGetLabelIdIS(dm, labelName, &idIS);
783:   ISGetLocalSize(idIS, &numFS);
784:   ISGetIndices(idIS, &ids);

786:   user->numSplitFaces = 0;
787:   for (fs = 0; fs < numFS; ++fs) {
788:     PetscInt numBdFaces;

790:     DMGetStratumSize(dm, labelName, ids[fs], &numBdFaces);
791:     user->numSplitFaces += numBdFaces;
792:   }
793:   DMPlexGetChart(dm, &pStart, &pEnd);
794:   pEnd += user->numSplitFaces;
795:   DMPlexSetChart(sdm, pStart, pEnd);
796:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
797:   DMPlexSetGhostCellStratum(sdm, cEndInterior, PETSC_DETERMINE);
798:   DMPlexGetHeightStratum(dm, 0, NULL, &cEnd);
799:   numGhostCells = cEnd - cEndInterior;
800:   /* Set cone and support sizes */
801:   DMPlexGetDepth(dm, &depth);
802:   for (d = 0; d <= depth; ++d) {
803:     DMPlexGetDepthStratum(dm, d, &pStart, &pEnd);
804:     for (p = pStart; p < pEnd; ++p) {
805:       PetscInt newp = p;
806:       PetscInt size;

808:       DMPlexGetConeSize(dm, p, &size);
809:       DMPlexSetConeSize(sdm, newp, size);
810:       DMPlexGetSupportSize(dm, p, &size);
811:       DMPlexSetSupportSize(sdm, newp, size);
812:     }
813:   }
814:   DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd);
815:   for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
816:     IS              faceIS;
817:     const PetscInt *faces;
818:     PetscInt        numFaces, f;

820:     DMGetStratumIS(dm, labelName, ids[fs], &faceIS);
821:     ISGetLocalSize(faceIS, &numFaces);
822:     ISGetIndices(faceIS, &faces);
823:     for (f = 0; f < numFaces; ++f, ++newf) {
824:       PetscInt size;

826:       /* Right now I think that both faces should see both cells */
827:       DMPlexGetConeSize(dm, faces[f], &size);
828:       DMPlexSetConeSize(sdm, newf, size);
829:       DMPlexGetSupportSize(dm, faces[f], &size);
830:       DMPlexSetSupportSize(sdm, newf, size);
831:     }
832:     ISRestoreIndices(faceIS, &faces);
833:     ISDestroy(&faceIS);
834:   }
835:   DMSetUp(sdm);
836:   /* Set cones and supports */
837:   DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize);
838:   PetscMalloc1(PetscMax(maxConeSize, maxSupportSize), &newpoints);
839:   DMPlexGetChart(dm, &pStart, &pEnd);
840:   for (p = pStart; p < pEnd; ++p) {
841:     const PetscInt *points, *orientations;
842:     PetscInt        size, i, newp = p;

844:     DMPlexGetConeSize(dm, p, &size);
845:     DMPlexGetCone(dm, p, &points);
846:     DMPlexGetConeOrientation(dm, p, &orientations);
847:     for (i = 0; i < size; ++i) newpoints[i] = points[i];
848:     DMPlexSetCone(sdm, newp, newpoints);
849:     DMPlexSetConeOrientation(sdm, newp, orientations);
850:     DMPlexGetSupportSize(dm, p, &size);
851:     DMPlexGetSupport(dm, p, &points);
852:     for (i = 0; i < size; ++i) newpoints[i] = points[i];
853:     DMPlexSetSupport(sdm, newp, newpoints);
854:   }
855:   PetscFree(newpoints);
856:   for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
857:     IS              faceIS;
858:     const PetscInt *faces;
859:     PetscInt        numFaces, f;

861:     DMGetStratumIS(dm, labelName, ids[fs], &faceIS);
862:     ISGetLocalSize(faceIS, &numFaces);
863:     ISGetIndices(faceIS, &faces);
864:     for (f = 0; f < numFaces; ++f, ++newf) {
865:       const PetscInt *points;

867:       DMPlexGetCone(dm, faces[f], &points);
868:       DMPlexSetCone(sdm, newf, points);
869:       DMPlexGetSupport(dm, faces[f], &points);
870:       DMPlexSetSupport(sdm, newf, points);
871:     }
872:     ISRestoreIndices(faceIS, &faces);
873:     ISDestroy(&faceIS);
874:   }
875:   ISRestoreIndices(idIS, &ids);
876:   ISDestroy(&idIS);
877:   DMPlexStratify(sdm);
878:   /* Convert coordinates */
879:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
880:   DMGetCoordinateSection(dm, &coordSection);
881:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection);
882:   PetscSectionSetNumFields(newCoordSection, 1);
883:   PetscSectionSetFieldComponents(newCoordSection, 0, dim);
884:   PetscSectionSetChart(newCoordSection, vStart, vEnd);
885:   for (v = vStart; v < vEnd; ++v) {
886:     PetscSectionSetDof(newCoordSection, v, dim);
887:     PetscSectionSetFieldDof(newCoordSection, v, 0, dim);
888:   }
889:   PetscSectionSetUp(newCoordSection);
890:   DMSetCoordinateSection(sdm, PETSC_DETERMINE, newCoordSection);
891:   PetscSectionDestroy(&newCoordSection); /* relinquish our reference */
892:   DMGetCoordinatesLocal(dm, &coordinates);
893:   DMSetCoordinatesLocal(sdm, coordinates);
894:   /* Convert labels */
895:   DMGetNumLabels(dm, &numLabels);
896:   for (l = 0; l < numLabels; ++l) {
897:     const char *lname;
898:     PetscBool   isDepth;

900:     DMGetLabelName(dm, l, &lname);
901:     PetscStrcmp(lname, "depth", &isDepth);
902:     if (isDepth) continue;
903:     DMCreateLabel(sdm, lname);
904:     DMGetLabelIdIS(dm, lname, &idIS);
905:     ISGetLocalSize(idIS, &numFS);
906:     ISGetIndices(idIS, &ids);
907:     for (fs = 0; fs < numFS; ++fs) {
908:       IS              pointIS;
909:       const PetscInt *points;
910:       PetscInt        numPoints;

912:       DMGetStratumIS(dm, lname, ids[fs], &pointIS);
913:       ISGetLocalSize(pointIS, &numPoints);
914:       ISGetIndices(pointIS, &points);
915:       for (p = 0; p < numPoints; ++p) {
916:         PetscInt newpoint = points[p];

918:         DMSetLabelValue(sdm, lname, newpoint, ids[fs]);
919:       }
920:       ISRestoreIndices(pointIS, &points);
921:       ISDestroy(&pointIS);
922:     }
923:     ISRestoreIndices(idIS, &ids);
924:     ISDestroy(&idIS);
925:   }
926:   /* Convert pointSF */
927:   const PetscSFNode *remotePoints;
928:   PetscSFNode       *gremotePoints;
929:   const PetscInt    *localPoints;
930:   PetscInt          *glocalPoints, *newLocation, *newRemoteLocation;
931:   PetscInt           numRoots, numLeaves;
932:   PetscMPIInt        size;

934:   MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size);
935:   DMGetPointSF(dm, &sfPoint);
936:   DMGetPointSF(sdm, &gsfPoint);
937:   DMPlexGetChart(dm, &pStart, &pEnd);
938:   PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints);
939:   if (numRoots >= 0) {
940:     PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation);
941:     for (l = 0; l < numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? numGhostCells : 0); */
942:     PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
943:     PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE);
944:     PetscMalloc1(numLeaves, &glocalPoints);
945:     PetscMalloc1(numLeaves, &gremotePoints);
946:     for (l = 0; l < numLeaves; ++l) {
947:       glocalPoints[l]        = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + numGhostCells : localPoints[l]; */
948:       gremotePoints[l].rank  = remotePoints[l].rank;
949:       gremotePoints[l].index = newRemoteLocation[localPoints[l]];
950:     }
951:     PetscFree2(newLocation, newRemoteLocation);
952:     PetscSFSetGraph(gsfPoint, numRoots + numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER);
953:   }
954:   DMDestroy(dmSplit);
955:   *dmSplit = sdm;
956:   return 0;
957: }

959: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
960: {
961:   PetscSF      sfPoint;
962:   PetscSection coordSection;
963:   Vec          coordinates;
964:   PetscSection sectionCell;
965:   PetscScalar *part;
966:   PetscInt     cStart, cEnd, c;
967:   PetscMPIInt  rank;

970:   DMGetCoordinateSection(dm, &coordSection);
971:   DMGetCoordinatesLocal(dm, &coordinates);
972:   DMClone(dm, dmCell);
973:   DMGetPointSF(dm, &sfPoint);
974:   DMSetPointSF(*dmCell, sfPoint);
975:   DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection);
976:   DMSetCoordinatesLocal(*dmCell, coordinates);
977:   MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank);
978:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell);
979:   DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd);
980:   PetscSectionSetChart(sectionCell, cStart, cEnd);
981:   for (c = cStart; c < cEnd; ++c) PetscSectionSetDof(sectionCell, c, 1);
982:   PetscSectionSetUp(sectionCell);
983:   DMSetLocalSection(*dmCell, sectionCell);
984:   PetscSectionDestroy(&sectionCell);
985:   DMCreateLocalVector(*dmCell, partition);
986:   PetscObjectSetName((PetscObject)*partition, "partition");
987:   VecGetArray(*partition, &part);
988:   for (c = cStart; c < cEnd; ++c) {
989:     PetscScalar *p;

991:     DMPlexPointLocalRef(*dmCell, c, part, &p);
992:     p[0] = rank;
993:   }
994:   VecRestoreArray(*partition, &part);
995:   return 0;
996: }

998: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
999: {
1000:   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
1001:   PetscSection       coordSection;
1002:   Vec                coordinates, facegeom, cellgeom;
1003:   PetscSection       sectionMass;
1004:   PetscScalar       *m;
1005:   const PetscScalar *fgeom, *cgeom, *coords;
1006:   PetscInt           vStart, vEnd, v;

1009:   DMConvert(dm, DMPLEX, &plex);
1010:   DMGetCoordinateSection(dm, &coordSection);
1011:   DMGetCoordinatesLocal(dm, &coordinates);
1012:   DMClone(dm, &dmMass);
1013:   DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection);
1014:   DMSetCoordinatesLocal(dmMass, coordinates);
1015:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass);
1016:   DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd);
1017:   PetscSectionSetChart(sectionMass, vStart, vEnd);
1018:   for (v = vStart; v < vEnd; ++v) {
1019:     PetscInt numFaces;

1021:     DMPlexGetSupportSize(dmMass, v, &numFaces);
1022:     PetscSectionSetDof(sectionMass, v, numFaces * numFaces);
1023:   }
1024:   PetscSectionSetUp(sectionMass);
1025:   DMSetLocalSection(dmMass, sectionMass);
1026:   PetscSectionDestroy(&sectionMass);
1027:   DMGetLocalVector(dmMass, massMatrix);
1028:   VecGetArray(*massMatrix, &m);
1029:   DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL);
1030:   VecGetDM(facegeom, &dmFace);
1031:   VecGetArrayRead(facegeom, &fgeom);
1032:   VecGetDM(cellgeom, &dmCell);
1033:   VecGetArrayRead(cellgeom, &cgeom);
1034:   DMGetCoordinateDM(dm, &dmCoord);
1035:   VecGetArrayRead(coordinates, &coords);
1036:   for (v = vStart; v < vEnd; ++v) {
1037:     const PetscInt        *faces;
1038:     const PetscFVFaceGeom *fgA, *fgB, *cg;
1039:     const PetscScalar     *vertex;
1040:     PetscInt               numFaces, sides[2], f, g;

1042:     DMPlexPointLocalRead(dmCoord, v, coords, &vertex);
1043:     DMPlexGetSupportSize(dmMass, v, &numFaces);
1044:     DMPlexGetSupport(dmMass, v, &faces);
1045:     for (f = 0; f < numFaces; ++f) {
1046:       sides[0] = faces[f];
1047:       DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA);
1048:       for (g = 0; g < numFaces; ++g) {
1049:         const PetscInt *cells = NULL;
1050:         PetscReal       area  = 0.0;
1051:         PetscInt        numCells;

1053:         sides[1] = faces[g];
1054:         DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB);
1055:         DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells);
1057:         DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg);
1058:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
1059:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
1060:         m[f * numFaces + g] = Dot2(fgA->normal, fgB->normal) * area * 0.5;
1061:         DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells);
1062:       }
1063:     }
1064:   }
1065:   VecRestoreArrayRead(facegeom, &fgeom);
1066:   VecRestoreArrayRead(cellgeom, &cgeom);
1067:   VecRestoreArrayRead(coordinates, &coords);
1068:   VecRestoreArray(*massMatrix, &m);
1069:   DMDestroy(&dmMass);
1070:   DMDestroy(&plex);
1071:   return 0;
1072: }

1074: PetscErrorCode SetUpLocalSpace(DM dm, User user)
1075: {
1076:   PetscSection stateSection;
1077:   Physics      phys;
1078:   PetscInt     dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, cEndInterior, c, i;

1081:   DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd);
1082:   DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL);
1083:   PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection);
1084:   phys = user->model->physics;
1085:   PetscSectionSetNumFields(stateSection, phys->nfields);
1086:   for (i = 0; i < phys->nfields; i++) {
1087:     PetscSectionSetFieldName(stateSection, i, phys->field_desc[i].name);
1088:     PetscSectionSetFieldComponents(stateSection, i, phys->field_desc[i].dof);
1089:   }
1090:   PetscSectionSetChart(stateSection, cStart, cEnd);
1091:   for (c = cStart; c < cEnd; ++c) {
1092:     for (i = 0; i < phys->nfields; i++) PetscSectionSetFieldDof(stateSection, c, i, phys->field_desc[i].dof);
1093:     PetscSectionSetDof(stateSection, c, dof);
1094:   }
1095:   for (c = cEndInterior; c < cEnd; ++c) PetscSectionSetConstraintDof(stateSection, c, dof);
1096:   PetscSectionSetUp(stateSection);
1097:   PetscMalloc1(dof, &cind);
1098:   for (d = 0; d < dof; ++d) cind[d] = d;
1099: #if 0
1100:   for (c = cStart; c < cEnd; ++c) {
1101:     PetscInt val;

1103:     DMGetLabelValue(dm, "vtk", c, &val);
1104:     if (val < 0) PetscSectionSetConstraintIndices(stateSection, c, cind);
1105:   }
1106: #endif
1107:   for (c = cEndInterior; c < cEnd; ++c) PetscSectionSetConstraintIndices(stateSection, c, cind);
1108:   PetscFree(cind);
1109:   PetscSectionGetStorageSize(stateSection, &stateSize);
1110:   DMSetLocalSection(dm, stateSection);
1111:   PetscSectionDestroy(&stateSection);
1112:   return 0;
1113: }

1115: #if 0
1116: PetscErrorCode SetUpBoundaries(DM dm, User user)
1117: {
1118:   Model          mod = user->model;
1119:   BoundaryLink   b;

1122:   PetscOptionsBegin(PetscObjectComm((PetscObject)dm),NULL,"Boundary condition options","");
1123:   for (b = mod->boundary; b; b=b->next) {
1124:     char      optname[512];
1125:     PetscInt  ids[512],len = 512;
1126:     PetscBool flg;
1127:     PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name);
1128:     PetscMemzero(ids,sizeof(ids));
1129:     PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg);
1130:     if (flg) {
1131:       /* TODO: check all IDs to make sure they exist in the mesh */
1132:       PetscFree(b->ids);
1133:       b->numids = len;
1134:       PetscMalloc1(len,&b->ids);
1135:       PetscArraycpy(b->ids,ids,len);
1136:     }
1137:   }
1138:   PetscOptionsEnd();
1139:   return 0;
1140: }
1141: #endif

1143: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
1144: static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
1145: {
1147:   mod->solution    = func;
1148:   mod->solutionctx = ctx;
1149:   return 0;
1150: }

1152: static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
1153: {
1154:   FunctionalLink link, *ptr;
1155:   PetscInt       lastoffset = -1;

1158:   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1159:   PetscNew(&link);
1160:   PetscStrallocpy(name, &link->name);
1161:   link->offset = lastoffset + 1;
1162:   link->func   = func;
1163:   link->ctx    = ctx;
1164:   link->next   = NULL;
1165:   *ptr         = link;
1166:   *offset      = link->offset;
1167:   return 0;
1168: }

1170: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptions *PetscOptionsObject)
1171: {
1172:   PetscInt       i, j;
1173:   FunctionalLink link;
1174:   char          *names[256];

1177:   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
1178:   PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL);
1179:   /* Create list of functionals that will be computed somehow */
1180:   PetscMalloc1(mod->numMonitored, &mod->functionalMonitored);
1181:   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1182:   PetscMalloc1(mod->numMonitored, &mod->functionalCall);
1183:   mod->numCall = 0;
1184:   for (i = 0; i < mod->numMonitored; i++) {
1185:     for (link = mod->functionalRegistry; link; link = link->next) {
1186:       PetscBool match;
1187:       PetscStrcasecmp(names[i], link->name, &match);
1188:       if (match) break;
1189:     }
1191:     mod->functionalMonitored[i] = link;
1192:     for (j = 0; j < i; j++) {
1193:       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
1194:     }
1195:     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1196:   next_name:
1197:     PetscFree(names[i]);
1198:   }

1200:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1201:   mod->maxComputed = -1;
1202:   for (link = mod->functionalRegistry; link; link = link->next) {
1203:     for (i = 0; i < mod->numCall; i++) {
1204:       FunctionalLink call = mod->functionalCall[i];
1205:       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
1206:     }
1207:   }
1208:   return 0;
1209: }

1211: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
1212: {
1213:   FunctionalLink l, next;

1216:   if (!link) return 0;
1217:   l     = *link;
1218:   *link = NULL;
1219:   for (; l; l = next) {
1220:     next = l->next;
1221:     PetscFree(l->name);
1222:     PetscFree(l);
1223:   }
1224:   return 0;
1225: }

1227: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
1228: {
1229:   DM                 plex, dmCell;
1230:   Model              mod = user->model;
1231:   Vec                cellgeom;
1232:   const PetscScalar *cgeom;
1233:   PetscScalar       *x;
1234:   PetscInt           cStart, cEnd, c;

1237:   DMConvert(dm, DMPLEX, &plex);
1238:   DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL);
1239:   VecGetDM(cellgeom, &dmCell);
1240:   DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1241:   VecGetArrayRead(cellgeom, &cgeom);
1242:   VecGetArray(X, &x);
1243:   for (c = cStart; c < cEnd; ++c) {
1244:     const PetscFVCellGeom *cg;
1245:     PetscScalar           *xc;

1247:     DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1248:     DMPlexPointGlobalRef(dm, c, x, &xc);
1249:     if (xc) (*mod->solution)(mod, 0.0, cg->centroid, xc, mod->solutionctx);
1250:   }
1251:   VecRestoreArrayRead(cellgeom, &cgeom);
1252:   VecRestoreArray(X, &x);
1253:   DMDestroy(&plex);
1254:   return 0;
1255: }

1257: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1258: {
1260:   PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);
1261:   PetscViewerSetType(*viewer, PETSCVIEWERVTK);
1262:   PetscViewerFileSetName(*viewer, filename);
1263:   return 0;
1264: }

1266: static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
1267: {
1268:   User        user = (User)ctx;
1269:   DM          dm, plex;
1270:   PetscViewer viewer;
1271:   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
1272:   PetscReal   xnorm;

1275:   PetscObjectSetName((PetscObject)X, "solution");
1276:   VecGetDM(X, &dm);
1277:   VecNorm(X, NORM_INFINITY, &xnorm);
1278:   if (stepnum >= 0) { /* No summary for final time */
1279:     Model              mod = user->model;
1280:     Vec                cellgeom;
1281:     PetscInt           c, cStart, cEnd, fcount, i;
1282:     size_t             ftableused, ftablealloc;
1283:     const PetscScalar *cgeom, *x;
1284:     DM                 dmCell;
1285:     PetscReal         *fmin, *fmax, *fintegral, *ftmp;

1287:     DMConvert(dm, DMPLEX, &plex);
1288:     DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL);
1289:     fcount = mod->maxComputed + 1;
1290:     PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp);
1291:     for (i = 0; i < fcount; i++) {
1292:       fmin[i]      = PETSC_MAX_REAL;
1293:       fmax[i]      = PETSC_MIN_REAL;
1294:       fintegral[i] = 0;
1295:     }
1296:     DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd);
1297:     VecGetDM(cellgeom, &dmCell);
1298:     VecGetArrayRead(cellgeom, &cgeom);
1299:     VecGetArrayRead(X, &x);
1300:     for (c = cStart; c < cEnd; ++c) {
1301:       const PetscFVCellGeom *cg;
1302:       const PetscScalar     *cx;
1303:       DMPlexPointLocalRead(dmCell, c, cgeom, &cg);
1304:       DMPlexPointGlobalRead(dm, c, x, &cx);
1305:       if (!cx) continue; /* not a global cell */
1306:       for (i = 0; i < mod->numCall; i++) {
1307:         FunctionalLink flink = mod->functionalCall[i];
1308:         (*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx);
1309:       }
1310:       for (i = 0; i < fcount; i++) {
1311:         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1312:         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1313:         fintegral[i] += cg->volume * ftmp[i];
1314:       }
1315:     }
1316:     VecRestoreArrayRead(cellgeom, &cgeom);
1317:     VecRestoreArrayRead(X, &x);
1318:     DMDestroy(&plex);
1319:     MPI_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts));
1320:     MPI_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts));
1321:     MPI_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts));

1323:     ftablealloc = fcount * 100;
1324:     ftableused  = 0;
1325:     PetscMalloc1(ftablealloc, &ftable);
1326:     for (i = 0; i < mod->numMonitored; i++) {
1327:       size_t         countused;
1328:       char           buffer[256], *p;
1329:       FunctionalLink flink = mod->functionalMonitored[i];
1330:       PetscInt       id    = flink->offset;
1331:       if (i % 3) {
1332:         PetscArraycpy(buffer, "  ", 2);
1333:         p = buffer + 2;
1334:       } else if (i) {
1335:         char newline[] = "\n";
1336:         PetscArraycpy(buffer, newline, sizeof(newline) - 1);
1337:         p = buffer + sizeof(newline) - 1;
1338:       } else {
1339:         p = buffer;
1340:       }
1341:       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]);
1342:       countused += p - buffer;
1343:       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1344:         char *ftablenew;
1345:         ftablealloc = 2 * ftablealloc + countused;
1346:         PetscMalloc(ftablealloc, &ftablenew);
1347:         PetscArraycpy(ftablenew, ftable, ftableused);
1348:         PetscFree(ftable);
1349:         ftable = ftablenew;
1350:       }
1351:       PetscArraycpy(ftable + ftableused, buffer, countused);
1352:       ftableused += countused;
1353:       ftable[ftableused] = 0;
1354:     }
1355:     PetscFree4(fmin, fmax, fintegral, ftmp);

1357:     PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : "");
1358:     PetscFree(ftable);
1359:   }
1360:   if (user->vtkInterval < 1) return 0;
1361:   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1362:     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
1363:       TSGetStepNumber(ts, &stepnum);
1364:     }
1365:     PetscSNPrintf(filename, sizeof filename, "ex11-%03" PetscInt_FMT ".vtu", stepnum);
1366:     OutputVTK(dm, filename, &viewer);
1367:     VecView(X, viewer);
1368:     PetscViewerDestroy(&viewer);
1369:   }
1370:   return 0;
1371: }

1373: static PetscErrorCode OutputBIN(DM dm, const char *filename, PetscViewer *viewer)
1374: {
1376:   PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer);
1377:   PetscViewerSetType(*viewer, PETSCVIEWERBINARY);
1378:   PetscViewerFileSetMode(*viewer, FILE_MODE_WRITE);
1379:   PetscViewerFileSetName(*viewer, filename);
1380:   return 0;
1381: }

1383: static PetscErrorCode TestMonitor(DM dm, const char *filename, Vec X, PetscReal time)
1384: {
1385:   Vec         odesolution;
1386:   PetscInt    Nr;
1387:   PetscBool   equal;
1388:   PetscReal   timeread;
1389:   PetscViewer viewer;

1392:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer);
1393:   VecCreate(PETSC_COMM_WORLD, &odesolution);
1394:   VecLoad(odesolution, viewer);
1395:   VecEqual(X, odesolution, &equal);
1397:   else
1398:   {
1399:     PetscPrintf(PETSC_COMM_WORLD, "IO test OK for Vec\n");
1400:   }
1401:   /*Nr   = 1;
1402:    PetscRealLoad(Nr,&Nr,&timeread,viewer);*/
1403:   PetscViewerBinaryRead(viewer, &timeread, 1, NULL, PETSC_REAL);

1406:   else
1407:   {
1408:     PetscPrintf(PETSC_COMM_WORLD, "IO test OK for PetscReal\n");
1409:   }

1411:   PetscViewerDestroy(&viewer);
1412:   VecDestroy(&odesolution);
1413:   return 0;
1414: }

1416: static PetscErrorCode MonitorBIN(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
1417: {
1418:   User        user = (User)ctx;
1419:   DM          dm;
1420:   PetscViewer viewer;
1421:   char        filename[PETSC_MAX_PATH_LEN];

1424:   VecGetDM(X, &dm);
1425:   PetscSNPrintf(filename, sizeof filename, "ex11-SA-%06d.bin", stepnum);
1426:   OutputBIN(dm, filename, &viewer);
1427:   VecView(X, viewer);
1428:   PetscRealView(1, &time, viewer);
1429:   /* PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL);*/
1430:   PetscViewerDestroy(&viewer);
1431:   TestMonitor(dm, filename, X, time);
1432:   return 0;
1433: }

1435: int main(int argc, char **argv)
1436: {
1437:   MPI_Comm          comm;
1438:   PetscFV           fvm;
1439:   User              user;
1440:   Model             mod;
1441:   Physics           phys;
1442:   DM                dm, plex;
1443:   PetscReal         ftime, cfl, dt, minRadius;
1444:   PetscInt          dim, nsteps;
1445:   TS                ts;
1446:   TSConvergedReason reason;
1447:   Vec               X;
1448:   PetscViewer       viewer;
1449:   PetscMPIInt       rank;
1450:   PetscBool         vtkCellGeom, splitFaces;
1451:   PetscInt          overlap, f;
1452:   char              filename[PETSC_MAX_PATH_LEN] = "sevenside.exo";

1455:   PetscInitialize(&argc, &argv, (char *)0, help);
1456:   comm = PETSC_COMM_WORLD;
1457:   MPI_Comm_rank(comm, &rank);

1459:   PetscNew(&user);
1460:   PetscNew(&user->model);
1461:   PetscNew(&user->model->physics);
1462:   mod       = user->model;
1463:   phys      = mod->physics;
1464:   mod->comm = comm;

1466:   /* Register physical models to be available on the command line */
1467:   PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect);
1468:   PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW);
1469:   PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler);

1471:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1472:   {
1473:     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1474:     PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL);
1475:     PetscOptionsString("-f", "Exodus.II filename to read", "", filename, filename, sizeof(filename), NULL);
1476:     splitFaces = PETSC_FALSE;
1477:     PetscOptionsBool("-ufv_split_faces", "Split faces between cell sets", "", splitFaces, &splitFaces, NULL);
1478:     overlap = 1;
1479:     PetscOptionsInt("-ufv_mesh_overlap", "Number of cells to overlap partitions", "", overlap, &overlap, NULL);
1480:     user->vtkInterval = 1;
1481:     PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL);
1482:     vtkCellGeom = PETSC_FALSE;
1483:     PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL);
1484:   }
1485:   PetscOptionsEnd();
1486:   DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, &dm);
1487:   DMViewFromOptions(dm, NULL, "-dm_view");
1488:   DMGetDimension(dm, &dim);

1490:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1491:   {
1492:     PetscDS prob;
1493:     PetscErrorCode (*physcreate)(PetscDS, Model, Physics);
1494:     char physname[256] = "advect";

1496:     DMCreateLabel(dm, "Face Sets");
1497:     PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL);
1498:     PetscFunctionListFind(PhysicsList, physname, &physcreate);
1499:     PetscMemzero(phys, sizeof(struct _n_Physics));
1500:     DMGetDS(dm, &prob);
1501:     (*physcreate)(prob, mod, phys);
1502:     mod->maxspeed = phys->maxspeed;
1503:     /* Count number of fields and dofs */
1504:     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;

1508:     ModelFunctionalSetFromOptions(mod, PetscOptionsObject);
1509:   }
1510:   PetscOptionsEnd();
1511:   {
1512:     DM dmDist;

1514:     DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE);
1515:     DMPlexDistribute(dm, overlap, NULL, &dmDist);
1516:     if (dmDist) {
1517:       DMDestroy(&dm);
1518:       dm = dmDist;
1519:     }
1520:   }
1521:   DMSetFromOptions(dm);
1522:   {
1523:     DM gdm;

1525:     DMPlexConstructGhostCells(dm, NULL, NULL, &gdm);
1526:     DMDestroy(&dm);
1527:     dm = gdm;
1528:     DMViewFromOptions(dm, NULL, "-dm_view");
1529:   }
1530:   if (splitFaces) ConstructCellBoundary(dm, user);
1531:   SplitFaces(&dm, "split faces", user);

1533:   PetscFVCreate(comm, &fvm);
1534:   PetscFVSetFromOptions(fvm);
1535:   DMSetNumFields(dm, phys->nfields);
1536:   for (f = 0; f < phys->nfields; ++f) DMSetField(dm, f, (PetscObject)fvm);
1537:   PetscFVSetNumComponents(fvm, phys->dof);
1538:   PetscFVSetSpatialDimension(fvm, dim);

1540:   /* Set up DM with section describing local vector and configure local vector. */
1541:   SetUpLocalSpace(dm, user);

1543:   TSCreate(comm, &ts);
1544:   TSSetType(ts, TSSSP);
1545:   /* TSSetType(ts, TSRK); */
1546:   TSSetDM(ts, dm);
1547:   /* TSMonitorSet(ts,MonitorVTK,user,NULL); */
1548:   TSMonitorSet(ts, MonitorBIN, user, NULL);
1549:   DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user);
1550:   DMCreateGlobalVector(dm, &X);
1551:   PetscObjectSetName((PetscObject)X, "solution");
1552:   SetInitialCondition(dm, X, user);
1553:   DMConvert(dm, DMPLEX, &plex);
1554:   if (vtkCellGeom) {
1555:     DM  dmCell;
1556:     Vec cellgeom, partition;

1558:     DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL);
1559:     OutputVTK(dm, "ex11-cellgeom.vtk", &viewer);
1560:     VecView(cellgeom, viewer);
1561:     PetscViewerDestroy(&viewer);
1562:     CreatePartitionVec(dm, &dmCell, &partition);
1563:     OutputVTK(dmCell, "ex11-partition.vtk", &viewer);
1564:     VecView(partition, viewer);
1565:     PetscViewerDestroy(&viewer);
1566:     VecDestroy(&partition);
1567:     DMDestroy(&dmCell);
1568:   }

1570:   DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius);
1571:   DMDestroy(&plex);
1572:   TSSetMaxTime(ts, 2.0);
1573:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1574:   dt = cfl * minRadius / user->model->maxspeed;
1575:   TSSetTimeStep(ts, dt);
1576:   TSSetFromOptions(ts);
1577:   TSSolve(ts, X);
1578:   TSGetSolveTime(ts, &ftime);
1579:   TSGetStepNumber(ts, &nsteps);
1580:   TSGetConvergedReason(ts, &reason);
1581:   PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps);
1582:   TSDestroy(&ts);

1584:   PetscFunctionListDestroy(&PhysicsList);
1585:   FunctionalLinkDestroy(&user->model->functionalRegistry);
1586:   PetscFree(user->model->functionalMonitored);
1587:   PetscFree(user->model->functionalCall);
1588:   PetscFree(user->model->physics->data);
1589:   PetscFree(user->model->physics);
1590:   PetscFree(user->model);
1591:   PetscFree(user);
1592:   VecDestroy(&X);
1593:   PetscFVDestroy(&fvm);
1594:   DMDestroy(&dm);
1595:   PetscFinalize();
1596:   return 0;
1597: }