xref: /petsc/src/ts/tutorials/ex11.c (revision 4d86920da9ee67c835173a5dfffa1b3a52fd24ca) 
1c4762a1bSJed Brown static char help[] = "Second Order TVD Finite Volume Example.\n";
2c4762a1bSJed Brown /*F
3c4762a1bSJed Brown 
4c4762a1bSJed Brown We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops
5c4762a1bSJed Brown over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values,
6c4762a1bSJed Brown \begin{equation}
7c4762a1bSJed Brown   f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
8c4762a1bSJed Brown \end{equation}
9c4762a1bSJed Brown and then update the cell values given the cell volume.
10c4762a1bSJed Brown \begin{eqnarray}
11c4762a1bSJed Brown     f^L_i &-=& \frac{f_i}{vol^L} \\
12c4762a1bSJed Brown     f^R_i &+=& \frac{f_i}{vol^R}
13c4762a1bSJed Brown \end{eqnarray}
14c4762a1bSJed Brown 
15c4762a1bSJed Brown As an example, we can consider the shallow water wave equation,
16c4762a1bSJed Brown \begin{eqnarray}
17c4762a1bSJed Brown      h_t + \nabla\cdot \left( uh                              \right) &=& 0 \\
18c4762a1bSJed Brown   (uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0
19c4762a1bSJed Brown \end{eqnarray}
20c4762a1bSJed Brown where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity.
21c4762a1bSJed Brown 
22c4762a1bSJed Brown A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
23c4762a1bSJed Brown \begin{eqnarray}
24c4762a1bSJed Brown   f^{L,R}_h    &=& uh^{L,R} \cdot \hat n \\
25c4762a1bSJed Brown   f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
26c4762a1bSJed Brown   c^{L,R}      &=& \sqrt{g h^{L,R}} \\
27c4762a1bSJed Brown   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) \\
28c4762a1bSJed Brown   f_i          &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
29c4762a1bSJed Brown \end{eqnarray}
30c4762a1bSJed Brown where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.
31c4762a1bSJed Brown 
32c4762a1bSJed Brown The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial
33c4762a1bSJed Brown over a neighborhood of the given element.
34c4762a1bSJed Brown 
35c4762a1bSJed Brown The mesh is read in from an ExodusII file, usually generated by Cubit.
3652115386SStefano Zampini 
3752115386SStefano Zampini The example also shows how to handle AMR in a time-dependent TS solver.
38c4762a1bSJed Brown F*/
39c4762a1bSJed Brown #include <petscdmplex.h>
40c4762a1bSJed Brown #include <petscdmforest.h>
41c4762a1bSJed Brown #include <petscds.h>
42c4762a1bSJed Brown #include <petscts.h>
43c4762a1bSJed Brown 
44de8de29fSMatthew G. Knepley #include "ex11.h"
45c4762a1bSJed Brown 
46de8de29fSMatthew G. Knepley static PetscFunctionList PhysicsList, PhysicsRiemannList_SW, PhysicsRiemannList_Euler;
47c4762a1bSJed Brown 
48c4762a1bSJed Brown /* 'User' implements a discretization of a continuous model. */
49c4762a1bSJed Brown typedef struct _n_User *User;
50c4762a1bSJed Brown typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
5145480ffeSMatthew G. Knepley typedef PetscErrorCode (*SetUpBCFunction)(DM, PetscDS, Physics);
52c4762a1bSJed Brown typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
53c4762a1bSJed Brown typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
54c4762a1bSJed Brown static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
55c4762a1bSJed Brown static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
56c4762a1bSJed Brown static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);
57c4762a1bSJed Brown 
58c4762a1bSJed Brown typedef struct _n_FunctionalLink *FunctionalLink;
59c4762a1bSJed Brown struct _n_FunctionalLink {
60c4762a1bSJed Brown   char              *name;
61c4762a1bSJed Brown   FunctionalFunction func;
62c4762a1bSJed Brown   void              *ctx;
63c4762a1bSJed Brown   PetscInt           offset;
64c4762a1bSJed Brown   FunctionalLink     next;
65c4762a1bSJed Brown };
66c4762a1bSJed Brown 
67c4762a1bSJed Brown struct _n_Model {
68da81f932SPierre Jolivet   MPI_Comm         comm; /* Does not do collective communication, but some error conditions can be collective */
69c4762a1bSJed Brown   Physics          physics;
70c4762a1bSJed Brown   FunctionalLink   functionalRegistry;
71c4762a1bSJed Brown   PetscInt         maxComputed;
72c4762a1bSJed Brown   PetscInt         numMonitored;
73c4762a1bSJed Brown   FunctionalLink  *functionalMonitored;
74c4762a1bSJed Brown   PetscInt         numCall;
75c4762a1bSJed Brown   FunctionalLink  *functionalCall;
76c4762a1bSJed Brown   SolutionFunction solution;
77c4762a1bSJed Brown   SetUpBCFunction  setupbc;
78c4762a1bSJed Brown   void            *solutionctx;
79c4762a1bSJed Brown   PetscReal        maxspeed; /* estimate of global maximum speed (for CFL calculation) */
80c4762a1bSJed Brown   PetscReal        bounds[2 * DIM];
81c4762a1bSJed Brown   PetscErrorCode (*errorIndicator)(PetscInt, PetscReal, PetscInt, const PetscScalar[], const PetscScalar[], PetscReal *, void *);
82c4762a1bSJed Brown   void *errorCtx;
83de8de29fSMatthew G. Knepley   PetscErrorCode (*setupCEED)(DM, Physics);
84c4762a1bSJed Brown };
85c4762a1bSJed Brown 
86c4762a1bSJed Brown struct _n_User {
87c4762a1bSJed Brown   PetscInt  vtkInterval;                        /* For monitor */
88c4762a1bSJed Brown   char      outputBasename[PETSC_MAX_PATH_LEN]; /* Basename for output files */
89c4762a1bSJed Brown   PetscInt  monitorStepOffset;
90c4762a1bSJed Brown   Model     model;
91c4762a1bSJed Brown   PetscBool vtkmon;
92c4762a1bSJed Brown };
93c4762a1bSJed Brown 
94de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
95de8de29fSMatthew G. Knepley // Free a plain data context that was allocated using PETSc; returning libCEED error codes
96de8de29fSMatthew G. Knepley static int FreeContextPetsc(void *data)
97d71ae5a4SJacob Faibussowitsch {
98de8de29fSMatthew G. Knepley   if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed");
99de8de29fSMatthew G. Knepley   return CEED_ERROR_SUCCESS;
100c4762a1bSJed Brown }
101de8de29fSMatthew G. Knepley #endif
102c4762a1bSJed Brown 
103c4762a1bSJed Brown /******************* Advect ********************/
1049371c9d4SSatish Balay typedef enum {
1059371c9d4SSatish Balay   ADVECT_SOL_TILTED,
1069371c9d4SSatish Balay   ADVECT_SOL_BUMP,
1079371c9d4SSatish Balay   ADVECT_SOL_BUMP_CAVITY
1089371c9d4SSatish Balay } AdvectSolType;
109c4762a1bSJed Brown static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "BUMP_CAVITY", "AdvectSolType", "ADVECT_SOL_", 0};
1109371c9d4SSatish Balay typedef enum {
1119371c9d4SSatish Balay   ADVECT_SOL_BUMP_CONE,
1129371c9d4SSatish Balay   ADVECT_SOL_BUMP_COS
1139371c9d4SSatish Balay } AdvectSolBumpType;
114c4762a1bSJed Brown static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};
115c4762a1bSJed Brown 
116c4762a1bSJed Brown typedef struct {
117c4762a1bSJed Brown   PetscReal wind[DIM];
118c4762a1bSJed Brown } Physics_Advect_Tilted;
119c4762a1bSJed Brown typedef struct {
120c4762a1bSJed Brown   PetscReal         center[DIM];
121c4762a1bSJed Brown   PetscReal         radius;
122c4762a1bSJed Brown   AdvectSolBumpType type;
123c4762a1bSJed Brown } Physics_Advect_Bump;
124c4762a1bSJed Brown 
125c4762a1bSJed Brown typedef struct {
126c4762a1bSJed Brown   PetscReal     inflowState;
127c4762a1bSJed Brown   AdvectSolType soltype;
1289371c9d4SSatish Balay   union
1299371c9d4SSatish Balay   {
130c4762a1bSJed Brown     Physics_Advect_Tilted tilted;
131c4762a1bSJed Brown     Physics_Advect_Bump   bump;
132c4762a1bSJed Brown   } sol;
133c4762a1bSJed Brown   struct {
134c4762a1bSJed Brown     PetscInt Solution;
135c4762a1bSJed Brown     PetscInt Error;
136c4762a1bSJed Brown   } functional;
137c4762a1bSJed Brown } Physics_Advect;
138c4762a1bSJed Brown 
1399371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_Advect[] = {
1409371c9d4SSatish Balay   {"U",  1},
1419371c9d4SSatish Balay   {NULL, 0}
1429371c9d4SSatish Balay };
143c4762a1bSJed Brown 
144d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
145d71ae5a4SJacob Faibussowitsch {
146c4762a1bSJed Brown   Physics         phys   = (Physics)ctx;
147c4762a1bSJed Brown   Physics_Advect *advect = (Physics_Advect *)phys->data;
148c4762a1bSJed Brown 
149c4762a1bSJed Brown   PetscFunctionBeginUser;
150c4762a1bSJed Brown   xG[0] = advect->inflowState;
1513ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
152c4762a1bSJed Brown }
153c4762a1bSJed Brown 
154d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
155d71ae5a4SJacob Faibussowitsch {
156c4762a1bSJed Brown   PetscFunctionBeginUser;
157c4762a1bSJed Brown   xG[0] = xI[0];
1583ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
159c4762a1bSJed Brown }
160c4762a1bSJed Brown 
161d71ae5a4SJacob Faibussowitsch static void PhysicsRiemann_Advect(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
162d71ae5a4SJacob Faibussowitsch {
163c4762a1bSJed Brown   Physics_Advect *advect = (Physics_Advect *)phys->data;
164c4762a1bSJed Brown   PetscReal       wind[DIM], wn;
165c4762a1bSJed Brown 
166c4762a1bSJed Brown   switch (advect->soltype) {
167c4762a1bSJed Brown   case ADVECT_SOL_TILTED: {
168c4762a1bSJed Brown     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
169c4762a1bSJed Brown     wind[0]                       = tilted->wind[0];
170c4762a1bSJed Brown     wind[1]                       = tilted->wind[1];
171c4762a1bSJed Brown   } break;
172c4762a1bSJed Brown   case ADVECT_SOL_BUMP:
173c4762a1bSJed Brown     wind[0] = -qp[1];
174c4762a1bSJed Brown     wind[1] = qp[0];
175c4762a1bSJed Brown     break;
1769371c9d4SSatish Balay   case ADVECT_SOL_BUMP_CAVITY: {
177c4762a1bSJed Brown     PetscInt  i;
178c4762a1bSJed Brown     PetscReal comp2[3] = {0., 0., 0.}, rad2;
179c4762a1bSJed Brown 
180c4762a1bSJed Brown     rad2 = 0.;
181c4762a1bSJed Brown     for (i = 0; i < dim; i++) {
182c4762a1bSJed Brown       comp2[i] = qp[i] * qp[i];
183c4762a1bSJed Brown       rad2 += comp2[i];
184c4762a1bSJed Brown     }
185c4762a1bSJed Brown 
186c4762a1bSJed Brown     wind[0] = -qp[1];
187c4762a1bSJed Brown     wind[1] = qp[0];
188c4762a1bSJed Brown     if (rad2 > 1.) {
189c4762a1bSJed Brown       PetscInt  maxI     = 0;
190c4762a1bSJed Brown       PetscReal maxComp2 = comp2[0];
191c4762a1bSJed Brown 
192c4762a1bSJed Brown       for (i = 1; i < dim; i++) {
193c4762a1bSJed Brown         if (comp2[i] > maxComp2) {
194c4762a1bSJed Brown           maxI     = i;
195c4762a1bSJed Brown           maxComp2 = comp2[i];
196c4762a1bSJed Brown         }
197c4762a1bSJed Brown       }
198c4762a1bSJed Brown       wind[maxI] = 0.;
199c4762a1bSJed Brown     }
2009371c9d4SSatish Balay   } break;
2019371c9d4SSatish Balay   default: {
202c4762a1bSJed Brown     PetscInt i;
203c4762a1bSJed Brown     for (i = 0; i < DIM; ++i) wind[i] = 0.0;
204c4762a1bSJed Brown   }
20598921bdaSJacob Faibussowitsch     /* default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
206c4762a1bSJed Brown   }
207c4762a1bSJed Brown   wn      = Dot2Real(wind, n);
208c4762a1bSJed Brown   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
209c4762a1bSJed Brown }
210c4762a1bSJed Brown 
211d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
212d71ae5a4SJacob Faibussowitsch {
213c4762a1bSJed Brown   Physics         phys   = (Physics)ctx;
214c4762a1bSJed Brown   Physics_Advect *advect = (Physics_Advect *)phys->data;
215c4762a1bSJed Brown 
216c4762a1bSJed Brown   PetscFunctionBeginUser;
217c4762a1bSJed Brown   switch (advect->soltype) {
218c4762a1bSJed Brown   case ADVECT_SOL_TILTED: {
219c4762a1bSJed Brown     PetscReal              x0[DIM];
220c4762a1bSJed Brown     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
221c4762a1bSJed Brown     Waxpy2Real(-time, tilted->wind, x, x0);
222c4762a1bSJed Brown     if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
223c4762a1bSJed Brown     else u[0] = advect->inflowState;
224c4762a1bSJed Brown   } break;
225c4762a1bSJed Brown   case ADVECT_SOL_BUMP_CAVITY:
226c4762a1bSJed Brown   case ADVECT_SOL_BUMP: {
227c4762a1bSJed Brown     Physics_Advect_Bump *bump = &advect->sol.bump;
228c4762a1bSJed Brown     PetscReal            x0[DIM], v[DIM], r, cost, sint;
229c4762a1bSJed Brown     cost  = PetscCosReal(time);
230c4762a1bSJed Brown     sint  = PetscSinReal(time);
231c4762a1bSJed Brown     x0[0] = cost * x[0] + sint * x[1];
232c4762a1bSJed Brown     x0[1] = -sint * x[0] + cost * x[1];
233c4762a1bSJed Brown     Waxpy2Real(-1, bump->center, x0, v);
234c4762a1bSJed Brown     r = Norm2Real(v);
235c4762a1bSJed Brown     switch (bump->type) {
236d71ae5a4SJacob Faibussowitsch     case ADVECT_SOL_BUMP_CONE:
237d71ae5a4SJacob Faibussowitsch       u[0] = PetscMax(1 - r / bump->radius, 0);
238d71ae5a4SJacob Faibussowitsch       break;
239d71ae5a4SJacob Faibussowitsch     case ADVECT_SOL_BUMP_COS:
240d71ae5a4SJacob Faibussowitsch       u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
241d71ae5a4SJacob Faibussowitsch       break;
242c4762a1bSJed Brown     }
243c4762a1bSJed Brown   } break;
244d71ae5a4SJacob Faibussowitsch   default:
245d71ae5a4SJacob Faibussowitsch     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
246c4762a1bSJed Brown   }
2473ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
248c4762a1bSJed Brown }
249c4762a1bSJed Brown 
250d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx)
251d71ae5a4SJacob Faibussowitsch {
252c4762a1bSJed Brown   Physics         phys      = (Physics)ctx;
253c4762a1bSJed Brown   Physics_Advect *advect    = (Physics_Advect *)phys->data;
254391faea4SSatish Balay   PetscScalar     yexact[1] = {0.0};
255c4762a1bSJed Brown 
256c4762a1bSJed Brown   PetscFunctionBeginUser;
2579566063dSJacob Faibussowitsch   PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
258c4762a1bSJed Brown   f[advect->functional.Solution] = PetscRealPart(y[0]);
259c4762a1bSJed Brown   f[advect->functional.Error]    = PetscAbsScalar(y[0] - yexact[0]);
2603ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
261c4762a1bSJed Brown }
262c4762a1bSJed Brown 
263d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_Advect(DM dm, PetscDS prob, Physics phys)
264d71ae5a4SJacob Faibussowitsch {
265c4762a1bSJed Brown   const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
26645480ffeSMatthew G. Knepley   DMLabel        label;
267c4762a1bSJed Brown 
268c4762a1bSJed Brown   PetscFunctionBeginUser;
269c4762a1bSJed Brown   /* Register "canned" boundary conditions and defaults for where to apply. */
2709566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
271dd39110bSPierre Jolivet   PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Inflow, NULL, phys, NULL));
272dd39110bSPierre Jolivet   PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Outflow, NULL, phys, NULL));
2733ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
274c4762a1bSJed Brown }
275c4762a1bSJed Brown 
276d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_Advect(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
277d71ae5a4SJacob Faibussowitsch {
278c4762a1bSJed Brown   Physics_Advect *advect;
279c4762a1bSJed Brown 
280c4762a1bSJed Brown   PetscFunctionBeginUser;
281c4762a1bSJed Brown   phys->field_desc = PhysicsFields_Advect;
282c4762a1bSJed Brown   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Advect;
2839566063dSJacob Faibussowitsch   PetscCall(PetscNew(&advect));
284c4762a1bSJed Brown   phys->data   = advect;
285c4762a1bSJed Brown   mod->setupbc = SetUpBC_Advect;
286c4762a1bSJed Brown 
287d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
288c4762a1bSJed Brown   {
289c4762a1bSJed Brown     PetscInt two = 2, dof = 1;
290c4762a1bSJed Brown     advect->soltype = ADVECT_SOL_TILTED;
2919566063dSJacob Faibussowitsch     PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL));
292c4762a1bSJed Brown     switch (advect->soltype) {
293c4762a1bSJed Brown     case ADVECT_SOL_TILTED: {
294c4762a1bSJed Brown       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
295c4762a1bSJed Brown       two                           = 2;
296c4762a1bSJed Brown       tilted->wind[0]               = 0.0;
297c4762a1bSJed Brown       tilted->wind[1]               = 1.0;
2989566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
299c4762a1bSJed Brown       advect->inflowState = -2.0;
3009566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL));
301c4762a1bSJed Brown       phys->maxspeed = Norm2Real(tilted->wind);
302c4762a1bSJed Brown     } break;
303c4762a1bSJed Brown     case ADVECT_SOL_BUMP_CAVITY:
304c4762a1bSJed Brown     case ADVECT_SOL_BUMP: {
305c4762a1bSJed Brown       Physics_Advect_Bump *bump = &advect->sol.bump;
306c4762a1bSJed Brown       two                       = 2;
307c4762a1bSJed Brown       bump->center[0]           = 2.;
308c4762a1bSJed Brown       bump->center[1]           = 0.;
3099566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
310c4762a1bSJed Brown       bump->radius = 0.9;
3119566063dSJacob Faibussowitsch       PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL));
312c4762a1bSJed Brown       bump->type = ADVECT_SOL_BUMP_CONE;
3139566063dSJacob Faibussowitsch       PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL));
314c4762a1bSJed Brown       phys->maxspeed = 3.; /* radius of mesh, kludge */
315c4762a1bSJed Brown     } break;
316c4762a1bSJed Brown     }
317c4762a1bSJed Brown   }
318d0609cedSBarry Smith   PetscOptionsHeadEnd();
319c4762a1bSJed Brown   /* Initial/transient solution with default boundary conditions */
3209566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys));
321c4762a1bSJed Brown   /* Register "canned" functionals */
3229566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Solution", &advect->functional.Solution, PhysicsFunctional_Advect, phys));
3239566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys));
3243ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
325c4762a1bSJed Brown }
326c4762a1bSJed Brown 
327c4762a1bSJed Brown /******************* Shallow Water ********************/
3289371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_SW[] = {
3299371c9d4SSatish Balay   {"Height",   1  },
3309371c9d4SSatish Balay   {"Momentum", DIM},
3319371c9d4SSatish Balay   {NULL,       0  }
3329371c9d4SSatish Balay };
333c4762a1bSJed Brown 
334d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
335d71ae5a4SJacob Faibussowitsch {
336c4762a1bSJed Brown   PetscFunctionBeginUser;
337c4762a1bSJed Brown   xG[0] = xI[0];
338c4762a1bSJed Brown   xG[1] = -xI[1];
339c4762a1bSJed Brown   xG[2] = -xI[2];
3403ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
341c4762a1bSJed Brown }
342c4762a1bSJed Brown 
343d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_SW(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
344d71ae5a4SJacob Faibussowitsch {
345c4762a1bSJed Brown   PetscReal dx[2], r, sigma;
346c4762a1bSJed Brown 
347c4762a1bSJed Brown   PetscFunctionBeginUser;
34854c59aa7SJacob Faibussowitsch   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
349c4762a1bSJed Brown   dx[0] = x[0] - 1.5;
350c4762a1bSJed Brown   dx[1] = x[1] - 1.0;
351c4762a1bSJed Brown   r     = Norm2Real(dx);
352c4762a1bSJed Brown   sigma = 0.5;
353c4762a1bSJed Brown   u[0]  = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma)));
354c4762a1bSJed Brown   u[1]  = 0.0;
355c4762a1bSJed Brown   u[2]  = 0.0;
3563ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
357c4762a1bSJed Brown }
358c4762a1bSJed Brown 
359d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
360d71ae5a4SJacob Faibussowitsch {
361c4762a1bSJed Brown   Physics       phys = (Physics)ctx;
362c4762a1bSJed Brown   Physics_SW   *sw   = (Physics_SW *)phys->data;
363c4762a1bSJed Brown   const SWNode *x    = (const SWNode *)xx;
364c4762a1bSJed Brown   PetscReal     u[2];
365c4762a1bSJed Brown   PetscReal     h;
366c4762a1bSJed Brown 
367c4762a1bSJed Brown   PetscFunctionBeginUser;
368c4762a1bSJed Brown   h = x->h;
369c4762a1bSJed Brown   Scale2Real(1. / x->h, x->uh, u);
370c4762a1bSJed Brown   f[sw->functional.Height] = h;
371c4762a1bSJed Brown   f[sw->functional.Speed]  = Norm2Real(u) + PetscSqrtReal(sw->gravity * h);
372c4762a1bSJed Brown   f[sw->functional.Energy] = 0.5 * (Dot2Real(x->uh, u) + sw->gravity * PetscSqr(h));
3733ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
374c4762a1bSJed Brown }
375c4762a1bSJed Brown 
376d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_SW(DM dm, PetscDS prob, Physics phys)
377d71ae5a4SJacob Faibussowitsch {
378c4762a1bSJed Brown   const PetscInt wallids[] = {100, 101, 200, 300};
37945480ffeSMatthew G. Knepley   DMLabel        label;
38045480ffeSMatthew G. Knepley 
381c4762a1bSJed Brown   PetscFunctionBeginUser;
3829566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
383dd39110bSPierre Jolivet   PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_SW_Wall, NULL, phys, NULL));
3843ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
385c4762a1bSJed Brown }
386c4762a1bSJed Brown 
387de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
388de8de29fSMatthew G. Knepley static PetscErrorCode CreateQFunctionContext_SW(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
389de8de29fSMatthew G. Knepley {
390de8de29fSMatthew G. Knepley   Physics_SW *in = (Physics_SW *)phys->data;
391de8de29fSMatthew G. Knepley   Physics_SW *sw;
392de8de29fSMatthew G. Knepley 
393de8de29fSMatthew G. Knepley   PetscFunctionBeginUser;
394de8de29fSMatthew G. Knepley   PetscCall(PetscCalloc1(1, &sw));
395de8de29fSMatthew G. Knepley 
396de8de29fSMatthew G. Knepley   sw->gravity = in->gravity;
397de8de29fSMatthew G. Knepley 
398de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
399de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*sw), sw));
400de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
401de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gravity", offsetof(Physics_SW, gravity), 1, "Accelaration due to gravity"));
402de8de29fSMatthew G. Knepley   PetscFunctionReturn(0);
403de8de29fSMatthew G. Knepley }
404de8de29fSMatthew G. Knepley #endif
405de8de29fSMatthew G. Knepley 
406de8de29fSMatthew G. Knepley static PetscErrorCode SetupCEED_SW(DM dm, Physics physics)
407de8de29fSMatthew G. Knepley {
408de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
409de8de29fSMatthew G. Knepley   Ceed                 ceed;
410de8de29fSMatthew G. Knepley   CeedQFunctionContext qfCtx;
411de8de29fSMatthew G. Knepley #endif
412de8de29fSMatthew G. Knepley 
413de8de29fSMatthew G. Knepley   PetscFunctionBegin;
414de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
415de8de29fSMatthew G. Knepley   PetscCall(DMGetCeed(dm, &ceed));
416de8de29fSMatthew G. Knepley   PetscCall(CreateQFunctionContext_SW(physics, ceed, &qfCtx));
417de8de29fSMatthew G. Knepley   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_SW_Rusanov_CEED, PhysicsRiemann_SW_Rusanov_CEED_loc, qfCtx));
418de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
419de8de29fSMatthew G. Knepley #endif
420de8de29fSMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
421de8de29fSMatthew G. Knepley }
422de8de29fSMatthew G. Knepley 
423d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
424d71ae5a4SJacob Faibussowitsch {
425c4762a1bSJed Brown   Physics_SW *sw;
4268d2c18caSMukkund Sunjii   char        sw_riemann[64] = "rusanov";
427c4762a1bSJed Brown 
428c4762a1bSJed Brown   PetscFunctionBeginUser;
429c4762a1bSJed Brown   phys->field_desc = PhysicsFields_SW;
4309566063dSJacob Faibussowitsch   PetscCall(PetscNew(&sw));
431c4762a1bSJed Brown   phys->data     = sw;
432c4762a1bSJed Brown   mod->setupbc   = SetUpBC_SW;
433de8de29fSMatthew G. Knepley   mod->setupCEED = SetupCEED_SW;
434c4762a1bSJed Brown 
4353ba16761SJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov));
4363ba16761SJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL));
437de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
438de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov_ceed", PhysicsRiemann_SW_Rusanov_CEED));
439de8de29fSMatthew G. Knepley #endif
4408d2c18caSMukkund Sunjii 
441d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
442c4762a1bSJed Brown   {
4438d2c18caSMukkund Sunjii     void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
444c4762a1bSJed Brown     sw->gravity = 1.0;
4459566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
4469566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL));
4479566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW));
4488d2c18caSMukkund Sunjii     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_SW;
449c4762a1bSJed Brown   }
450d0609cedSBarry Smith   PetscOptionsHeadEnd();
451c4762a1bSJed Brown   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */
452c4762a1bSJed Brown 
4539566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
4549566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
4559566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
4569566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
4573ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
458c4762a1bSJed Brown }
459c4762a1bSJed Brown 
460c4762a1bSJed Brown /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
461c4762a1bSJed Brown /* An initial-value and self-similar solutions of the compressible Euler equations */
462c4762a1bSJed Brown /* Ravi Samtaney and D. I. Pullin */
463c4762a1bSJed Brown /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
4649371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_Euler[] = {
4659371c9d4SSatish Balay   {"Density",  1  },
4669371c9d4SSatish Balay   {"Momentum", DIM},
4679371c9d4SSatish Balay   {"Energy",   1  },
4689371c9d4SSatish Balay   {NULL,       0  }
4699371c9d4SSatish Balay };
470c4762a1bSJed Brown 
471c4762a1bSJed Brown /* initial condition */
472c4762a1bSJed Brown int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
473d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
474d71ae5a4SJacob Faibussowitsch {
475c4762a1bSJed Brown   PetscInt       i;
476c4762a1bSJed Brown   Physics        phys = (Physics)ctx;
477c4762a1bSJed Brown   Physics_Euler *eu   = (Physics_Euler *)phys->data;
478c4762a1bSJed Brown   EulerNode     *uu   = (EulerNode *)u;
479de8de29fSMatthew G. Knepley   PetscReal      p0, gamma, c = 0.;
480*4d86920dSPierre Jolivet 
481c4762a1bSJed Brown   PetscFunctionBeginUser;
48254c59aa7SJacob Faibussowitsch   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
483c4762a1bSJed Brown 
484c4762a1bSJed Brown   for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
485c4762a1bSJed Brown   /* set E and rho */
486de8de29fSMatthew G. Knepley   gamma = eu->gamma;
487c4762a1bSJed Brown 
488c4762a1bSJed Brown   if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) {
489c4762a1bSJed Brown     /******************* Euler Density Shock ********************/
490c4762a1bSJed Brown     /* On initial-value and self-similar solutions of the compressible Euler equations */
491c4762a1bSJed Brown     /* Ravi Samtaney and D. I. Pullin */
492c4762a1bSJed Brown     /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
493c4762a1bSJed Brown     /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity,  */
494c4762a1bSJed Brown     p0 = 1.;
495de8de29fSMatthew G. Knepley     if (x[0] < 0.0 + x[1] * eu->itana) {
496c4762a1bSJed Brown       if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */
497c4762a1bSJed Brown         PetscReal amach, rho, press, gas1, p1;
498de8de29fSMatthew G. Knepley         amach     = eu->amach;
499c4762a1bSJed Brown         rho       = 1.;
500c4762a1bSJed Brown         press     = p0;
501c4762a1bSJed Brown         p1        = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0));
502c4762a1bSJed Brown         gas1      = (gamma - 1.0) / (gamma + 1.0);
503c4762a1bSJed Brown         uu->r     = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0);
504c4762a1bSJed Brown         uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach);
505c4762a1bSJed Brown         uu->E     = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0];
5069371c9d4SSatish Balay       } else {      /* left of discontinuity (0) */
507c4762a1bSJed Brown         uu->r = 1.; /* rho = 1 */
508c4762a1bSJed Brown         uu->E = p0 / (gamma - 1.0);
509c4762a1bSJed Brown       }
5109371c9d4SSatish Balay     } else { /* right of discontinuity (2) */
511de8de29fSMatthew G. Knepley       uu->r = eu->rhoR;
512c4762a1bSJed Brown       uu->E = p0 / (gamma - 1.0);
513c4762a1bSJed Brown     }
5149371c9d4SSatish Balay   } else if (eu->type == EULER_SHOCK_TUBE) {
515c4762a1bSJed Brown     /* For (x<x0) set (rho,u,p)=(8,0,10) and for (x>x0) set (rho,u,p)=(1,0,1). Choose x0 to the midpoint of the domain in the x-direction. */
516c4762a1bSJed Brown     if (x[0] < 0.0) {
517c4762a1bSJed Brown       uu->r = 8.;
518c4762a1bSJed Brown       uu->E = 10. / (gamma - 1.);
5199371c9d4SSatish Balay     } else {
520c4762a1bSJed Brown       uu->r = 1.;
521c4762a1bSJed Brown       uu->E = 1. / (gamma - 1.);
522c4762a1bSJed Brown     }
5239371c9d4SSatish Balay   } else if (eu->type == EULER_LINEAR_WAVE) {
524c4762a1bSJed Brown     initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
5259371c9d4SSatish Balay   } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type);
526c4762a1bSJed Brown 
527c4762a1bSJed Brown   /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */
528de8de29fSMatthew G. Knepley   PetscCall(SpeedOfSound_PG(gamma, uu, &c));
529c4762a1bSJed Brown   c = (uu->ru[0] / uu->r) + c;
530c4762a1bSJed Brown   if (c > phys->maxspeed) phys->maxspeed = c;
5313ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
532c4762a1bSJed Brown }
533c4762a1bSJed Brown 
534c4762a1bSJed Brown /* PetscReal* => EulerNode* conversion */
535d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
536d71ae5a4SJacob Faibussowitsch {
537c4762a1bSJed Brown   PetscInt         i;
538c4762a1bSJed Brown   const EulerNode *xI   = (const EulerNode *)a_xI;
539c4762a1bSJed Brown   EulerNode       *xG   = (EulerNode *)a_xG;
540c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
541c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
542*4d86920dSPierre Jolivet 
543c4762a1bSJed Brown   PetscFunctionBeginUser;
544c4762a1bSJed Brown   xG->r = xI->r;                                     /* ghost cell density - same */
545c4762a1bSJed Brown   xG->E = xI->E;                                     /* ghost cell energy - same */
546c4762a1bSJed Brown   if (n[1] != 0.) {                                  /* top and bottom */
547c4762a1bSJed Brown     xG->ru[0] = xI->ru[0];                           /* copy tang to wall */
548c4762a1bSJed Brown     xG->ru[1] = -xI->ru[1];                          /* reflect perp to t/b wall */
5499371c9d4SSatish Balay   } else {                                           /* sides */
550c4762a1bSJed Brown     for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
551c4762a1bSJed Brown   }
552c4762a1bSJed Brown   if (eu->type == EULER_LINEAR_WAVE) { /* debug */
553c4762a1bSJed Brown #if 0
55463a3b9bcSJacob Faibussowitsch     PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]);
555c4762a1bSJed Brown #endif
556c4762a1bSJed Brown   }
5573ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
558c4762a1bSJed Brown }
559c4762a1bSJed Brown 
560d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
561d71ae5a4SJacob Faibussowitsch {
562c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
563c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
564c4762a1bSJed Brown   const EulerNode *x    = (const EulerNode *)xx;
565c4762a1bSJed Brown   PetscReal        p;
566c4762a1bSJed Brown 
567c4762a1bSJed Brown   PetscFunctionBeginUser;
568c4762a1bSJed Brown   f[eu->monitor.Density]  = x->r;
569c4762a1bSJed Brown   f[eu->monitor.Momentum] = NormDIM(x->ru);
570c4762a1bSJed Brown   f[eu->monitor.Energy]   = x->E;
571c4762a1bSJed Brown   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
572de8de29fSMatthew G. Knepley   PetscCall(Pressure_PG(eu->gamma, x, &p));
573c4762a1bSJed Brown   f[eu->monitor.Pressure] = p;
5743ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
575c4762a1bSJed Brown }
576c4762a1bSJed Brown 
577d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys)
578d71ae5a4SJacob Faibussowitsch {
579c4762a1bSJed Brown   Physics_Euler *eu = (Physics_Euler *)phys->data;
58045480ffeSMatthew G. Knepley   DMLabel        label;
58145480ffeSMatthew G. Knepley 
582362febeeSStefano Zampini   PetscFunctionBeginUser;
5839566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
584c4762a1bSJed Brown   if (eu->type == EULER_LINEAR_WAVE) {
585c4762a1bSJed Brown     const PetscInt wallids[] = {100, 101};
586dd39110bSPierre Jolivet     PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
5879371c9d4SSatish Balay   } else {
588c4762a1bSJed Brown     const PetscInt wallids[] = {100, 101, 200, 300};
589dd39110bSPierre Jolivet     PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
590c4762a1bSJed Brown   }
5913ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
592c4762a1bSJed Brown }
593c4762a1bSJed Brown 
594de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
595de8de29fSMatthew G. Knepley static PetscErrorCode CreateQFunctionContext_Euler(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
596de8de29fSMatthew G. Knepley {
597de8de29fSMatthew G. Knepley   Physics_Euler *in = (Physics_Euler *)phys->data;
598de8de29fSMatthew G. Knepley   Physics_Euler *eu;
599de8de29fSMatthew G. Knepley 
600de8de29fSMatthew G. Knepley   PetscFunctionBeginUser;
601de8de29fSMatthew G. Knepley   PetscCall(PetscCalloc1(1, &eu));
602de8de29fSMatthew G. Knepley 
603de8de29fSMatthew G. Knepley   eu->gamma = in->gamma;
604de8de29fSMatthew G. Knepley 
605de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
606de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*eu), eu));
607de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
608de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gamma", offsetof(Physics_Euler, gamma), 1, "Heat capacity ratio"));
609de8de29fSMatthew G. Knepley   PetscFunctionReturn(0);
610de8de29fSMatthew G. Knepley }
611de8de29fSMatthew G. Knepley #endif
612de8de29fSMatthew G. Knepley 
613de8de29fSMatthew G. Knepley static PetscErrorCode SetupCEED_Euler(DM dm, Physics physics)
614de8de29fSMatthew G. Knepley {
615de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
616de8de29fSMatthew G. Knepley   Ceed                 ceed;
617de8de29fSMatthew G. Knepley   CeedQFunctionContext qfCtx;
618de8de29fSMatthew G. Knepley #endif
619de8de29fSMatthew G. Knepley 
620de8de29fSMatthew G. Knepley   PetscFunctionBegin;
621de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
622de8de29fSMatthew G. Knepley   PetscCall(DMGetCeed(dm, &ceed));
623de8de29fSMatthew G. Knepley   PetscCall(CreateQFunctionContext_Euler(physics, ceed, &qfCtx));
624de8de29fSMatthew G. Knepley   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_Euler_Godunov_CEED, PhysicsRiemann_Euler_Godunov_CEED_loc, qfCtx));
625de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
626de8de29fSMatthew G. Knepley #endif
627de8de29fSMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
628de8de29fSMatthew G. Knepley }
629de8de29fSMatthew G. Knepley 
630d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
631d71ae5a4SJacob Faibussowitsch {
632c4762a1bSJed Brown   Physics_Euler *eu;
633c4762a1bSJed Brown 
634c4762a1bSJed Brown   PetscFunctionBeginUser;
635c4762a1bSJed Brown   phys->field_desc = PhysicsFields_Euler;
636c4762a1bSJed Brown   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Euler_Godunov;
6379566063dSJacob Faibussowitsch   PetscCall(PetscNew(&eu));
638c4762a1bSJed Brown   phys->data     = eu;
639c4762a1bSJed Brown   mod->setupbc   = SetUpBC_Euler;
640de8de29fSMatthew G. Knepley   mod->setupCEED = SetupCEED_Euler;
641de8de29fSMatthew G. Knepley 
642de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov", PhysicsRiemann_Euler_Godunov));
643de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
644de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov_ceed", PhysicsRiemann_Euler_Godunov_CEED));
645de8de29fSMatthew G. Knepley #endif
646de8de29fSMatthew G. Knepley 
647d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
648c4762a1bSJed Brown   {
649de8de29fSMatthew G. Knepley     void (*PhysicsRiemann_Euler)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
650c4762a1bSJed Brown     PetscReal alpha;
651c4762a1bSJed Brown     char      type[64]       = "linear_wave";
652de8de29fSMatthew G. Knepley     char      eu_riemann[64] = "godunov";
653c4762a1bSJed Brown     PetscBool is;
654de8de29fSMatthew G. Knepley     eu->gamma = 1.4;
655de8de29fSMatthew G. Knepley     eu->amach = 2.02;
656de8de29fSMatthew G. Knepley     eu->rhoR  = 3.0;
657de8de29fSMatthew G. Knepley     eu->itana = 0.57735026918963; /* angle of Euler self similar (SS) shock */
658de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsFList("-eu_riemann", "Riemann solver", "", PhysicsRiemannList_Euler, eu_riemann, eu_riemann, sizeof eu_riemann, NULL));
659de8de29fSMatthew G. Knepley     PetscCall(PetscFunctionListFind(PhysicsRiemannList_Euler, eu_riemann, &PhysicsRiemann_Euler));
660de8de29fSMatthew G. Knepley     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Euler;
661de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->gamma, &eu->gamma, NULL));
662de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->amach, &eu->amach, NULL));
663de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->rhoR, &eu->rhoR, NULL));
664c4762a1bSJed Brown     alpha = 60.;
66552ce0ab5SPierre Jolivet     PetscCall(PetscOptionsRangeReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL, 0.0, 90.0));
666de8de29fSMatthew G. Knepley     eu->itana = 1. / PetscTanReal(alpha * PETSC_PI / 180.0);
6679566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL));
6689566063dSJacob Faibussowitsch     PetscCall(PetscStrcmp(type, "linear_wave", &is));
669c4762a1bSJed Brown     if (is) {
67030602db0SMatthew G. Knepley       /* Remember this should be periodic */
671c4762a1bSJed Brown       eu->type = EULER_LINEAR_WAVE;
6729566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave"));
6739371c9d4SSatish Balay     } else {
67454c59aa7SJacob Faibussowitsch       PetscCheck(DIM == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "DIM must be 2 unless linear wave test %s", type);
6759566063dSJacob Faibussowitsch       PetscCall(PetscStrcmp(type, "iv_shock", &is));
676c4762a1bSJed Brown       if (is) {
677c4762a1bSJed Brown         eu->type = EULER_IV_SHOCK;
6789566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock"));
6799371c9d4SSatish Balay       } else {
6809566063dSJacob Faibussowitsch         PetscCall(PetscStrcmp(type, "ss_shock", &is));
681c4762a1bSJed Brown         if (is) {
682c4762a1bSJed Brown           eu->type = EULER_SS_SHOCK;
6839566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock"));
6849371c9d4SSatish Balay         } else {
6859566063dSJacob Faibussowitsch           PetscCall(PetscStrcmp(type, "shock_tube", &is));
686c4762a1bSJed Brown           if (is) eu->type = EULER_SHOCK_TUBE;
68798921bdaSJacob Faibussowitsch           else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type);
6889566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube"));
689c4762a1bSJed Brown         }
690c4762a1bSJed Brown       }
691c4762a1bSJed Brown     }
692c4762a1bSJed Brown   }
693d0609cedSBarry Smith   PetscOptionsHeadEnd();
694c4762a1bSJed Brown   phys->maxspeed = 0.; /* will get set in solution */
6959566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
6969566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
6979566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
6989566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
6999566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
7009566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
7013ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
702c4762a1bSJed Brown }
703c4762a1bSJed Brown 
704d71ae5a4SJacob Faibussowitsch static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx)
705d71ae5a4SJacob Faibussowitsch {
706c4762a1bSJed Brown   PetscReal err = 0.;
707c4762a1bSJed Brown   PetscInt  i, j;
708c4762a1bSJed Brown 
709c4762a1bSJed Brown   PetscFunctionBeginUser;
710c4762a1bSJed Brown   for (i = 0; i < numComps; i++) {
711ad540459SPierre Jolivet     for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
712c4762a1bSJed Brown   }
713c4762a1bSJed Brown   *error = volume * err;
7143ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
715c4762a1bSJed Brown }
716c4762a1bSJed Brown 
717d71ae5a4SJacob Faibussowitsch PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
718d71ae5a4SJacob Faibussowitsch {
719c4762a1bSJed Brown   PetscSF      sfPoint;
720c4762a1bSJed Brown   PetscSection coordSection;
721c4762a1bSJed Brown   Vec          coordinates;
722c4762a1bSJed Brown   PetscSection sectionCell;
723c4762a1bSJed Brown   PetscScalar *part;
724c4762a1bSJed Brown   PetscInt     cStart, cEnd, c;
725c4762a1bSJed Brown   PetscMPIInt  rank;
726c4762a1bSJed Brown 
727c4762a1bSJed Brown   PetscFunctionBeginUser;
7289566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7299566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7309566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, dmCell));
7319566063dSJacob Faibussowitsch   PetscCall(DMGetPointSF(dm, &sfPoint));
7329566063dSJacob Faibussowitsch   PetscCall(DMSetPointSF(*dmCell, sfPoint));
7339566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
7349566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
7359566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
7369566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
7379566063dSJacob Faibussowitsch   PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
7389566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
73948a46eb9SPierre Jolivet   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
7409566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionCell));
7419566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(*dmCell, sectionCell));
7429566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionCell));
7439566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(*dmCell, partition));
7449566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
7459566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*partition, &part));
746c4762a1bSJed Brown   for (c = cStart; c < cEnd; ++c) {
747c4762a1bSJed Brown     PetscScalar *p;
748c4762a1bSJed Brown 
7499566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
750c4762a1bSJed Brown     p[0] = rank;
751c4762a1bSJed Brown   }
7529566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*partition, &part));
7533ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
754c4762a1bSJed Brown }
755c4762a1bSJed Brown 
756d71ae5a4SJacob Faibussowitsch PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
757d71ae5a4SJacob Faibussowitsch {
7583e9753d6SMatthew G. Knepley   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
759c4762a1bSJed Brown   PetscSection       coordSection;
760c4762a1bSJed Brown   Vec                coordinates, facegeom, cellgeom;
761c4762a1bSJed Brown   PetscSection       sectionMass;
762c4762a1bSJed Brown   PetscScalar       *m;
763c4762a1bSJed Brown   const PetscScalar *fgeom, *cgeom, *coords;
764c4762a1bSJed Brown   PetscInt           vStart, vEnd, v;
765c4762a1bSJed Brown 
766c4762a1bSJed Brown   PetscFunctionBeginUser;
7679566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
7689566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7699566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7709566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, &dmMass));
7719566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
7729566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
7739566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
7749566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
7759566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
776c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
777c4762a1bSJed Brown     PetscInt numFaces;
778c4762a1bSJed Brown 
7799566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
7809566063dSJacob Faibussowitsch     PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
781c4762a1bSJed Brown   }
7829566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionMass));
7839566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(dmMass, sectionMass));
7849566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionMass));
7859566063dSJacob Faibussowitsch   PetscCall(DMGetLocalVector(dmMass, massMatrix));
7869566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*massMatrix, &m));
7879566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
7889566063dSJacob Faibussowitsch   PetscCall(VecGetDM(facegeom, &dmFace));
7899566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(facegeom, &fgeom));
7909566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellgeom, &dmCell));
7919566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
7929566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
7939566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(coordinates, &coords));
794c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
795c4762a1bSJed Brown     const PetscInt  *faces;
796c4762a1bSJed Brown     PetscFVFaceGeom *fgA, *fgB, *cg;
797c4762a1bSJed Brown     PetscScalar     *vertex;
798c4762a1bSJed Brown     PetscInt         numFaces, sides[2], f, g;
799c4762a1bSJed Brown 
8009566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
8019566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
8029566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupport(dmMass, v, &faces));
803c4762a1bSJed Brown     for (f = 0; f < numFaces; ++f) {
804c4762a1bSJed Brown       sides[0] = faces[f];
8059566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
806c4762a1bSJed Brown       for (g = 0; g < numFaces; ++g) {
807c4762a1bSJed Brown         const PetscInt *cells = NULL;
808c4762a1bSJed Brown         PetscReal       area  = 0.0;
809c4762a1bSJed Brown         PetscInt        numCells;
810c4762a1bSJed Brown 
811c4762a1bSJed Brown         sides[1] = faces[g];
8129566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
8139566063dSJacob Faibussowitsch         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
81454c59aa7SJacob Faibussowitsch         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
8159566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
816c4762a1bSJed Brown         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
817c4762a1bSJed Brown         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
818c4762a1bSJed Brown         m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
8199566063dSJacob Faibussowitsch         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
820c4762a1bSJed Brown       }
821c4762a1bSJed Brown     }
822c4762a1bSJed Brown   }
8239566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
8249566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
8259566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(coordinates, &coords));
8269566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*massMatrix, &m));
8279566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmMass));
8289566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
8293ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
830c4762a1bSJed Brown }
831c4762a1bSJed Brown 
832c4762a1bSJed Brown /* Behavior will be different for multi-physics or when using non-default boundary conditions */
833d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
834d71ae5a4SJacob Faibussowitsch {
835c4762a1bSJed Brown   PetscFunctionBeginUser;
836c4762a1bSJed Brown   mod->solution    = func;
837c4762a1bSJed Brown   mod->solutionctx = ctx;
8383ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
839c4762a1bSJed Brown }
840c4762a1bSJed Brown 
841d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
842d71ae5a4SJacob Faibussowitsch {
843c4762a1bSJed Brown   FunctionalLink link, *ptr;
844c4762a1bSJed Brown   PetscInt       lastoffset = -1;
845c4762a1bSJed Brown 
846c4762a1bSJed Brown   PetscFunctionBeginUser;
847c4762a1bSJed Brown   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
8489566063dSJacob Faibussowitsch   PetscCall(PetscNew(&link));
8499566063dSJacob Faibussowitsch   PetscCall(PetscStrallocpy(name, &link->name));
850c4762a1bSJed Brown   link->offset = lastoffset + 1;
851c4762a1bSJed Brown   link->func   = func;
852c4762a1bSJed Brown   link->ctx    = ctx;
853c4762a1bSJed Brown   link->next   = NULL;
854c4762a1bSJed Brown   *ptr         = link;
855c4762a1bSJed Brown   *offset      = link->offset;
8563ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
857c4762a1bSJed Brown }
858c4762a1bSJed Brown 
859d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject)
860d71ae5a4SJacob Faibussowitsch {
861c4762a1bSJed Brown   PetscInt       i, j;
862c4762a1bSJed Brown   FunctionalLink link;
863c4762a1bSJed Brown   char          *names[256];
864c4762a1bSJed Brown 
865c4762a1bSJed Brown   PetscFunctionBeginUser;
866dd39110bSPierre Jolivet   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
8679566063dSJacob Faibussowitsch   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
868c4762a1bSJed Brown   /* Create list of functionals that will be computed somehow */
8699566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
870c4762a1bSJed Brown   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
8719566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
872c4762a1bSJed Brown   mod->numCall = 0;
873c4762a1bSJed Brown   for (i = 0; i < mod->numMonitored; i++) {
874c4762a1bSJed Brown     for (link = mod->functionalRegistry; link; link = link->next) {
875c4762a1bSJed Brown       PetscBool match;
8769566063dSJacob Faibussowitsch       PetscCall(PetscStrcasecmp(names[i], link->name, &match));
877c4762a1bSJed Brown       if (match) break;
878c4762a1bSJed Brown     }
87954c59aa7SJacob Faibussowitsch     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
880c4762a1bSJed Brown     mod->functionalMonitored[i] = link;
881c4762a1bSJed Brown     for (j = 0; j < i; j++) {
882c4762a1bSJed Brown       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
883c4762a1bSJed Brown     }
884c4762a1bSJed Brown     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
885c4762a1bSJed Brown   next_name:
8869566063dSJacob Faibussowitsch     PetscCall(PetscFree(names[i]));
887c4762a1bSJed Brown   }
888c4762a1bSJed Brown 
889c4762a1bSJed Brown   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
890c4762a1bSJed Brown   mod->maxComputed = -1;
891c4762a1bSJed Brown   for (link = mod->functionalRegistry; link; link = link->next) {
892c4762a1bSJed Brown     for (i = 0; i < mod->numCall; i++) {
893c4762a1bSJed Brown       FunctionalLink call = mod->functionalCall[i];
894ad540459SPierre Jolivet       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
895c4762a1bSJed Brown     }
896c4762a1bSJed Brown   }
8973ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
898c4762a1bSJed Brown }
899c4762a1bSJed Brown 
900d71ae5a4SJacob Faibussowitsch static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
901d71ae5a4SJacob Faibussowitsch {
902c4762a1bSJed Brown   FunctionalLink l, next;
903c4762a1bSJed Brown 
904c4762a1bSJed Brown   PetscFunctionBeginUser;
9053ba16761SJacob Faibussowitsch   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
906c4762a1bSJed Brown   l     = *link;
907c4762a1bSJed Brown   *link = NULL;
908c4762a1bSJed Brown   for (; l; l = next) {
909c4762a1bSJed Brown     next = l->next;
9109566063dSJacob Faibussowitsch     PetscCall(PetscFree(l->name));
9119566063dSJacob Faibussowitsch     PetscCall(PetscFree(l));
912c4762a1bSJed Brown   }
9133ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
914c4762a1bSJed Brown }
915c4762a1bSJed Brown 
916c4762a1bSJed Brown /* put the solution callback into a functional callback */
917d71ae5a4SJacob Faibussowitsch static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
918d71ae5a4SJacob Faibussowitsch {
919c4762a1bSJed Brown   Model mod;
920*4d86920dSPierre Jolivet 
9217510d9b0SBarry Smith   PetscFunctionBeginUser;
922c4762a1bSJed Brown   mod = (Model)modctx;
9239566063dSJacob Faibussowitsch   PetscCall((*mod->solution)(mod, time, x, u, mod->solutionctx));
9243ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
925c4762a1bSJed Brown }
926c4762a1bSJed Brown 
927d71ae5a4SJacob Faibussowitsch PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
928d71ae5a4SJacob Faibussowitsch {
929c4762a1bSJed Brown   PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
930c4762a1bSJed Brown   void *ctx[1];
931c4762a1bSJed Brown   Model mod = user->model;
932c4762a1bSJed Brown 
933c4762a1bSJed Brown   PetscFunctionBeginUser;
934c4762a1bSJed Brown   func[0] = SolutionFunctional;
935c4762a1bSJed Brown   ctx[0]  = (void *)mod;
9369566063dSJacob Faibussowitsch   PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X));
9373ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
938c4762a1bSJed Brown }
939c4762a1bSJed Brown 
940d71ae5a4SJacob Faibussowitsch static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
941d71ae5a4SJacob Faibussowitsch {
942c4762a1bSJed Brown   PetscFunctionBeginUser;
9439566063dSJacob Faibussowitsch   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
9449566063dSJacob Faibussowitsch   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
9459566063dSJacob Faibussowitsch   PetscCall(PetscViewerFileSetName(*viewer, filename));
9463ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
947c4762a1bSJed Brown }
948c4762a1bSJed Brown 
949d71ae5a4SJacob Faibussowitsch static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
950d71ae5a4SJacob Faibussowitsch {
951c4762a1bSJed Brown   User        user = (User)ctx;
9523e9753d6SMatthew G. Knepley   DM          dm, plex;
953c4762a1bSJed Brown   PetscViewer viewer;
954c4762a1bSJed Brown   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
955c4762a1bSJed Brown   PetscReal   xnorm;
956c4762a1bSJed Brown 
957c4762a1bSJed Brown   PetscFunctionBeginUser;
9589566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "u"));
9599566063dSJacob Faibussowitsch   PetscCall(VecGetDM(X, &dm));
9609566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));
961c4762a1bSJed Brown 
962ad540459SPierre Jolivet   if (stepnum >= 0) stepnum += user->monitorStepOffset;
963c4762a1bSJed Brown   if (stepnum >= 0) { /* No summary for final time */
964c4762a1bSJed Brown     Model              mod = user->model;
9653e9753d6SMatthew G. Knepley     Vec                cellgeom;
966c4762a1bSJed Brown     PetscInt           c, cStart, cEnd, fcount, i;
967c4762a1bSJed Brown     size_t             ftableused, ftablealloc;
968c4762a1bSJed Brown     const PetscScalar *cgeom, *x;
969c4762a1bSJed Brown     DM                 dmCell;
970c4762a1bSJed Brown     DMLabel            vtkLabel;
971c4762a1bSJed Brown     PetscReal         *fmin, *fmax, *fintegral, *ftmp;
9723e9753d6SMatthew G. Knepley 
9739566063dSJacob Faibussowitsch     PetscCall(DMConvert(dm, DMPLEX, &plex));
9749566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
975c4762a1bSJed Brown     fcount = mod->maxComputed + 1;
9769566063dSJacob Faibussowitsch     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
977c4762a1bSJed Brown     for (i = 0; i < fcount; i++) {
978c4762a1bSJed Brown       fmin[i]      = PETSC_MAX_REAL;
979c4762a1bSJed Brown       fmax[i]      = PETSC_MIN_REAL;
980c4762a1bSJed Brown       fintegral[i] = 0;
981c4762a1bSJed Brown     }
9829566063dSJacob Faibussowitsch     PetscCall(VecGetDM(cellgeom, &dmCell));
9839566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd));
9849566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
9859566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(X, &x));
9869566063dSJacob Faibussowitsch     PetscCall(DMGetLabel(dm, "vtk", &vtkLabel));
987c4762a1bSJed Brown     for (c = cStart; c < cEnd; ++c) {
988c4762a1bSJed Brown       PetscFVCellGeom   *cg;
989c4762a1bSJed Brown       const PetscScalar *cx     = NULL;
990c4762a1bSJed Brown       PetscInt           vtkVal = 0;
991c4762a1bSJed Brown 
992c4762a1bSJed Brown       /* not that these two routines as currently implemented work for any dm with a
993c4762a1bSJed Brown        * localSection/globalSection */
9949566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
9959566063dSJacob Faibussowitsch       PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
9969566063dSJacob Faibussowitsch       if (vtkLabel) PetscCall(DMLabelGetValue(vtkLabel, c, &vtkVal));
997c4762a1bSJed Brown       if (!vtkVal || !cx) continue; /* ghost, or not a global cell */
998c4762a1bSJed Brown       for (i = 0; i < mod->numCall; i++) {
999c4762a1bSJed Brown         FunctionalLink flink = mod->functionalCall[i];
10009566063dSJacob Faibussowitsch         PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
1001c4762a1bSJed Brown       }
1002c4762a1bSJed Brown       for (i = 0; i < fcount; i++) {
1003c4762a1bSJed Brown         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1004c4762a1bSJed Brown         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1005c4762a1bSJed Brown         fintegral[i] += cg->volume * ftmp[i];
1006c4762a1bSJed Brown       }
1007c4762a1bSJed Brown     }
10089566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
10099566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(X, &x));
10109566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&plex));
1011712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1012712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1013712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
1014c4762a1bSJed Brown 
1015c4762a1bSJed Brown     ftablealloc = fcount * 100;
1016c4762a1bSJed Brown     ftableused  = 0;
10179566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(ftablealloc, &ftable));
1018c4762a1bSJed Brown     for (i = 0; i < mod->numMonitored; i++) {
1019c4762a1bSJed Brown       size_t         countused;
1020c4762a1bSJed Brown       char           buffer[256], *p;
1021c4762a1bSJed Brown       FunctionalLink flink = mod->functionalMonitored[i];
1022c4762a1bSJed Brown       PetscInt       id    = flink->offset;
1023c4762a1bSJed Brown       if (i % 3) {
10249566063dSJacob Faibussowitsch         PetscCall(PetscArraycpy(buffer, "  ", 2));
1025c4762a1bSJed Brown         p = buffer + 2;
1026c4762a1bSJed Brown       } else if (i) {
1027c4762a1bSJed Brown         char newline[] = "\n";
10289566063dSJacob Faibussowitsch         PetscCall(PetscMemcpy(buffer, newline, sizeof(newline) - 1));
1029c4762a1bSJed Brown         p = buffer + sizeof(newline) - 1;
1030c4762a1bSJed Brown       } else {
1031c4762a1bSJed Brown         p = buffer;
1032c4762a1bSJed Brown       }
10339566063dSJacob Faibussowitsch       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]));
1034c4762a1bSJed Brown       countused--;
1035c4762a1bSJed Brown       countused += p - buffer;
1036c4762a1bSJed Brown       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1037c4762a1bSJed Brown         char *ftablenew;
1038c4762a1bSJed Brown         ftablealloc = 2 * ftablealloc + countused;
10399566063dSJacob Faibussowitsch         PetscCall(PetscMalloc(ftablealloc, &ftablenew));
10409566063dSJacob Faibussowitsch         PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
10419566063dSJacob Faibussowitsch         PetscCall(PetscFree(ftable));
1042c4762a1bSJed Brown         ftable = ftablenew;
1043c4762a1bSJed Brown       }
10449566063dSJacob Faibussowitsch       PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
1045c4762a1bSJed Brown       ftableused += countused;
1046c4762a1bSJed Brown       ftable[ftableused] = 0;
1047c4762a1bSJed Brown     }
10489566063dSJacob Faibussowitsch     PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));
1049c4762a1bSJed Brown 
105063a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
10519566063dSJacob Faibussowitsch     PetscCall(PetscFree(ftable));
1052c4762a1bSJed Brown   }
10533ba16761SJacob Faibussowitsch   if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
1054c4762a1bSJed Brown   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1055c4762a1bSJed Brown     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
10569566063dSJacob Faibussowitsch       PetscCall(TSGetStepNumber(ts, &stepnum));
1057c4762a1bSJed Brown     }
105863a3b9bcSJacob Faibussowitsch     PetscCall(PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum));
10599566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, filename, &viewer));
10609566063dSJacob Faibussowitsch     PetscCall(VecView(X, viewer));
10619566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
1062c4762a1bSJed Brown   }
10633ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1064c4762a1bSJed Brown }
1065c4762a1bSJed Brown 
1066d71ae5a4SJacob Faibussowitsch static PetscErrorCode initializeTS(DM dm, User user, TS *ts)
1067d71ae5a4SJacob Faibussowitsch {
1068de8de29fSMatthew G. Knepley   PetscBool useCeed;
1069de8de29fSMatthew G. Knepley 
10707510d9b0SBarry Smith   PetscFunctionBeginUser;
10719566063dSJacob Faibussowitsch   PetscCall(TSCreate(PetscObjectComm((PetscObject)dm), ts));
10729566063dSJacob Faibussowitsch   PetscCall(TSSetType(*ts, TSSSP));
10739566063dSJacob Faibussowitsch   PetscCall(TSSetDM(*ts, dm));
10741baa6e33SBarry Smith   if (user->vtkmon) PetscCall(TSMonitorSet(*ts, MonitorVTK, user, NULL));
1075de8de29fSMatthew G. Knepley   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
10769566063dSJacob Faibussowitsch   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user));
1077de8de29fSMatthew G. Knepley   if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user));
1078de8de29fSMatthew G. Knepley   else PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
10799566063dSJacob Faibussowitsch   PetscCall(TSSetMaxTime(*ts, 2.0));
10809566063dSJacob Faibussowitsch   PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
10813ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1082c4762a1bSJed Brown }
1083c4762a1bSJed Brown 
108452115386SStefano Zampini typedef struct {
108552115386SStefano Zampini   PetscFV      fvm;
108652115386SStefano Zampini   VecTagger    refineTag;
108752115386SStefano Zampini   VecTagger    coarsenTag;
108852115386SStefano Zampini   DM           adaptedDM;
108952115386SStefano Zampini   User         user;
109052115386SStefano Zampini   PetscReal    cfl;
109152115386SStefano Zampini   PetscLimiter limiter;
109252115386SStefano Zampini   PetscLimiter noneLimiter;
109352115386SStefano Zampini } TransferCtx;
109452115386SStefano Zampini 
109552115386SStefano Zampini static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
1096d71ae5a4SJacob Faibussowitsch {
109752115386SStefano Zampini   TransferCtx       *tctx       = (TransferCtx *)ctx;
109852115386SStefano Zampini   PetscFV            fvm        = tctx->fvm;
109952115386SStefano Zampini   VecTagger          refineTag  = tctx->refineTag;
110052115386SStefano Zampini   VecTagger          coarsenTag = tctx->coarsenTag;
110152115386SStefano Zampini   User               user       = tctx->user;
1102c4762a1bSJed Brown   DM                 dm, gradDM, plex, cellDM, adaptedDM = NULL;
1103c4762a1bSJed Brown   Vec                cellGeom, faceGeom;
1104c4762a1bSJed Brown   PetscBool          isForest, computeGradient;
1105c4762a1bSJed Brown   Vec                grad, locGrad, locX, errVec;
1106c4762a1bSJed Brown   PetscInt           cStart, cEnd, c, dim, nRefine, nCoarsen;
110752115386SStefano Zampini   PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1108c4762a1bSJed Brown   PetscScalar       *errArray;
1109c4762a1bSJed Brown   const PetscScalar *pointVals;
1110c4762a1bSJed Brown   const PetscScalar *pointGrads;
1111c4762a1bSJed Brown   const PetscScalar *pointGeom;
1112c4762a1bSJed Brown   DMLabel            adaptLabel = NULL;
1113c4762a1bSJed Brown   IS                 refineIS, coarsenIS;
1114c4762a1bSJed Brown 
11157510d9b0SBarry Smith   PetscFunctionBeginUser;
111652115386SStefano Zampini   *resize = PETSC_FALSE;
11179566063dSJacob Faibussowitsch   PetscCall(VecGetDM(sol, &dm));
11189566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
111952115386SStefano Zampini   PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
11209566063dSJacob Faibussowitsch   PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
11219566063dSJacob Faibussowitsch   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
11229566063dSJacob Faibussowitsch   PetscCall(DMIsForest(dm, &isForest));
11239566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
11249566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
11259566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(plex, &locX));
11269566063dSJacob Faibussowitsch   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
11279566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
11289566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
11299566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(gradDM, &grad));
11309566063dSJacob Faibussowitsch   PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
11319566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(gradDM, &locGrad));
11329566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
11339566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
11349566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&grad));
11359566063dSJacob Faibussowitsch   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
11369566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(locGrad, &pointGrads));
11379566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
11389566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(locX, &pointVals));
11399566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellGeom, &cellDM));
11409566063dSJacob Faibussowitsch   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
114177433607SBarry Smith   PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
11429566063dSJacob Faibussowitsch   PetscCall(VecSetUp(errVec));
11439566063dSJacob Faibussowitsch   PetscCall(VecGetArray(errVec, &errArray));
1144c4762a1bSJed Brown   for (c = cStart; c < cEnd; c++) {
1145c4762a1bSJed Brown     PetscReal        errInd = 0.;
1146c4762a1bSJed Brown     PetscScalar     *pointGrad;
1147c4762a1bSJed Brown     PetscScalar     *pointVal;
1148c4762a1bSJed Brown     PetscFVCellGeom *cg;
1149c4762a1bSJed Brown 
11509566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
11519566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
11529566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));
1153c4762a1bSJed Brown 
11549566063dSJacob Faibussowitsch     PetscCall((user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1155c4762a1bSJed Brown     errArray[c - cStart] = errInd;
1156c4762a1bSJed Brown     minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
1157c4762a1bSJed Brown     minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
1158c4762a1bSJed Brown   }
11599566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(errVec, &errArray));
11609566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(locX, &pointVals));
11619566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
11629566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
11639566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&locGrad));
11649566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&locX));
11659566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1166c4762a1bSJed Brown 
11679566063dSJacob Faibussowitsch   PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
11689566063dSJacob Faibussowitsch   PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
11699566063dSJacob Faibussowitsch   PetscCall(ISGetSize(refineIS, &nRefine));
11709566063dSJacob Faibussowitsch   PetscCall(ISGetSize(coarsenIS, &nCoarsen));
11719566063dSJacob Faibussowitsch   if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
11729566063dSJacob Faibussowitsch   if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
11739566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&coarsenIS));
11749566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&refineIS));
11759566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&errVec));
1176c4762a1bSJed Brown 
11779566063dSJacob Faibussowitsch   PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
117852115386SStefano Zampini   PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1179c4762a1bSJed Brown   minMaxInd[1] = -minMaxInd[1];
1180712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
118152115386SStefano Zampini   PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1182c4762a1bSJed Brown   if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
11839566063dSJacob Faibussowitsch     PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1184c4762a1bSJed Brown   }
11859566063dSJacob Faibussowitsch   PetscCall(DMLabelDestroy(&adaptLabel));
1186c4762a1bSJed Brown   if (adaptedDM) {
118752115386SStefano Zampini     tctx->adaptedDM = adaptedDM;
118852115386SStefano Zampini     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nstep, nRefine, nCoarsen));
118952115386SStefano Zampini     *resize = PETSC_TRUE;
1190c4762a1bSJed Brown   } else {
119152115386SStefano Zampini     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1192c4762a1bSJed Brown   }
11933ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1194c4762a1bSJed Brown }
1195c4762a1bSJed Brown 
119652115386SStefano Zampini static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
119752115386SStefano Zampini {
119852115386SStefano Zampini   TransferCtx *tctx = (TransferCtx *)ctx;
119952115386SStefano Zampini   DM           dm;
120052115386SStefano Zampini   PetscReal    time;
120152115386SStefano Zampini 
120252115386SStefano Zampini   PetscFunctionBeginUser;
120352115386SStefano Zampini   PetscCall(TSGetDM(ts, &dm));
120452115386SStefano Zampini   PetscCall(TSGetTime(ts, &time));
120552115386SStefano Zampini   PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
120652115386SStefano Zampini   for (PetscInt i = 0; i < nv; i++) {
120752115386SStefano Zampini     const char *name;
120852115386SStefano Zampini 
120952115386SStefano Zampini     PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
121052115386SStefano Zampini     PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
121152115386SStefano Zampini     PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
121252115386SStefano Zampini     PetscCall(PetscObjectSetName((PetscObject)vecsout[i], name));
121352115386SStefano Zampini   }
121452115386SStefano Zampini   PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */
121552115386SStefano Zampini 
121652115386SStefano Zampini   Model     mod  = tctx->user->model;
121752115386SStefano Zampini   Physics   phys = mod->physics;
121852115386SStefano Zampini   PetscReal minRadius;
121952115386SStefano Zampini 
122052115386SStefano Zampini   PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius));
122152115386SStefano Zampini   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
122252115386SStefano Zampini   PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed");
122352115386SStefano Zampini 
122452115386SStefano Zampini   PetscReal dt = tctx->cfl * minRadius / mod->maxspeed;
122552115386SStefano Zampini   PetscCall(TSSetTimeStep(ts, dt));
122652115386SStefano Zampini 
122752115386SStefano Zampini   PetscCall(TSSetDM(ts, tctx->adaptedDM));
122852115386SStefano Zampini   PetscCall(DMDestroy(&tctx->adaptedDM));
122952115386SStefano Zampini   PetscFunctionReturn(PETSC_SUCCESS);
123052115386SStefano Zampini }
123152115386SStefano Zampini 
1232d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
1233d71ae5a4SJacob Faibussowitsch {
1234c4762a1bSJed Brown   MPI_Comm          comm;
1235c4762a1bSJed Brown   PetscDS           prob;
1236c4762a1bSJed Brown   PetscFV           fvm;
1237c4762a1bSJed Brown   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1238c4762a1bSJed Brown   User              user;
1239c4762a1bSJed Brown   Model             mod;
1240c4762a1bSJed Brown   Physics           phys;
12413e9753d6SMatthew G. Knepley   DM                dm, plex;
1242c4762a1bSJed Brown   PetscReal         ftime, cfl, dt, minRadius;
1243c4762a1bSJed Brown   PetscInt          dim, nsteps;
1244c4762a1bSJed Brown   TS                ts;
1245c4762a1bSJed Brown   TSConvergedReason reason;
1246c4762a1bSJed Brown   Vec               X;
1247c4762a1bSJed Brown   PetscViewer       viewer;
1248b5a892a1SMatthew G. Knepley   PetscBool         vtkCellGeom, useAMR;
124930602db0SMatthew G. Knepley   PetscInt          adaptInterval;
1250c4762a1bSJed Brown   char              physname[256] = "advect";
1251c4762a1bSJed Brown   VecTagger         refineTag = NULL, coarsenTag = NULL;
125252115386SStefano Zampini   TransferCtx       tctx;
1253c4762a1bSJed Brown 
1254327415f7SBarry Smith   PetscFunctionBeginUser;
12559566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1256c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
1257c4762a1bSJed Brown 
12589566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user));
12599566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model));
12609566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model->physics));
1261c4762a1bSJed Brown   mod           = user->model;
1262c4762a1bSJed Brown   phys          = mod->physics;
1263c4762a1bSJed Brown   mod->comm     = comm;
1264c4762a1bSJed Brown   useAMR        = PETSC_FALSE;
1265c4762a1bSJed Brown   adaptInterval = 1;
1266c4762a1bSJed Brown 
1267c4762a1bSJed Brown   /* Register physical models to be available on the command line */
12689566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
12699566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
12709566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1271c4762a1bSJed Brown 
1272d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1273c4762a1bSJed Brown   {
1274c4762a1bSJed Brown     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
12759566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1276c4762a1bSJed Brown     user->vtkInterval = 1;
12779566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1278c4762a1bSJed Brown     user->vtkmon = PETSC_TRUE;
12799566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1280c4762a1bSJed Brown     vtkCellGeom = PETSC_FALSE;
1281c6a7a370SJeremy L Thompson     PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
12829566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
12839566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
12849566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
12859566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1286c4762a1bSJed Brown   }
1287d0609cedSBarry Smith   PetscOptionsEnd();
1288c4762a1bSJed Brown 
1289c4762a1bSJed Brown   if (useAMR) {
1290c4762a1bSJed Brown     VecTaggerBox refineBox, coarsenBox;
1291c4762a1bSJed Brown 
1292c4762a1bSJed Brown     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1293c4762a1bSJed Brown     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;
1294c4762a1bSJed Brown 
12959566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &refineTag));
12969566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
12979566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
12989566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
12999566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(refineTag));
13009566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(refineTag));
13019566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
1302c4762a1bSJed Brown 
13039566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &coarsenTag));
13049566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
13059566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
13069566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
13079566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(coarsenTag));
13089566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(coarsenTag));
13099566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1310c4762a1bSJed Brown   }
1311c4762a1bSJed Brown 
1312d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1313c4762a1bSJed Brown   {
1314c4762a1bSJed Brown     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
13159566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
13169566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
13179566063dSJacob Faibussowitsch     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
13189566063dSJacob Faibussowitsch     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1319c4762a1bSJed Brown     /* Count number of fields and dofs */
1320c4762a1bSJed Brown     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
132154c59aa7SJacob Faibussowitsch     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
13229566063dSJacob Faibussowitsch     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1323c4762a1bSJed Brown   }
1324d0609cedSBarry Smith   PetscOptionsEnd();
1325c4762a1bSJed Brown 
1326c4762a1bSJed Brown   /* Create mesh */
1327c4762a1bSJed Brown   {
132830602db0SMatthew G. Knepley     PetscInt i;
132930602db0SMatthew G. Knepley 
13309566063dSJacob Faibussowitsch     PetscCall(DMCreate(comm, &dm));
13319566063dSJacob Faibussowitsch     PetscCall(DMSetType(dm, DMPLEX));
13329566063dSJacob Faibussowitsch     PetscCall(DMSetFromOptions(dm));
13339371c9d4SSatish Balay     for (i = 0; i < DIM; i++) {
13349371c9d4SSatish Balay       mod->bounds[2 * i]     = 0.;
13359371c9d4SSatish Balay       mod->bounds[2 * i + 1] = 1.;
13369371c9d4SSatish Balay     };
1337c4762a1bSJed Brown     dim = DIM;
133830602db0SMatthew G. Knepley     { /* a null name means just do a hex box */
133930602db0SMatthew G. Knepley       PetscInt  cells[3] = {1, 1, 1}, n = 3;
134030602db0SMatthew G. Knepley       PetscBool flg2, skew              = PETSC_FALSE;
1341c4762a1bSJed Brown       PetscInt  nret2 = 2 * DIM;
1342d0609cedSBarry Smith       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
13439566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-grid_bounds", "bounds of the mesh in each direction (i.e., x_min,x_max,y_min,y_max", "", mod->bounds, &nret2, &flg2));
13449566063dSJacob Faibussowitsch       PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
13459566063dSJacob Faibussowitsch       PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1346d0609cedSBarry Smith       PetscOptionsEnd();
134730602db0SMatthew G. Knepley       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1348c4762a1bSJed Brown       if (flg2) {
1349c4762a1bSJed Brown         PetscInt     dimEmbed, i;
1350c4762a1bSJed Brown         PetscInt     nCoords;
1351c4762a1bSJed Brown         PetscScalar *coords;
1352c4762a1bSJed Brown         Vec          coordinates;
1353c4762a1bSJed Brown 
13549566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
13559566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
13569566063dSJacob Faibussowitsch         PetscCall(VecGetLocalSize(coordinates, &nCoords));
135754c59aa7SJacob Faibussowitsch         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
13589566063dSJacob Faibussowitsch         PetscCall(VecGetArray(coordinates, &coords));
1359c4762a1bSJed Brown         for (i = 0; i < nCoords; i += dimEmbed) {
1360c4762a1bSJed Brown           PetscInt j;
1361c4762a1bSJed Brown 
1362c4762a1bSJed Brown           PetscScalar *coord = &coords[i];
1363c4762a1bSJed Brown           for (j = 0; j < dimEmbed; j++) {
1364c4762a1bSJed Brown             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1365c4762a1bSJed Brown             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1366c4762a1bSJed Brown               if (cells[0] == 2 && i == 8) {
1367c4762a1bSJed Brown                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
136830602db0SMatthew G. Knepley               } else if (cells[0] == 3) {
1369c4762a1bSJed Brown                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1370c4762a1bSJed Brown                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1371c4762a1bSJed Brown                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1372c4762a1bSJed Brown               }
1373c4762a1bSJed Brown             }
1374c4762a1bSJed Brown           }
1375c4762a1bSJed Brown         }
13769566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(coordinates, &coords));
13779566063dSJacob Faibussowitsch         PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1378c4762a1bSJed Brown       }
1379c4762a1bSJed Brown     }
1380c4762a1bSJed Brown   }
13819566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
13829566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
1383c4762a1bSJed Brown 
1384c4762a1bSJed Brown   /* set up BCs, functions, tags */
13859566063dSJacob Faibussowitsch   PetscCall(DMCreateLabel(dm, "Face Sets"));
1386c4762a1bSJed Brown   mod->errorIndicator = ErrorIndicator_Simple;
1387c4762a1bSJed Brown 
1388c4762a1bSJed Brown   {
1389c4762a1bSJed Brown     DM gdm;
1390c4762a1bSJed Brown 
13919566063dSJacob Faibussowitsch     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
13929566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dm));
1393c4762a1bSJed Brown     dm = gdm;
13949566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1395c4762a1bSJed Brown   }
1396c4762a1bSJed Brown 
13979566063dSJacob Faibussowitsch   PetscCall(PetscFVCreate(comm, &fvm));
13989566063dSJacob Faibussowitsch   PetscCall(PetscFVSetFromOptions(fvm));
13999566063dSJacob Faibussowitsch   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
14009566063dSJacob Faibussowitsch   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
14019566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1402c4762a1bSJed Brown   {
1403c4762a1bSJed Brown     PetscInt f, dof;
1404c4762a1bSJed Brown     for (f = 0, dof = 0; f < phys->nfields; f++) {
1405c4762a1bSJed Brown       PetscInt newDof = phys->field_desc[f].dof;
1406c4762a1bSJed Brown 
1407c4762a1bSJed Brown       if (newDof == 1) {
14089566063dSJacob Faibussowitsch         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
14099371c9d4SSatish Balay       } else {
1410c4762a1bSJed Brown         PetscInt j;
1411c4762a1bSJed Brown 
1412c4762a1bSJed Brown         for (j = 0; j < newDof; j++) {
1413c4762a1bSJed Brown           char compName[256] = "Unknown";
1414c4762a1bSJed Brown 
141563a3b9bcSJacob Faibussowitsch           PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
14169566063dSJacob Faibussowitsch           PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1417c4762a1bSJed Brown         }
1418c4762a1bSJed Brown       }
1419c4762a1bSJed Brown       dof += newDof;
1420c4762a1bSJed Brown     }
1421c4762a1bSJed Brown   }
1422c4762a1bSJed Brown   /* FV is now structured with one field having all physics as components */
14239566063dSJacob Faibussowitsch   PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
14249566063dSJacob Faibussowitsch   PetscCall(DMCreateDS(dm));
14259566063dSJacob Faibussowitsch   PetscCall(DMGetDS(dm, &prob));
14269566063dSJacob Faibussowitsch   PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
14279566063dSJacob Faibussowitsch   PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
14289566063dSJacob Faibussowitsch   PetscCall((*mod->setupbc)(dm, prob, phys));
14299566063dSJacob Faibussowitsch   PetscCall(PetscDSSetFromOptions(prob));
1430c4762a1bSJed Brown   {
1431c4762a1bSJed Brown     char      convType[256];
1432c4762a1bSJed Brown     PetscBool flg;
1433c4762a1bSJed Brown 
1434d0609cedSBarry Smith     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
14359566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1436d0609cedSBarry Smith     PetscOptionsEnd();
1437c4762a1bSJed Brown     if (flg) {
1438c4762a1bSJed Brown       DM dmConv;
1439c4762a1bSJed Brown 
14409566063dSJacob Faibussowitsch       PetscCall(DMConvert(dm, convType, &dmConv));
1441c4762a1bSJed Brown       if (dmConv) {
14429566063dSJacob Faibussowitsch         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
14439566063dSJacob Faibussowitsch         PetscCall(DMDestroy(&dm));
1444c4762a1bSJed Brown         dm = dmConv;
14459566063dSJacob Faibussowitsch         PetscCall(DMSetFromOptions(dm));
1446c4762a1bSJed Brown       }
1447c4762a1bSJed Brown     }
1448c4762a1bSJed Brown   }
1449de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
1450de8de29fSMatthew G. Knepley   {
1451de8de29fSMatthew G. Knepley     PetscBool useCeed;
1452de8de29fSMatthew G. Knepley     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1453de8de29fSMatthew G. Knepley     if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1454de8de29fSMatthew G. Knepley   }
1455de8de29fSMatthew G. Knepley #endif
1456c4762a1bSJed Brown 
14579566063dSJacob Faibussowitsch   PetscCall(initializeTS(dm, user, &ts));
1458c4762a1bSJed Brown 
14599566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(dm, &X));
14609566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
14619566063dSJacob Faibussowitsch   PetscCall(SetInitialCondition(dm, X, user));
1462c4762a1bSJed Brown   if (useAMR) {
1463c4762a1bSJed Brown     PetscInt adaptIter;
1464c4762a1bSJed Brown 
1465c4762a1bSJed Brown     /* use no limiting when reconstructing gradients for adaptivity */
14669566063dSJacob Faibussowitsch     PetscCall(PetscFVGetLimiter(fvm, &limiter));
14679566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)limiter));
14689566063dSJacob Faibussowitsch     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
14699566063dSJacob Faibussowitsch     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
1470c4762a1bSJed Brown 
147152115386SStefano Zampini     /* Refinement context */
147252115386SStefano Zampini     tctx.fvm         = fvm;
147352115386SStefano Zampini     tctx.refineTag   = refineTag;
147452115386SStefano Zampini     tctx.coarsenTag  = coarsenTag;
147552115386SStefano Zampini     tctx.adaptedDM   = NULL;
147652115386SStefano Zampini     tctx.user        = user;
147752115386SStefano Zampini     tctx.noneLimiter = noneLimiter;
147852115386SStefano Zampini     tctx.limiter     = limiter;
147952115386SStefano Zampini     tctx.cfl         = cfl;
148052115386SStefano Zampini 
148152115386SStefano Zampini     /* Do some initial refinement steps */
1482c4762a1bSJed Brown     for (adaptIter = 0;; ++adaptIter) {
1483c4762a1bSJed Brown       PetscLogDouble bytes;
148452115386SStefano Zampini       PetscBool      resize;
1485c4762a1bSJed Brown 
14869566063dSJacob Faibussowitsch       PetscCall(PetscMemoryGetCurrentUsage(&bytes));
148763a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
14889566063dSJacob Faibussowitsch       PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
14899566063dSJacob Faibussowitsch       PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));
1490c4762a1bSJed Brown 
149152115386SStefano Zampini       PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
149252115386SStefano Zampini       if (!resize) break;
14939566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dm));
14949566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&X));
149552115386SStefano Zampini       dm             = tctx.adaptedDM;
149652115386SStefano Zampini       tctx.adaptedDM = NULL;
149752115386SStefano Zampini       PetscCall(TSSetDM(ts, dm));
14989566063dSJacob Faibussowitsch       PetscCall(DMCreateGlobalVector(dm, &X));
14999566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
15009566063dSJacob Faibussowitsch       PetscCall(SetInitialCondition(dm, X, user));
1501c4762a1bSJed Brown     }
1502c4762a1bSJed Brown   }
1503c4762a1bSJed Brown 
15049566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
1505c4762a1bSJed Brown   if (vtkCellGeom) {
1506c4762a1bSJed Brown     DM  dmCell;
1507c4762a1bSJed Brown     Vec cellgeom, partition;
1508c4762a1bSJed Brown 
15099566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
15109566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
15119566063dSJacob Faibussowitsch     PetscCall(VecView(cellgeom, viewer));
15129566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15139566063dSJacob Faibussowitsch     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
15149566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
15159566063dSJacob Faibussowitsch     PetscCall(VecView(partition, viewer));
15169566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15179566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&partition));
15189566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dmCell));
1519c4762a1bSJed Brown   }
1520c4762a1bSJed Brown   /* collect max maxspeed from all processes -- todo */
15219566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
15229566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1523712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
152454c59aa7SJacob Faibussowitsch   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1525c4762a1bSJed Brown   dt = cfl * minRadius / mod->maxspeed;
15269566063dSJacob Faibussowitsch   PetscCall(TSSetTimeStep(ts, dt));
15279566063dSJacob Faibussowitsch   PetscCall(TSSetFromOptions(ts));
1528c4762a1bSJed Brown 
152952115386SStefano Zampini   /* When using adaptive mesh refinement
153052115386SStefano Zampini      specify callbacks to refine the solution
153152115386SStefano Zampini      and interpolate data from old to new mesh */
153252115386SStefano Zampini   if (useAMR) { PetscCall(TSSetResize(ts, adaptToleranceFVMSetUp, Transfer, &tctx)); }
153352115386SStefano Zampini   PetscCall(TSSetSolution(ts, X));
15349566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&X));
153552115386SStefano Zampini   PetscCall(TSSolve(ts, NULL));
15369566063dSJacob Faibussowitsch   PetscCall(TSGetSolveTime(ts, &ftime));
15379566063dSJacob Faibussowitsch   PetscCall(TSGetStepNumber(ts, &nsteps));
15389566063dSJacob Faibussowitsch   PetscCall(TSGetConvergedReason(ts, &reason));
153963a3b9bcSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
15409566063dSJacob Faibussowitsch   PetscCall(TSDestroy(&ts));
1541c4762a1bSJed Brown 
15429566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&refineTag));
15439566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&coarsenTag));
15449566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsList));
15459566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1546de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
15479566063dSJacob Faibussowitsch   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
15489566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalMonitored));
15499566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalCall));
15509566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics->data));
15519566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics));
15529566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model));
15539566063dSJacob Faibussowitsch   PetscCall(PetscFree(user));
15549566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&limiter));
15559566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&noneLimiter));
15569566063dSJacob Faibussowitsch   PetscCall(PetscFVDestroy(&fvm));
15579566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dm));
15589566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
1559b122ec5aSJacob Faibussowitsch   return 0;
1560c4762a1bSJed Brown }
1561c4762a1bSJed Brown 
1562c4762a1bSJed Brown /* Subroutine to set up the initial conditions for the */
1563c4762a1bSJed Brown /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1564c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1565d71ae5a4SJacob Faibussowitsch int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1566d71ae5a4SJacob Faibussowitsch {
1567c4762a1bSJed Brown   int j, k;
1568c4762a1bSJed Brown   /*      Wc=matmul(lv,Ueq) 3 vars */
1569c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1570c4762a1bSJed Brown     wc[k] = 0.;
1571ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1572c4762a1bSJed Brown   }
1573c4762a1bSJed Brown   return 0;
1574c4762a1bSJed Brown }
1575c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1576d71ae5a4SJacob Faibussowitsch int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1577d71ae5a4SJacob Faibussowitsch {
1578c4762a1bSJed Brown   int k, j;
1579c4762a1bSJed Brown   /*      V=matmul(rv,WC) 3 vars */
1580c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1581c4762a1bSJed Brown     v[k] = 0.;
1582ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1583c4762a1bSJed Brown   }
1584c4762a1bSJed Brown   return 0;
1585c4762a1bSJed Brown }
1586c4762a1bSJed Brown /* ---------------------------------------------------------------------- */
1587d71ae5a4SJacob Faibussowitsch int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1588d71ae5a4SJacob Faibussowitsch {
1589c4762a1bSJed Brown   int       j, k;
1590c4762a1bSJed Brown   PetscReal rho, csnd, p0;
1591c4762a1bSJed Brown   /* PetscScalar u; */
1592c4762a1bSJed Brown 
15939371c9d4SSatish Balay   for (k = 0; k < 3; ++k)
15949371c9d4SSatish Balay     for (j = 0; j < 3; ++j) {
15959371c9d4SSatish Balay       lv[k][j] = 0.;
15969371c9d4SSatish Balay       rv[k][j] = 0.;
15979371c9d4SSatish Balay     }
1598c4762a1bSJed Brown   rho = ueq[0];
1599c4762a1bSJed Brown   /* u = ueq[1]; */
1600c4762a1bSJed Brown   p0       = ueq[2];
1601c4762a1bSJed Brown   csnd     = PetscSqrtReal(gamma * p0 / rho);
1602c4762a1bSJed Brown   lv[0][1] = rho * .5;
1603c4762a1bSJed Brown   lv[0][2] = -.5 / csnd;
1604c4762a1bSJed Brown   lv[1][0] = csnd;
1605c4762a1bSJed Brown   lv[1][2] = -1. / csnd;
1606c4762a1bSJed Brown   lv[2][1] = rho * .5;
1607c4762a1bSJed Brown   lv[2][2] = .5 / csnd;
1608c4762a1bSJed Brown   rv[0][0] = -1. / csnd;
1609c4762a1bSJed Brown   rv[1][0] = 1. / rho;
1610c4762a1bSJed Brown   rv[2][0] = -csnd;
1611c4762a1bSJed Brown   rv[0][1] = 1. / csnd;
1612c4762a1bSJed Brown   rv[0][2] = 1. / csnd;
1613c4762a1bSJed Brown   rv[1][2] = 1. / rho;
1614c4762a1bSJed Brown   rv[2][2] = csnd;
1615c4762a1bSJed Brown   return 0;
1616c4762a1bSJed Brown }
1617c4762a1bSJed Brown 
1618d71ae5a4SJacob Faibussowitsch int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1619d71ae5a4SJacob Faibussowitsch {
1620c4762a1bSJed Brown   PetscReal p0, u0, wcp[3], wc[3];
1621c4762a1bSJed Brown   PetscReal lv[3][3];
1622c4762a1bSJed Brown   PetscReal vp[3];
1623c4762a1bSJed Brown   PetscReal rv[3][3];
1624c4762a1bSJed Brown   PetscReal eps, ueq[3], rho0, twopi;
1625c4762a1bSJed Brown 
1626c4762a1bSJed Brown   /* Function Body */
1627c4762a1bSJed Brown   twopi  = 2. * PETSC_PI;
1628c4762a1bSJed Brown   eps    = 1e-4;    /* perturbation */
1629c4762a1bSJed Brown   rho0   = 1e3;     /* density of water */
1630c4762a1bSJed Brown   p0     = 101325.; /* init pressure of 1 atm (?) */
1631c4762a1bSJed Brown   u0     = 0.;
1632c4762a1bSJed Brown   ueq[0] = rho0;
1633c4762a1bSJed Brown   ueq[1] = u0;
1634c4762a1bSJed Brown   ueq[2] = p0;
1635c4762a1bSJed Brown   /* Project initial state to characteristic variables */
1636c4762a1bSJed Brown   eigenvectors(rv, lv, ueq, gamma);
1637c4762a1bSJed Brown   projecteqstate(wc, ueq, lv);
1638c4762a1bSJed Brown   wcp[0] = wc[0];
1639c4762a1bSJed Brown   wcp[1] = wc[1];
1640c4762a1bSJed Brown   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1641c4762a1bSJed Brown   projecttoprim(vp, wcp, rv);
1642c4762a1bSJed Brown   ux->r     = vp[0];         /* density */
1643c4762a1bSJed Brown   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1644c4762a1bSJed Brown   ux->ru[1] = 0.;
1645c4762a1bSJed Brown #if defined DIM > 2
1646c4762a1bSJed Brown   if (dim > 2) ux->ru[2] = 0.;
1647c4762a1bSJed Brown #endif
1648c4762a1bSJed Brown   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1649c4762a1bSJed Brown   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1650c4762a1bSJed Brown   return 0;
1651c4762a1bSJed Brown }
1652c4762a1bSJed Brown 
1653c4762a1bSJed Brown /*TEST
1654c4762a1bSJed Brown 
165530602db0SMatthew G. Knepley   testset:
165630602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0
165730602db0SMatthew G. Knepley 
165830602db0SMatthew G. Knepley     test:
165930602db0SMatthew G. Knepley       suffix: adv_2d_tri_0
166030602db0SMatthew G. Knepley       requires: triangle
166130602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166230602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -simplex -dm_refine 3 -dm_plex_faces 1,1 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
166330602db0SMatthew G. Knepley 
166430602db0SMatthew G. Knepley     test:
166530602db0SMatthew G. Knepley       suffix: adv_2d_tri_1
166630602db0SMatthew G. Knepley       requires: triangle
166730602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166830602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -simplex -dm_refine 5 -dm_plex_faces 1,1 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
166930602db0SMatthew G. Knepley 
167030602db0SMatthew G. Knepley     test:
167130602db0SMatthew G. Knepley       suffix: tut_1
167230602db0SMatthew G. Knepley       requires: exodusii
167330602db0SMatthew G. Knepley       nsize: 1
167446ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -
167530602db0SMatthew G. Knepley 
167630602db0SMatthew G. Knepley     test:
167730602db0SMatthew G. Knepley       suffix: tut_2
167830602db0SMatthew G. Knepley       requires: exodusii
167930602db0SMatthew G. Knepley       nsize: 1
168046ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw
168130602db0SMatthew G. Knepley 
168230602db0SMatthew G. Knepley     test:
168330602db0SMatthew G. Knepley       suffix: tut_3
168430602db0SMatthew G. Knepley       requires: exodusii
168530602db0SMatthew G. Knepley       nsize: 4
1686e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -monitor Error -advect_sol_type bump -petscfv_type leastsquares -petsclimiter_type sin
168730602db0SMatthew G. Knepley 
168830602db0SMatthew G. Knepley     test:
168930602db0SMatthew G. Knepley       suffix: tut_4
169030602db0SMatthew G. Knepley       requires: exodusii
169130602db0SMatthew G. Knepley       nsize: 4
1692e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -physics sw -monitor Height,Energy -petscfv_type leastsquares -petsclimiter_type minmod
169330602db0SMatthew G. Knepley 
169430602db0SMatthew G. Knepley   testset:
169530602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
169630602db0SMatthew G. Knepley 
1697c4762a1bSJed Brown     # 2D Advection 0-10
1698c4762a1bSJed Brown     test:
1699c4762a1bSJed Brown       suffix: 0
1700c4762a1bSJed Brown       requires: exodusii
170146ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1702c4762a1bSJed Brown 
1703c4762a1bSJed Brown     test:
1704c4762a1bSJed Brown       suffix: 1
1705c4762a1bSJed Brown       requires: exodusii
170630602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1707c4762a1bSJed Brown 
1708c4762a1bSJed Brown     test:
1709c4762a1bSJed Brown       suffix: 2
1710c4762a1bSJed Brown       requires: exodusii
1711c4762a1bSJed Brown       nsize: 2
171246ac1a18SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1713c4762a1bSJed Brown 
1714c4762a1bSJed Brown     test:
1715c4762a1bSJed Brown       suffix: 3
1716c4762a1bSJed Brown       requires: exodusii
1717c4762a1bSJed Brown       nsize: 2
1718e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1719c4762a1bSJed Brown 
1720c4762a1bSJed Brown     test:
1721c4762a1bSJed Brown       suffix: 4
1722c4762a1bSJed Brown       requires: exodusii
17236c2b77d5SStefano Zampini       nsize: 4
17246c2b77d5SStefano Zampini       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple
1725c4762a1bSJed Brown 
1726c4762a1bSJed Brown     test:
1727c4762a1bSJed Brown       suffix: 5
1728c4762a1bSJed Brown       requires: exodusii
172946ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1
1730c4762a1bSJed Brown 
1731c4762a1bSJed Brown     test:
1732c4762a1bSJed Brown       suffix: 7
1733c4762a1bSJed Brown       requires: exodusii
173430602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1735c4762a1bSJed Brown 
1736c4762a1bSJed Brown     test:
1737c4762a1bSJed Brown       suffix: 8
1738c4762a1bSJed Brown       requires: exodusii
1739c4762a1bSJed Brown       nsize: 2
1740e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1741c4762a1bSJed Brown 
1742c4762a1bSJed Brown     test:
1743c4762a1bSJed Brown       suffix: 9
1744c4762a1bSJed Brown       requires: exodusii
1745c4762a1bSJed Brown       nsize: 8
1746e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1747c4762a1bSJed Brown 
1748c4762a1bSJed Brown     test:
1749c4762a1bSJed Brown       suffix: 10
1750c4762a1bSJed Brown       requires: exodusii
175130602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo
1752c4762a1bSJed Brown 
1753c4762a1bSJed Brown   # 2D Shallow water
1754993a5519SMatthew G. Knepley   testset:
1755993a5519SMatthew G. Knepley     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0
1756993a5519SMatthew G. Knepley 
1757c4762a1bSJed Brown     test:
1758c4762a1bSJed Brown       suffix: sw_0
1759c4762a1bSJed Brown       requires: exodusii
1760993a5519SMatthew G. Knepley       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1761993a5519SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1762993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1763993a5519SMatthew G. Knepley             -monitor height,energy
1764993a5519SMatthew G. Knepley 
1765993a5519SMatthew G. Knepley     test:
1766de8de29fSMatthew G. Knepley       suffix: sw_ceed
1767de8de29fSMatthew G. Knepley       requires: exodusii libceed
1768de8de29fSMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1769de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1770de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1771de8de29fSMatthew G. Knepley             -monitor height,energy
1772de8de29fSMatthew G. Knepley 
1773de8de29fSMatthew G. Knepley     test:
1774a7d931c6SMatthew G. Knepley       suffix: sw_ceed_small
1775a7d931c6SMatthew G. Knepley       requires: exodusii libceed
1776a7d931c6SMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \
1777a7d931c6SMatthew G. Knepley             -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_lower 0,1 -dm_plex_box_upper 6.28,3 -dm_plex_box_faces 8,2 \
1778a7d931c6SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1779a7d931c6SMatthew G. Knepley             -monitor height,energy
1780a7d931c6SMatthew G. Knepley 
1781a7d931c6SMatthew G. Knepley     test:
1782993a5519SMatthew G. Knepley       suffix: sw_1
1783993a5519SMatthew G. Knepley       nsize: 2
1784993a5519SMatthew G. Knepley       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1785993a5519SMatthew G. Knepley             -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_faces 24,12 -dm_plex_box_lower 0,1 -dm_plex_box_upper 1,3 -dm_distribute_overlap 1 \
1786993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1787993a5519SMatthew G. Knepley             -monitor height,energy
1788c4762a1bSJed Brown 
17898d2c18caSMukkund Sunjii     test:
17908d2c18caSMukkund Sunjii       suffix: sw_hll
1791993a5519SMatthew G. Knepley       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1792993a5519SMatthew G. Knepley             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1793993a5519SMatthew G. Knepley             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1794993a5519SMatthew G. Knepley             -monitor height,energy
1795993a5519SMatthew G. Knepley 
1796de8de29fSMatthew G. Knepley   # 2D Euler
1797de8de29fSMatthew G. Knepley   testset:
1798de8de29fSMatthew G. Knepley     args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1799de8de29fSMatthew G. Knepley           -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy
1800de8de29fSMatthew G. Knepley 
1801de8de29fSMatthew G. Knepley     test:
1802277e32beSMatthew G. Knepley       suffix: euler_0
1803277e32beSMatthew G. Knepley       requires: exodusii !complex
1804277e32beSMatthew G. Knepley       args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1805277e32beSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1806277e32beSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1807277e32beSMatthew G. Knepley 
1808277e32beSMatthew G. Knepley     test:
1809de8de29fSMatthew G. Knepley       suffix: euler_ceed
1810de8de29fSMatthew G. Knepley       requires: exodusii libceed
1811de8de29fSMatthew G. Knepley       args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1812de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1813de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1814de8de29fSMatthew G. Knepley 
1815993a5519SMatthew G. Knepley   testset:
1816993a5519SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
18178d2c18caSMukkund Sunjii 
1818c4762a1bSJed Brown     # 2D Advection: p4est
1819c4762a1bSJed Brown     test:
1820c4762a1bSJed Brown       suffix: p4est_advec_2d
1821c4762a1bSJed Brown       requires: p4est
182230602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_type p4est -dm_forest_minimum_refinement 1 -dm_forest_initial_refinement 2 -dm_p4est_refine_pattern hash -dm_forest_maximum_refinement 5
1823c4762a1bSJed Brown 
1824c4762a1bSJed Brown     # Advection in a box
1825c4762a1bSJed Brown     test:
1826c4762a1bSJed Brown       suffix: adv_2d_quad_0
1827c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1828c4762a1bSJed Brown 
1829c4762a1bSJed Brown     test:
1830c4762a1bSJed Brown       suffix: adv_2d_quad_1
1831c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
1832c4762a1bSJed Brown       timeoutfactor: 3
1833c4762a1bSJed Brown 
1834c4762a1bSJed Brown     test:
1835c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_0
1836c4762a1bSJed Brown       requires: p4est
1837c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1838c4762a1bSJed Brown 
1839c4762a1bSJed Brown     test:
1840c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_1
1841c4762a1bSJed Brown       requires: p4est
1842c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
1843c4762a1bSJed Brown       timeoutfactor: 3
1844c4762a1bSJed Brown 
184553df731dSPierre Jolivet     test: # broken for quad precision
1846c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_adapt_0
184753df731dSPierre Jolivet       requires: p4est !__float128
1848c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 -ufv_use_amr -refine_vec_tagger_box 0.005,inf -coarsen_vec_tagger_box 0,1.e-5 -petscfv_type leastsquares -ts_max_time 0.01
1849c4762a1bSJed Brown       timeoutfactor: 3
1850c4762a1bSJed Brown 
1851c4762a1bSJed Brown     test:
1852c4762a1bSJed Brown       suffix: adv_0
1853c4762a1bSJed Brown       requires: exodusii
185446ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -bc_inflow 100,101,200 -bc_outflow 201
1855c4762a1bSJed Brown 
1856c4762a1bSJed Brown     test:
1857c4762a1bSJed Brown       suffix: shock_0
1858c4762a1bSJed Brown       requires: p4est !single !complex
185930602db0SMatthew G. Knepley       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
186030602db0SMatthew G. Knepley       -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 6 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
186130602db0SMatthew G. Knepley       -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
186230602db0SMatthew G. Knepley       -bc_wall 1,2,3,4 -physics euler -eu_type iv_shock -ufv_cfl 10 -eu_alpha 60. -eu_gamma 1.4 -eu_amach 2.02 -eu_rho2 3. \
186330602db0SMatthew G. Knepley       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
186430602db0SMatthew G. Knepley       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
186530602db0SMatthew G. Knepley       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
1866c4762a1bSJed Brown       timeoutfactor: 3
1867c4762a1bSJed Brown 
1868c4762a1bSJed Brown     # Test GLVis visualization of PetscFV fields
1869c4762a1bSJed Brown     test:
1870c4762a1bSJed Brown       suffix: glvis_adv_2d_tet
187130602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
187230602db0SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
187330602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1874c4762a1bSJed Brown 
1875c4762a1bSJed Brown     test:
1876c4762a1bSJed Brown       suffix: glvis_adv_2d_quad
187730602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
187830602db0SMatthew G. Knepley             -dm_refine 5 -dm_plex_separate_marker \
187930602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1880c4762a1bSJed Brown 
1881c4762a1bSJed Brown TEST*/
1882