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:
7:   f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
8:
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: };

72:   char              *name;
73:   FunctionalFunction func;
74:   void              *ctx;
75:   PetscInt           offset;
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 communication, but some error conditions can be collective */
90:   Physics          physics;
92:   PetscInt         maxComputed;
93:   PetscInt         numMonitored;
95:   PetscInt         numCall;
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: }

190: typedef enum {
195: typedef enum {

201: typedef struct {
202:   PetscReal wind[DIM];
204: typedef struct {
205:   PetscReal         center[DIM];

210: typedef struct {
211:   PetscReal     inflowState;
213:   union
214:   {
217:   } sol;
218:   struct {
219:     PetscInt Error;
220:   } functional;

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;

233:   PetscFunctionBeginUser;
235:   PetscFunctionReturn(PETSC_SUCCESS);
236: }

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

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

250:   PetscFunctionBeginUser;
254:     wind[0]                       = tilted->wind[0];
255:     wind[1]                       = tilted->wind[1];
256:   } break;
258:     wind[0] = -qp[1];
259:     wind[1] = qp[0];
260:     break;
261:   default:
263:   }
264:   wn      = Dot2(wind, n);
265:   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
266:   PetscFunctionReturn(PETSC_SUCCESS);
267: }

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

274:   PetscFunctionBeginUser;
277:     PetscReal              x0[DIM];
279:     Waxpy2(-time, tilted->wind, x, x0);
280:     if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
282:   } break;
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) {
294:       u[0] = PetscMax(1 - r / bump->radius, 0);
295:       break;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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;
311:   PetscScalar     yexact[1];

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

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

323:   PetscFunctionBeginUser;
329:   {
330:     PetscInt two = 2, dof = 1;
336:       two                           = 2;
337:       tilted->wind[0]               = 0.0;
338:       tilted->wind[1]               = 1.0;
339:       PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
342:       phys->maxspeed = Norm2(tilted->wind);
343:     } break;
346:       two                       = 2;
347:       bump->center[0]           = 2.;
348:       bump->center[1]           = 0.;
349:       PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
354:       phys->maxspeed = 3.; /* radius of mesh, kludge */
355:     } break;
356:     }
357:   }
359:   {
360:     const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
361:     DMLabel        label;

363:     PetscCall(DMGetLabel(dm, "Face Sets", &label));
364:     /* Register "canned" boundary conditions and defaults for where to apply. */
365:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Inflow, NULL, phys, NULL));
366:     PetscCall(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 */
369:     /* Register "canned" functionals */
371:   }
372:   PetscFunctionReturn(PETSC_SUCCESS);
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;

408:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
414: }

416: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
417: {
418:   PetscFunctionBeginUser;
419:   xG[0] = xI[0];
420:   xG[1] = -xI[1];
421:   xG[2] = -xI[2];
422:   PetscFunctionReturn(PETSC_SUCCESS);
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;

433:   PetscFunctionBeginUser;
434:   PetscCheck(uL->h >= 0 && uR->h >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
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:   PetscFunctionReturn(PETSC_SUCCESS);
445: }

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

451:   PetscFunctionBeginUser;
452:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
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:   PetscFunctionReturn(PETSC_SUCCESS);
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;

471:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
478: }

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

484:   PetscFunctionBeginUser;
485:   phys->field_desc = PhysicsFields_SW;
486:   phys->riemann    = (RiemannFunction)PhysicsRiemann_SW;
487:   PetscCall(PetscNew(&sw));
488:   phys->data = sw;
490:   {
491:     sw->gravity = 1.0;
492:     PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
493:   }
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:     PetscCall(DMGetLabel(dm, "Face Sets", &label));
502:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_SW_Wall, NULL, phys, NULL));
503:     PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
504:     PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
505:     PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
506:     PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
507:   }
508:   PetscFunctionReturn(PETSC_SUCCESS);
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;

544:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
550: }

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

556:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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;

578:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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];

596:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
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;

613:   PetscFunctionBeginUser;
614:   PetscCheck(uL->r >= 0 && uR->r >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
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:   PetscFunctionReturn(PETSC_SUCCESS);
622: }

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

628:   PetscFunctionBeginUser;
629:   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
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:   PetscFunctionReturn(PETSC_SUCCESS);
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;

643:   PetscFunctionBeginUser;
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:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

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

657:   PetscFunctionBeginUser;
658:   phys->field_desc = PhysicsFields_Euler;
659:   phys->riemann    = (RiemannFunction)PhysicsRiemann_Euler_Rusanov;
660:   PetscCall(PetscNew(&eu));
661:   phys->data = eu;
663:   {
664:     eu->pars[0] = 3.0;
665:     eu->pars[1] = 1.67;
666:     PetscCall(PetscOptionsReal("-eu_f", "Degrees of freedom", "", eu->pars[0], &eu->pars[0], NULL));
667:     PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->pars[1], &eu->pars[1], NULL));
668:   }
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:     PetscCall(DMGetLabel(dm, "Face Sets", &label));
678:     PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
679:     PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
680:     PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
681:     PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
682:     PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
683:     PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
684:     PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
685:   }
686:   PetscFunctionReturn(PETSC_SUCCESS);
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;

699:   PetscFunctionBeginUser;
700:   PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
701:   PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
702:   PetscCall(DMPlexGetCellTypeStratum(dm, DM_POLYTOPE_FV_GHOST, &cEndInterior, NULL));

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

724:     PetscCheck((cell >= cStart) && (cell < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", cell);
725:     PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
726:     PetscCall(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;

732:       PetscCheck((face >= fStart) && (face < fEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a face", face);
733:       PetscCall(DMPlexGetSupportSize(dm, face, &nC));
734:       if (nC != 2) continue;
735:       PetscCall(DMPlexGetSupport(dm, face, &neighbors));
736:       if ((neighbors[0] >= cEndInterior) || (neighbors[1] >= cEndInterior)) continue;
737:       PetscCheck((neighbors[0] >= cStart) && (neighbors[0] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[0]);
738:       PetscCheck((neighbors[1] >= cStart) && (neighbors[1] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[1]);
739:       PetscCall(DMGetLabelValue(dm, name, neighbors[0], &regionA));
740:       PetscCall(DMGetLabelValue(dm, name, neighbors[1], &regionB));
741:       PetscCheck(regionA >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[0]);
742:       PetscCheck(regionB >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[1]);
743:       if (regionA != regionB) PetscCall(DMSetLabelValue(dm, bdname, faces[f], 1));
744:     }
745:   }
746:   PetscCall(ISRestoreIndices(innerIS, &cells));
747:   PetscCall(ISDestroy(&innerIS));
748:   {
749:     DMLabel label;

751:     PetscCall(DMGetLabel(dm, bdname, &label));
752:     PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
753:   }
754:   PetscFunctionReturn(PETSC_SUCCESS);
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;

774:   PetscFunctionBeginUser;
775:   PetscCall(DMHasLabel(dm, labelName, &hasLabel));
776:   if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
777:   PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &sdm));
778:   PetscCall(DMSetType(sdm, DMPLEX));
779:   PetscCall(DMGetDimension(dm, &dim));
780:   PetscCall(DMSetDimension(sdm, dim));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1052:         sides[1] = faces[g];
1054:         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
1055:         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
1057:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
1058:         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
1059:         m[f * numFaces + g] = Dot2(fgA->normal, fgB->normal) * area * 0.5;
1060:         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
1061:       }
1062:     }
1063:   }
1067:   PetscCall(VecRestoreArray(*massMatrix, &m));
1068:   PetscCall(DMDestroy(&dmMass));
1069:   PetscCall(DMDestroy(&plex));
1070:   PetscFunctionReturn(PETSC_SUCCESS);
1071: }

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

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

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

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

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

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

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

1156:   PetscFunctionBeginUser;
1157:   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1160:   link->offset = lastoffset + 1;
1166:   PetscFunctionReturn(PETSC_SUCCESS);
1167: }

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

1175:   PetscFunctionBeginUser;
1176:   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
1177:   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
1178:   /* Create list of functionals that will be computed somehow */
1179:   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
1180:   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1181:   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
1182:   mod->numCall = 0;
1183:   for (i = 0; i < mod->numMonitored; i++) {
1185:       PetscBool match;
1187:       if (match) break;
1188:     }
1189:     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
1191:     for (j = 0; j < i; j++) {
1193:     }
1194:     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1195:   next_name:
1196:     PetscCall(PetscFree(names[i]));
1197:   }

1199:   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1200:   mod->maxComputed = -1;
1202:     for (i = 0; i < mod->numCall; i++) {
1205:     }
1206:   }
1207:   PetscFunctionReturn(PETSC_SUCCESS);
1208: }

1211: {

1214:   PetscFunctionBeginUser;
1218:   for (; l; l = next) {
1219:     next = l->next;
1220:     PetscCall(PetscFree(l->name));
1221:     PetscCall(PetscFree(l));
1222:   }
1223:   PetscFunctionReturn(PETSC_SUCCESS);
1224: }

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

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

1247:     PetscCall(DMPlexPointGlobalRef(dm, c, x, &xc));
1248:     if (xc) PetscCall((*mod->solution)(mod, 0.0, cg->centroid, xc, mod->solutionctx));
1249:   }
1251:   PetscCall(VecRestoreArray(X, &x));
1252:   PetscCall(DMDestroy(&plex));
1253:   PetscFunctionReturn(PETSC_SUCCESS);
1254: }

1256: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1257: {
1258:   PetscFunctionBeginUser;
1259:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1260:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
1261:   PetscCall(PetscViewerFileSetName(*viewer, filename));
1262:   PetscFunctionReturn(PETSC_SUCCESS);
1263: }

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

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

1286:     PetscCall(DMConvert(dm, DMPLEX, &plex));
1287:     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1288:     fcount = mod->maxComputed + 1;
1289:     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
1290:     for (i = 0; i < fcount; i++) {
1291:       fmin[i]      = PETSC_MAX_REAL;
1292:       fmax[i]      = PETSC_MIN_REAL;
1293:       fintegral[i] = 0;
1294:     }
1295:     PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1296:     PetscCall(VecGetDM(cellgeom, &dmCell));
1299:     for (c = cStart; c < cEnd; ++c) {
1300:       const PetscFVCellGeom *cg;
1301:       const PetscScalar     *cx;
1304:       if (!cx) continue; /* not a global cell */
1305:       for (i = 0; i < mod->numCall; i++) {
1308:       }
1309:       for (i = 0; i < fcount; i++) {
1310:         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1311:         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1312:         fintegral[i] += cg->volume * ftmp[i];
1313:       }
1314:     }
1317:     PetscCall(DMDestroy(&plex));
1318:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1319:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1320:     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));

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

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

1372: static PetscErrorCode OutputBIN(DM dm, const char *filename, PetscViewer *viewer)
1373: {
1374:   PetscFunctionBeginUser;
1375:   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1376:   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERBINARY));
1377:   PetscCall(PetscViewerFileSetMode(*viewer, FILE_MODE_WRITE));
1378:   PetscCall(PetscViewerFileSetName(*viewer, filename));
1379:   PetscFunctionReturn(PETSC_SUCCESS);
1380: }

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

1390:   PetscFunctionBeginUser;
1392:   PetscCall(VecCreate(PETSC_COMM_WORLD, &odesolution));
1394:   VecEqual(X, odesolution, &equal);
1395:   PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the vec data from file");
1396:   else
1397:   {
1398:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for Vec\n"));
1399:   }
1400:   /*Nr   = 1;

1404:   PetscCheck(timeread == time, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the scalar data from file");
1405:   else
1406:   {
1407:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for PetscReal\n"));
1408:   }

1410:   PetscCall(PetscViewerDestroy(&viewer));
1411:   PetscCall(VecDestroy(&odesolution));
1412:   PetscFunctionReturn(PETSC_SUCCESS);
1413: }

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

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

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

1453:   PetscFunctionBeginUser;
1454:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1455:   comm = PETSC_COMM_WORLD;
1456:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

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

1465:   /* Register physical models to be available on the command line */

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

1489:   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1490:   {
1491:     PetscDS prob;
1492:     PetscErrorCode (*physcreate)(PetscDS, Model, Physics);

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

1505:     PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1506:     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1507:     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1508:   }
1509:   PetscOptionsEnd();
1510:   {
1511:     DM dmDist;

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

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

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

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

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

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

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

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