xref: /petsc/src/ts/tutorials/ex11.c (revision a7d931c6983080384301a9270cd0ca2d035918fb) 
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));
457c4762a1bSJed Brown 
4583ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
459c4762a1bSJed Brown }
460c4762a1bSJed Brown 
461c4762a1bSJed Brown /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
462c4762a1bSJed Brown /* An initial-value and self-similar solutions of the compressible Euler equations */
463c4762a1bSJed Brown /* Ravi Samtaney and D. I. Pullin */
464c4762a1bSJed Brown /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
4659371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_Euler[] = {
4669371c9d4SSatish Balay   {"Density",  1  },
4679371c9d4SSatish Balay   {"Momentum", DIM},
4689371c9d4SSatish Balay   {"Energy",   1  },
4699371c9d4SSatish Balay   {NULL,       0  }
4709371c9d4SSatish Balay };
471c4762a1bSJed Brown 
472c4762a1bSJed Brown /* initial condition */
473c4762a1bSJed Brown int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
474d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
475d71ae5a4SJacob Faibussowitsch {
476c4762a1bSJed Brown   PetscInt       i;
477c4762a1bSJed Brown   Physics        phys = (Physics)ctx;
478c4762a1bSJed Brown   Physics_Euler *eu   = (Physics_Euler *)phys->data;
479c4762a1bSJed Brown   EulerNode     *uu   = (EulerNode *)u;
480de8de29fSMatthew G. Knepley   PetscReal      p0, gamma, c = 0.;
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;
531c4762a1bSJed Brown 
5323ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
533c4762a1bSJed Brown }
534c4762a1bSJed Brown 
535c4762a1bSJed Brown /* PetscReal* => EulerNode* conversion */
536d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
537d71ae5a4SJacob Faibussowitsch {
538c4762a1bSJed Brown   PetscInt         i;
539c4762a1bSJed Brown   const EulerNode *xI   = (const EulerNode *)a_xI;
540c4762a1bSJed Brown   EulerNode       *xG   = (EulerNode *)a_xG;
541c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
542c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
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));
701c4762a1bSJed Brown 
7023ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
703c4762a1bSJed Brown }
704c4762a1bSJed Brown 
705d71ae5a4SJacob Faibussowitsch static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx)
706d71ae5a4SJacob Faibussowitsch {
707c4762a1bSJed Brown   PetscReal err = 0.;
708c4762a1bSJed Brown   PetscInt  i, j;
709c4762a1bSJed Brown 
710c4762a1bSJed Brown   PetscFunctionBeginUser;
711c4762a1bSJed Brown   for (i = 0; i < numComps; i++) {
712ad540459SPierre Jolivet     for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
713c4762a1bSJed Brown   }
714c4762a1bSJed Brown   *error = volume * err;
7153ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
716c4762a1bSJed Brown }
717c4762a1bSJed Brown 
718d71ae5a4SJacob Faibussowitsch PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
719d71ae5a4SJacob Faibussowitsch {
720c4762a1bSJed Brown   PetscSF      sfPoint;
721c4762a1bSJed Brown   PetscSection coordSection;
722c4762a1bSJed Brown   Vec          coordinates;
723c4762a1bSJed Brown   PetscSection sectionCell;
724c4762a1bSJed Brown   PetscScalar *part;
725c4762a1bSJed Brown   PetscInt     cStart, cEnd, c;
726c4762a1bSJed Brown   PetscMPIInt  rank;
727c4762a1bSJed Brown 
728c4762a1bSJed Brown   PetscFunctionBeginUser;
7299566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7309566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7319566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, dmCell));
7329566063dSJacob Faibussowitsch   PetscCall(DMGetPointSF(dm, &sfPoint));
7339566063dSJacob Faibussowitsch   PetscCall(DMSetPointSF(*dmCell, sfPoint));
7349566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
7359566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
7369566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
7379566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
7389566063dSJacob Faibussowitsch   PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
7399566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
74048a46eb9SPierre Jolivet   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
7419566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionCell));
7429566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(*dmCell, sectionCell));
7439566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionCell));
7449566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(*dmCell, partition));
7459566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
7469566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*partition, &part));
747c4762a1bSJed Brown   for (c = cStart; c < cEnd; ++c) {
748c4762a1bSJed Brown     PetscScalar *p;
749c4762a1bSJed Brown 
7509566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
751c4762a1bSJed Brown     p[0] = rank;
752c4762a1bSJed Brown   }
7539566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*partition, &part));
7543ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
755c4762a1bSJed Brown }
756c4762a1bSJed Brown 
757d71ae5a4SJacob Faibussowitsch PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
758d71ae5a4SJacob Faibussowitsch {
7593e9753d6SMatthew G. Knepley   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
760c4762a1bSJed Brown   PetscSection       coordSection;
761c4762a1bSJed Brown   Vec                coordinates, facegeom, cellgeom;
762c4762a1bSJed Brown   PetscSection       sectionMass;
763c4762a1bSJed Brown   PetscScalar       *m;
764c4762a1bSJed Brown   const PetscScalar *fgeom, *cgeom, *coords;
765c4762a1bSJed Brown   PetscInt           vStart, vEnd, v;
766c4762a1bSJed Brown 
767c4762a1bSJed Brown   PetscFunctionBeginUser;
7689566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
7699566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7709566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7719566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, &dmMass));
7729566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
7739566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
7749566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
7759566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
7769566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
777c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
778c4762a1bSJed Brown     PetscInt numFaces;
779c4762a1bSJed Brown 
7809566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
7819566063dSJacob Faibussowitsch     PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
782c4762a1bSJed Brown   }
7839566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionMass));
7849566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(dmMass, sectionMass));
7859566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionMass));
7869566063dSJacob Faibussowitsch   PetscCall(DMGetLocalVector(dmMass, massMatrix));
7879566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*massMatrix, &m));
7889566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
7899566063dSJacob Faibussowitsch   PetscCall(VecGetDM(facegeom, &dmFace));
7909566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(facegeom, &fgeom));
7919566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellgeom, &dmCell));
7929566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
7939566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
7949566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(coordinates, &coords));
795c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
796c4762a1bSJed Brown     const PetscInt  *faces;
797c4762a1bSJed Brown     PetscFVFaceGeom *fgA, *fgB, *cg;
798c4762a1bSJed Brown     PetscScalar     *vertex;
799c4762a1bSJed Brown     PetscInt         numFaces, sides[2], f, g;
800c4762a1bSJed Brown 
8019566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
8029566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
8039566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupport(dmMass, v, &faces));
804c4762a1bSJed Brown     for (f = 0; f < numFaces; ++f) {
805c4762a1bSJed Brown       sides[0] = faces[f];
8069566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
807c4762a1bSJed Brown       for (g = 0; g < numFaces; ++g) {
808c4762a1bSJed Brown         const PetscInt *cells = NULL;
809c4762a1bSJed Brown         PetscReal       area  = 0.0;
810c4762a1bSJed Brown         PetscInt        numCells;
811c4762a1bSJed Brown 
812c4762a1bSJed Brown         sides[1] = faces[g];
8139566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
8149566063dSJacob Faibussowitsch         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
81554c59aa7SJacob Faibussowitsch         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
8169566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
817c4762a1bSJed 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]));
818c4762a1bSJed 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]));
819c4762a1bSJed Brown         m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
8209566063dSJacob Faibussowitsch         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
821c4762a1bSJed Brown       }
822c4762a1bSJed Brown     }
823c4762a1bSJed Brown   }
8249566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
8259566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
8269566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(coordinates, &coords));
8279566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*massMatrix, &m));
8289566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmMass));
8299566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
8303ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
831c4762a1bSJed Brown }
832c4762a1bSJed Brown 
833c4762a1bSJed Brown /* Behavior will be different for multi-physics or when using non-default boundary conditions */
834d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
835d71ae5a4SJacob Faibussowitsch {
836c4762a1bSJed Brown   PetscFunctionBeginUser;
837c4762a1bSJed Brown   mod->solution    = func;
838c4762a1bSJed Brown   mod->solutionctx = ctx;
8393ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
840c4762a1bSJed Brown }
841c4762a1bSJed Brown 
842d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
843d71ae5a4SJacob Faibussowitsch {
844c4762a1bSJed Brown   FunctionalLink link, *ptr;
845c4762a1bSJed Brown   PetscInt       lastoffset = -1;
846c4762a1bSJed Brown 
847c4762a1bSJed Brown   PetscFunctionBeginUser;
848c4762a1bSJed Brown   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
8499566063dSJacob Faibussowitsch   PetscCall(PetscNew(&link));
8509566063dSJacob Faibussowitsch   PetscCall(PetscStrallocpy(name, &link->name));
851c4762a1bSJed Brown   link->offset = lastoffset + 1;
852c4762a1bSJed Brown   link->func   = func;
853c4762a1bSJed Brown   link->ctx    = ctx;
854c4762a1bSJed Brown   link->next   = NULL;
855c4762a1bSJed Brown   *ptr         = link;
856c4762a1bSJed Brown   *offset      = link->offset;
8573ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
858c4762a1bSJed Brown }
859c4762a1bSJed Brown 
860d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject)
861d71ae5a4SJacob Faibussowitsch {
862c4762a1bSJed Brown   PetscInt       i, j;
863c4762a1bSJed Brown   FunctionalLink link;
864c4762a1bSJed Brown   char          *names[256];
865c4762a1bSJed Brown 
866c4762a1bSJed Brown   PetscFunctionBeginUser;
867dd39110bSPierre Jolivet   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
8689566063dSJacob Faibussowitsch   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
869c4762a1bSJed Brown   /* Create list of functionals that will be computed somehow */
8709566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
871c4762a1bSJed Brown   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
8729566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
873c4762a1bSJed Brown   mod->numCall = 0;
874c4762a1bSJed Brown   for (i = 0; i < mod->numMonitored; i++) {
875c4762a1bSJed Brown     for (link = mod->functionalRegistry; link; link = link->next) {
876c4762a1bSJed Brown       PetscBool match;
8779566063dSJacob Faibussowitsch       PetscCall(PetscStrcasecmp(names[i], link->name, &match));
878c4762a1bSJed Brown       if (match) break;
879c4762a1bSJed Brown     }
88054c59aa7SJacob Faibussowitsch     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
881c4762a1bSJed Brown     mod->functionalMonitored[i] = link;
882c4762a1bSJed Brown     for (j = 0; j < i; j++) {
883c4762a1bSJed Brown       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
884c4762a1bSJed Brown     }
885c4762a1bSJed Brown     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
886c4762a1bSJed Brown   next_name:
8879566063dSJacob Faibussowitsch     PetscCall(PetscFree(names[i]));
888c4762a1bSJed Brown   }
889c4762a1bSJed Brown 
890c4762a1bSJed Brown   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
891c4762a1bSJed Brown   mod->maxComputed = -1;
892c4762a1bSJed Brown   for (link = mod->functionalRegistry; link; link = link->next) {
893c4762a1bSJed Brown     for (i = 0; i < mod->numCall; i++) {
894c4762a1bSJed Brown       FunctionalLink call = mod->functionalCall[i];
895ad540459SPierre Jolivet       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
896c4762a1bSJed Brown     }
897c4762a1bSJed Brown   }
8983ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
899c4762a1bSJed Brown }
900c4762a1bSJed Brown 
901d71ae5a4SJacob Faibussowitsch static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
902d71ae5a4SJacob Faibussowitsch {
903c4762a1bSJed Brown   FunctionalLink l, next;
904c4762a1bSJed Brown 
905c4762a1bSJed Brown   PetscFunctionBeginUser;
9063ba16761SJacob Faibussowitsch   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
907c4762a1bSJed Brown   l     = *link;
908c4762a1bSJed Brown   *link = NULL;
909c4762a1bSJed Brown   for (; l; l = next) {
910c4762a1bSJed Brown     next = l->next;
9119566063dSJacob Faibussowitsch     PetscCall(PetscFree(l->name));
9129566063dSJacob Faibussowitsch     PetscCall(PetscFree(l));
913c4762a1bSJed Brown   }
9143ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
915c4762a1bSJed Brown }
916c4762a1bSJed Brown 
917c4762a1bSJed Brown /* put the solution callback into a functional callback */
918d71ae5a4SJacob Faibussowitsch static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
919d71ae5a4SJacob Faibussowitsch {
920c4762a1bSJed Brown   Model mod;
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 
123052115386SStefano Zampini   PetscFunctionReturn(PETSC_SUCCESS);
123152115386SStefano Zampini }
123252115386SStefano Zampini 
1233d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
1234d71ae5a4SJacob Faibussowitsch {
1235c4762a1bSJed Brown   MPI_Comm          comm;
1236c4762a1bSJed Brown   PetscDS           prob;
1237c4762a1bSJed Brown   PetscFV           fvm;
1238c4762a1bSJed Brown   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1239c4762a1bSJed Brown   User              user;
1240c4762a1bSJed Brown   Model             mod;
1241c4762a1bSJed Brown   Physics           phys;
12423e9753d6SMatthew G. Knepley   DM                dm, plex;
1243c4762a1bSJed Brown   PetscReal         ftime, cfl, dt, minRadius;
1244c4762a1bSJed Brown   PetscInt          dim, nsteps;
1245c4762a1bSJed Brown   TS                ts;
1246c4762a1bSJed Brown   TSConvergedReason reason;
1247c4762a1bSJed Brown   Vec               X;
1248c4762a1bSJed Brown   PetscViewer       viewer;
1249b5a892a1SMatthew G. Knepley   PetscBool         vtkCellGeom, useAMR;
125030602db0SMatthew G. Knepley   PetscInt          adaptInterval;
1251c4762a1bSJed Brown   char              physname[256] = "advect";
1252c4762a1bSJed Brown   VecTagger         refineTag = NULL, coarsenTag = NULL;
125352115386SStefano Zampini   TransferCtx       tctx;
1254c4762a1bSJed Brown 
1255327415f7SBarry Smith   PetscFunctionBeginUser;
12569566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1257c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
1258c4762a1bSJed Brown 
12599566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user));
12609566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model));
12619566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model->physics));
1262c4762a1bSJed Brown   mod           = user->model;
1263c4762a1bSJed Brown   phys          = mod->physics;
1264c4762a1bSJed Brown   mod->comm     = comm;
1265c4762a1bSJed Brown   useAMR        = PETSC_FALSE;
1266c4762a1bSJed Brown   adaptInterval = 1;
1267c4762a1bSJed Brown 
1268c4762a1bSJed Brown   /* Register physical models to be available on the command line */
12699566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
12709566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
12719566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1272c4762a1bSJed Brown 
1273d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1274c4762a1bSJed Brown   {
1275c4762a1bSJed Brown     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
12769566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1277c4762a1bSJed Brown     user->vtkInterval = 1;
12789566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1279c4762a1bSJed Brown     user->vtkmon = PETSC_TRUE;
12809566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1281c4762a1bSJed Brown     vtkCellGeom = PETSC_FALSE;
1282c6a7a370SJeremy L Thompson     PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
12839566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
12849566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
12859566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
12869566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1287c4762a1bSJed Brown   }
1288d0609cedSBarry Smith   PetscOptionsEnd();
1289c4762a1bSJed Brown 
1290c4762a1bSJed Brown   if (useAMR) {
1291c4762a1bSJed Brown     VecTaggerBox refineBox, coarsenBox;
1292c4762a1bSJed Brown 
1293c4762a1bSJed Brown     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1294c4762a1bSJed Brown     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;
1295c4762a1bSJed Brown 
12969566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &refineTag));
12979566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
12989566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
12999566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
13009566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(refineTag));
13019566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(refineTag));
13029566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
1303c4762a1bSJed Brown 
13049566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &coarsenTag));
13059566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
13069566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
13079566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
13089566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(coarsenTag));
13099566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(coarsenTag));
13109566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1311c4762a1bSJed Brown   }
1312c4762a1bSJed Brown 
1313d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1314c4762a1bSJed Brown   {
1315c4762a1bSJed Brown     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
13169566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
13179566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
13189566063dSJacob Faibussowitsch     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
13199566063dSJacob Faibussowitsch     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1320c4762a1bSJed Brown     /* Count number of fields and dofs */
1321c4762a1bSJed Brown     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
132254c59aa7SJacob Faibussowitsch     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
13239566063dSJacob Faibussowitsch     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1324c4762a1bSJed Brown   }
1325d0609cedSBarry Smith   PetscOptionsEnd();
1326c4762a1bSJed Brown 
1327c4762a1bSJed Brown   /* Create mesh */
1328c4762a1bSJed Brown   {
132930602db0SMatthew G. Knepley     PetscInt i;
133030602db0SMatthew G. Knepley 
13319566063dSJacob Faibussowitsch     PetscCall(DMCreate(comm, &dm));
13329566063dSJacob Faibussowitsch     PetscCall(DMSetType(dm, DMPLEX));
13339566063dSJacob Faibussowitsch     PetscCall(DMSetFromOptions(dm));
13349371c9d4SSatish Balay     for (i = 0; i < DIM; i++) {
13359371c9d4SSatish Balay       mod->bounds[2 * i]     = 0.;
13369371c9d4SSatish Balay       mod->bounds[2 * i + 1] = 1.;
13379371c9d4SSatish Balay     };
1338c4762a1bSJed Brown     dim = DIM;
133930602db0SMatthew G. Knepley     { /* a null name means just do a hex box */
134030602db0SMatthew G. Knepley       PetscInt  cells[3] = {1, 1, 1}, n = 3;
134130602db0SMatthew G. Knepley       PetscBool flg2, skew              = PETSC_FALSE;
1342c4762a1bSJed Brown       PetscInt  nret2 = 2 * DIM;
1343d0609cedSBarry Smith       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
13449566063dSJacob 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));
13459566063dSJacob Faibussowitsch       PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
13469566063dSJacob Faibussowitsch       PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1347d0609cedSBarry Smith       PetscOptionsEnd();
134830602db0SMatthew G. Knepley       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1349c4762a1bSJed Brown       if (flg2) {
1350c4762a1bSJed Brown         PetscInt     dimEmbed, i;
1351c4762a1bSJed Brown         PetscInt     nCoords;
1352c4762a1bSJed Brown         PetscScalar *coords;
1353c4762a1bSJed Brown         Vec          coordinates;
1354c4762a1bSJed Brown 
13559566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
13569566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
13579566063dSJacob Faibussowitsch         PetscCall(VecGetLocalSize(coordinates, &nCoords));
135854c59aa7SJacob Faibussowitsch         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
13599566063dSJacob Faibussowitsch         PetscCall(VecGetArray(coordinates, &coords));
1360c4762a1bSJed Brown         for (i = 0; i < nCoords; i += dimEmbed) {
1361c4762a1bSJed Brown           PetscInt j;
1362c4762a1bSJed Brown 
1363c4762a1bSJed Brown           PetscScalar *coord = &coords[i];
1364c4762a1bSJed Brown           for (j = 0; j < dimEmbed; j++) {
1365c4762a1bSJed Brown             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1366c4762a1bSJed Brown             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1367c4762a1bSJed Brown               if (cells[0] == 2 && i == 8) {
1368c4762a1bSJed Brown                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
136930602db0SMatthew G. Knepley               } else if (cells[0] == 3) {
1370c4762a1bSJed Brown                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1371c4762a1bSJed Brown                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1372c4762a1bSJed Brown                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1373c4762a1bSJed Brown               }
1374c4762a1bSJed Brown             }
1375c4762a1bSJed Brown           }
1376c4762a1bSJed Brown         }
13779566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(coordinates, &coords));
13789566063dSJacob Faibussowitsch         PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1379c4762a1bSJed Brown       }
1380c4762a1bSJed Brown     }
1381c4762a1bSJed Brown   }
13829566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
13839566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
1384c4762a1bSJed Brown 
1385c4762a1bSJed Brown   /* set up BCs, functions, tags */
13869566063dSJacob Faibussowitsch   PetscCall(DMCreateLabel(dm, "Face Sets"));
1387c4762a1bSJed Brown   mod->errorIndicator = ErrorIndicator_Simple;
1388c4762a1bSJed Brown 
1389c4762a1bSJed Brown   {
1390c4762a1bSJed Brown     DM gdm;
1391c4762a1bSJed Brown 
13929566063dSJacob Faibussowitsch     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
13939566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dm));
1394c4762a1bSJed Brown     dm = gdm;
13959566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1396c4762a1bSJed Brown   }
1397c4762a1bSJed Brown 
13989566063dSJacob Faibussowitsch   PetscCall(PetscFVCreate(comm, &fvm));
13999566063dSJacob Faibussowitsch   PetscCall(PetscFVSetFromOptions(fvm));
14009566063dSJacob Faibussowitsch   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
14019566063dSJacob Faibussowitsch   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
14029566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1403c4762a1bSJed Brown   {
1404c4762a1bSJed Brown     PetscInt f, dof;
1405c4762a1bSJed Brown     for (f = 0, dof = 0; f < phys->nfields; f++) {
1406c4762a1bSJed Brown       PetscInt newDof = phys->field_desc[f].dof;
1407c4762a1bSJed Brown 
1408c4762a1bSJed Brown       if (newDof == 1) {
14099566063dSJacob Faibussowitsch         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
14109371c9d4SSatish Balay       } else {
1411c4762a1bSJed Brown         PetscInt j;
1412c4762a1bSJed Brown 
1413c4762a1bSJed Brown         for (j = 0; j < newDof; j++) {
1414c4762a1bSJed Brown           char compName[256] = "Unknown";
1415c4762a1bSJed Brown 
141663a3b9bcSJacob Faibussowitsch           PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
14179566063dSJacob Faibussowitsch           PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1418c4762a1bSJed Brown         }
1419c4762a1bSJed Brown       }
1420c4762a1bSJed Brown       dof += newDof;
1421c4762a1bSJed Brown     }
1422c4762a1bSJed Brown   }
1423c4762a1bSJed Brown   /* FV is now structured with one field having all physics as components */
14249566063dSJacob Faibussowitsch   PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
14259566063dSJacob Faibussowitsch   PetscCall(DMCreateDS(dm));
14269566063dSJacob Faibussowitsch   PetscCall(DMGetDS(dm, &prob));
14279566063dSJacob Faibussowitsch   PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
14289566063dSJacob Faibussowitsch   PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
14299566063dSJacob Faibussowitsch   PetscCall((*mod->setupbc)(dm, prob, phys));
14309566063dSJacob Faibussowitsch   PetscCall(PetscDSSetFromOptions(prob));
1431c4762a1bSJed Brown   {
1432c4762a1bSJed Brown     char      convType[256];
1433c4762a1bSJed Brown     PetscBool flg;
1434c4762a1bSJed Brown 
1435d0609cedSBarry Smith     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
14369566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1437d0609cedSBarry Smith     PetscOptionsEnd();
1438c4762a1bSJed Brown     if (flg) {
1439c4762a1bSJed Brown       DM dmConv;
1440c4762a1bSJed Brown 
14419566063dSJacob Faibussowitsch       PetscCall(DMConvert(dm, convType, &dmConv));
1442c4762a1bSJed Brown       if (dmConv) {
14439566063dSJacob Faibussowitsch         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
14449566063dSJacob Faibussowitsch         PetscCall(DMDestroy(&dm));
1445c4762a1bSJed Brown         dm = dmConv;
14469566063dSJacob Faibussowitsch         PetscCall(DMSetFromOptions(dm));
1447c4762a1bSJed Brown       }
1448c4762a1bSJed Brown     }
1449c4762a1bSJed Brown   }
1450de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
1451de8de29fSMatthew G. Knepley   {
1452de8de29fSMatthew G. Knepley     PetscBool useCeed;
1453de8de29fSMatthew G. Knepley     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1454de8de29fSMatthew G. Knepley     if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1455de8de29fSMatthew G. Knepley   }
1456de8de29fSMatthew G. Knepley #endif
1457c4762a1bSJed Brown 
14589566063dSJacob Faibussowitsch   PetscCall(initializeTS(dm, user, &ts));
1459c4762a1bSJed Brown 
14609566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(dm, &X));
14619566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
14629566063dSJacob Faibussowitsch   PetscCall(SetInitialCondition(dm, X, user));
1463c4762a1bSJed Brown   if (useAMR) {
1464c4762a1bSJed Brown     PetscInt adaptIter;
1465c4762a1bSJed Brown 
1466c4762a1bSJed Brown     /* use no limiting when reconstructing gradients for adaptivity */
14679566063dSJacob Faibussowitsch     PetscCall(PetscFVGetLimiter(fvm, &limiter));
14689566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)limiter));
14699566063dSJacob Faibussowitsch     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
14709566063dSJacob Faibussowitsch     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
1471c4762a1bSJed Brown 
147252115386SStefano Zampini     /* Refinement context */
147352115386SStefano Zampini     tctx.fvm         = fvm;
147452115386SStefano Zampini     tctx.refineTag   = refineTag;
147552115386SStefano Zampini     tctx.coarsenTag  = coarsenTag;
147652115386SStefano Zampini     tctx.adaptedDM   = NULL;
147752115386SStefano Zampini     tctx.user        = user;
147852115386SStefano Zampini     tctx.noneLimiter = noneLimiter;
147952115386SStefano Zampini     tctx.limiter     = limiter;
148052115386SStefano Zampini     tctx.cfl         = cfl;
148152115386SStefano Zampini 
148252115386SStefano Zampini     /* Do some initial refinement steps */
1483c4762a1bSJed Brown     for (adaptIter = 0;; ++adaptIter) {
1484c4762a1bSJed Brown       PetscLogDouble bytes;
148552115386SStefano Zampini       PetscBool      resize;
1486c4762a1bSJed Brown 
14879566063dSJacob Faibussowitsch       PetscCall(PetscMemoryGetCurrentUsage(&bytes));
148863a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
14899566063dSJacob Faibussowitsch       PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
14909566063dSJacob Faibussowitsch       PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));
1491c4762a1bSJed Brown 
149252115386SStefano Zampini       PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
149352115386SStefano Zampini       if (!resize) break;
14949566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dm));
14959566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&X));
149652115386SStefano Zampini       dm             = tctx.adaptedDM;
149752115386SStefano Zampini       tctx.adaptedDM = NULL;
149852115386SStefano Zampini       PetscCall(TSSetDM(ts, dm));
14999566063dSJacob Faibussowitsch       PetscCall(DMCreateGlobalVector(dm, &X));
15009566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
15019566063dSJacob Faibussowitsch       PetscCall(SetInitialCondition(dm, X, user));
1502c4762a1bSJed Brown     }
1503c4762a1bSJed Brown   }
1504c4762a1bSJed Brown 
15059566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
1506c4762a1bSJed Brown   if (vtkCellGeom) {
1507c4762a1bSJed Brown     DM  dmCell;
1508c4762a1bSJed Brown     Vec cellgeom, partition;
1509c4762a1bSJed Brown 
15109566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
15119566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
15129566063dSJacob Faibussowitsch     PetscCall(VecView(cellgeom, viewer));
15139566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15149566063dSJacob Faibussowitsch     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
15159566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
15169566063dSJacob Faibussowitsch     PetscCall(VecView(partition, viewer));
15179566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15189566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&partition));
15199566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dmCell));
1520c4762a1bSJed Brown   }
1521c4762a1bSJed Brown   /* collect max maxspeed from all processes -- todo */
15229566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
15239566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1524712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
152554c59aa7SJacob Faibussowitsch   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1526c4762a1bSJed Brown   dt = cfl * minRadius / mod->maxspeed;
15279566063dSJacob Faibussowitsch   PetscCall(TSSetTimeStep(ts, dt));
15289566063dSJacob Faibussowitsch   PetscCall(TSSetFromOptions(ts));
1529c4762a1bSJed Brown 
153052115386SStefano Zampini   /* When using adaptive mesh refinement
153152115386SStefano Zampini      specify callbacks to refine the solution
153252115386SStefano Zampini      and interpolate data from old to new mesh */
153352115386SStefano Zampini   if (useAMR) { PetscCall(TSSetResize(ts, adaptToleranceFVMSetUp, Transfer, &tctx)); }
153452115386SStefano Zampini   PetscCall(TSSetSolution(ts, X));
15359566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&X));
153652115386SStefano Zampini   PetscCall(TSSolve(ts, NULL));
15379566063dSJacob Faibussowitsch   PetscCall(TSGetSolveTime(ts, &ftime));
15389566063dSJacob Faibussowitsch   PetscCall(TSGetStepNumber(ts, &nsteps));
15399566063dSJacob Faibussowitsch   PetscCall(TSGetConvergedReason(ts, &reason));
154063a3b9bcSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
15419566063dSJacob Faibussowitsch   PetscCall(TSDestroy(&ts));
1542c4762a1bSJed Brown 
15439566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&refineTag));
15449566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&coarsenTag));
15459566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsList));
15469566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1547de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
15489566063dSJacob Faibussowitsch   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
15499566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalMonitored));
15509566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalCall));
15519566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics->data));
15529566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics));
15539566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model));
15549566063dSJacob Faibussowitsch   PetscCall(PetscFree(user));
15559566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&limiter));
15569566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&noneLimiter));
15579566063dSJacob Faibussowitsch   PetscCall(PetscFVDestroy(&fvm));
15589566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dm));
15599566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
1560b122ec5aSJacob Faibussowitsch   return 0;
1561c4762a1bSJed Brown }
1562c4762a1bSJed Brown 
1563c4762a1bSJed Brown /* Subroutine to set up the initial conditions for the */
1564c4762a1bSJed Brown /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1565c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1566d71ae5a4SJacob Faibussowitsch int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1567d71ae5a4SJacob Faibussowitsch {
1568c4762a1bSJed Brown   int j, k;
1569c4762a1bSJed Brown   /*      Wc=matmul(lv,Ueq) 3 vars */
1570c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1571c4762a1bSJed Brown     wc[k] = 0.;
1572ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1573c4762a1bSJed Brown   }
1574c4762a1bSJed Brown   return 0;
1575c4762a1bSJed Brown }
1576c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1577d71ae5a4SJacob Faibussowitsch int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1578d71ae5a4SJacob Faibussowitsch {
1579c4762a1bSJed Brown   int k, j;
1580c4762a1bSJed Brown   /*      V=matmul(rv,WC) 3 vars */
1581c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1582c4762a1bSJed Brown     v[k] = 0.;
1583ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1584c4762a1bSJed Brown   }
1585c4762a1bSJed Brown   return 0;
1586c4762a1bSJed Brown }
1587c4762a1bSJed Brown /* ---------------------------------------------------------------------- */
1588d71ae5a4SJacob Faibussowitsch int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1589d71ae5a4SJacob Faibussowitsch {
1590c4762a1bSJed Brown   int       j, k;
1591c4762a1bSJed Brown   PetscReal rho, csnd, p0;
1592c4762a1bSJed Brown   /* PetscScalar u; */
1593c4762a1bSJed Brown 
15949371c9d4SSatish Balay   for (k = 0; k < 3; ++k)
15959371c9d4SSatish Balay     for (j = 0; j < 3; ++j) {
15969371c9d4SSatish Balay       lv[k][j] = 0.;
15979371c9d4SSatish Balay       rv[k][j] = 0.;
15989371c9d4SSatish Balay     }
1599c4762a1bSJed Brown   rho = ueq[0];
1600c4762a1bSJed Brown   /* u = ueq[1]; */
1601c4762a1bSJed Brown   p0       = ueq[2];
1602c4762a1bSJed Brown   csnd     = PetscSqrtReal(gamma * p0 / rho);
1603c4762a1bSJed Brown   lv[0][1] = rho * .5;
1604c4762a1bSJed Brown   lv[0][2] = -.5 / csnd;
1605c4762a1bSJed Brown   lv[1][0] = csnd;
1606c4762a1bSJed Brown   lv[1][2] = -1. / csnd;
1607c4762a1bSJed Brown   lv[2][1] = rho * .5;
1608c4762a1bSJed Brown   lv[2][2] = .5 / csnd;
1609c4762a1bSJed Brown   rv[0][0] = -1. / csnd;
1610c4762a1bSJed Brown   rv[1][0] = 1. / rho;
1611c4762a1bSJed Brown   rv[2][0] = -csnd;
1612c4762a1bSJed Brown   rv[0][1] = 1. / csnd;
1613c4762a1bSJed Brown   rv[0][2] = 1. / csnd;
1614c4762a1bSJed Brown   rv[1][2] = 1. / rho;
1615c4762a1bSJed Brown   rv[2][2] = csnd;
1616c4762a1bSJed Brown   return 0;
1617c4762a1bSJed Brown }
1618c4762a1bSJed Brown 
1619d71ae5a4SJacob Faibussowitsch int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1620d71ae5a4SJacob Faibussowitsch {
1621c4762a1bSJed Brown   PetscReal p0, u0, wcp[3], wc[3];
1622c4762a1bSJed Brown   PetscReal lv[3][3];
1623c4762a1bSJed Brown   PetscReal vp[3];
1624c4762a1bSJed Brown   PetscReal rv[3][3];
1625c4762a1bSJed Brown   PetscReal eps, ueq[3], rho0, twopi;
1626c4762a1bSJed Brown 
1627c4762a1bSJed Brown   /* Function Body */
1628c4762a1bSJed Brown   twopi  = 2. * PETSC_PI;
1629c4762a1bSJed Brown   eps    = 1e-4;    /* perturbation */
1630c4762a1bSJed Brown   rho0   = 1e3;     /* density of water */
1631c4762a1bSJed Brown   p0     = 101325.; /* init pressure of 1 atm (?) */
1632c4762a1bSJed Brown   u0     = 0.;
1633c4762a1bSJed Brown   ueq[0] = rho0;
1634c4762a1bSJed Brown   ueq[1] = u0;
1635c4762a1bSJed Brown   ueq[2] = p0;
1636c4762a1bSJed Brown   /* Project initial state to characteristic variables */
1637c4762a1bSJed Brown   eigenvectors(rv, lv, ueq, gamma);
1638c4762a1bSJed Brown   projecteqstate(wc, ueq, lv);
1639c4762a1bSJed Brown   wcp[0] = wc[0];
1640c4762a1bSJed Brown   wcp[1] = wc[1];
1641c4762a1bSJed Brown   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1642c4762a1bSJed Brown   projecttoprim(vp, wcp, rv);
1643c4762a1bSJed Brown   ux->r     = vp[0];         /* density */
1644c4762a1bSJed Brown   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1645c4762a1bSJed Brown   ux->ru[1] = 0.;
1646c4762a1bSJed Brown #if defined DIM > 2
1647c4762a1bSJed Brown   if (dim > 2) ux->ru[2] = 0.;
1648c4762a1bSJed Brown #endif
1649c4762a1bSJed Brown   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1650c4762a1bSJed Brown   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1651c4762a1bSJed Brown   return 0;
1652c4762a1bSJed Brown }
1653c4762a1bSJed Brown 
1654c4762a1bSJed Brown /*TEST
1655c4762a1bSJed Brown 
165630602db0SMatthew G. Knepley   testset:
165730602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0
165830602db0SMatthew G. Knepley 
165930602db0SMatthew G. Knepley     test:
166030602db0SMatthew G. Knepley       suffix: adv_2d_tri_0
166130602db0SMatthew G. Knepley       requires: triangle
166230602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166330602db0SMatthew 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
166430602db0SMatthew G. Knepley 
166530602db0SMatthew G. Knepley     test:
166630602db0SMatthew G. Knepley       suffix: adv_2d_tri_1
166730602db0SMatthew G. Knepley       requires: triangle
166830602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166930602db0SMatthew 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
167030602db0SMatthew G. Knepley 
167130602db0SMatthew G. Knepley     test:
167230602db0SMatthew G. Knepley       suffix: tut_1
167330602db0SMatthew G. Knepley       requires: exodusii
167430602db0SMatthew G. Knepley       nsize: 1
167546ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -
167630602db0SMatthew G. Knepley 
167730602db0SMatthew G. Knepley     test:
167830602db0SMatthew G. Knepley       suffix: tut_2
167930602db0SMatthew G. Knepley       requires: exodusii
168030602db0SMatthew G. Knepley       nsize: 1
168146ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw
168230602db0SMatthew G. Knepley 
168330602db0SMatthew G. Knepley     test:
168430602db0SMatthew G. Knepley       suffix: tut_3
168530602db0SMatthew G. Knepley       requires: exodusii
168630602db0SMatthew G. Knepley       nsize: 4
1687e600fa54SMatthew 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
168830602db0SMatthew G. Knepley 
168930602db0SMatthew G. Knepley     test:
169030602db0SMatthew G. Knepley       suffix: tut_4
169130602db0SMatthew G. Knepley       requires: exodusii
169230602db0SMatthew G. Knepley       nsize: 4
1693e600fa54SMatthew 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
169430602db0SMatthew G. Knepley 
169530602db0SMatthew G. Knepley   testset:
169630602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
169730602db0SMatthew G. Knepley 
1698c4762a1bSJed Brown     # 2D Advection 0-10
1699c4762a1bSJed Brown     test:
1700c4762a1bSJed Brown       suffix: 0
1701c4762a1bSJed Brown       requires: exodusii
170246ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1703c4762a1bSJed Brown 
1704c4762a1bSJed Brown     test:
1705c4762a1bSJed Brown       suffix: 1
1706c4762a1bSJed Brown       requires: exodusii
170730602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1708c4762a1bSJed Brown 
1709c4762a1bSJed Brown     test:
1710c4762a1bSJed Brown       suffix: 2
1711c4762a1bSJed Brown       requires: exodusii
1712c4762a1bSJed Brown       nsize: 2
171346ac1a18SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1714c4762a1bSJed Brown 
1715c4762a1bSJed Brown     test:
1716c4762a1bSJed Brown       suffix: 3
1717c4762a1bSJed Brown       requires: exodusii
1718c4762a1bSJed Brown       nsize: 2
1719e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1720c4762a1bSJed Brown 
1721c4762a1bSJed Brown     test:
1722c4762a1bSJed Brown       suffix: 4
1723c4762a1bSJed Brown       requires: exodusii
17246c2b77d5SStefano Zampini       nsize: 4
17256c2b77d5SStefano 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
1726c4762a1bSJed Brown 
1727c4762a1bSJed Brown     test:
1728c4762a1bSJed Brown       suffix: 5
1729c4762a1bSJed Brown       requires: exodusii
173046ac1a18SMatthew 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
1731c4762a1bSJed Brown 
1732c4762a1bSJed Brown     test:
1733c4762a1bSJed Brown       suffix: 7
1734c4762a1bSJed Brown       requires: exodusii
173530602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1736c4762a1bSJed Brown 
1737c4762a1bSJed Brown     test:
1738c4762a1bSJed Brown       suffix: 8
1739c4762a1bSJed Brown       requires: exodusii
1740c4762a1bSJed Brown       nsize: 2
1741e600fa54SMatthew 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
1742c4762a1bSJed Brown 
1743c4762a1bSJed Brown     test:
1744c4762a1bSJed Brown       suffix: 9
1745c4762a1bSJed Brown       requires: exodusii
1746c4762a1bSJed Brown       nsize: 8
1747e600fa54SMatthew 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
1748c4762a1bSJed Brown 
1749c4762a1bSJed Brown     test:
1750c4762a1bSJed Brown       suffix: 10
1751c4762a1bSJed Brown       requires: exodusii
175230602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo
1753c4762a1bSJed Brown 
1754c4762a1bSJed Brown   # 2D Shallow water
1755993a5519SMatthew G. Knepley   testset:
1756993a5519SMatthew G. Knepley     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0
1757993a5519SMatthew G. Knepley 
1758c4762a1bSJed Brown     test:
1759c4762a1bSJed Brown       suffix: sw_0
1760c4762a1bSJed Brown       requires: exodusii
1761993a5519SMatthew G. Knepley       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1762993a5519SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1763993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1764993a5519SMatthew G. Knepley             -monitor height,energy
1765993a5519SMatthew G. Knepley 
1766993a5519SMatthew G. Knepley     test:
1767de8de29fSMatthew G. Knepley       suffix: sw_ceed
1768de8de29fSMatthew G. Knepley       requires: exodusii libceed
1769de8de29fSMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1770de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1771de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1772de8de29fSMatthew G. Knepley             -monitor height,energy
1773de8de29fSMatthew G. Knepley 
1774de8de29fSMatthew G. Knepley     test:
1775*a7d931c6SMatthew G. Knepley       suffix: sw_ceed_small
1776*a7d931c6SMatthew G. Knepley       requires: exodusii libceed
1777*a7d931c6SMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \
1778*a7d931c6SMatthew 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 \
1779*a7d931c6SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1780*a7d931c6SMatthew G. Knepley             -monitor height,energy
1781*a7d931c6SMatthew G. Knepley 
1782*a7d931c6SMatthew G. Knepley     test:
1783993a5519SMatthew G. Knepley       suffix: sw_1
1784993a5519SMatthew G. Knepley       nsize: 2
1785993a5519SMatthew G. Knepley       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1786993a5519SMatthew 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 \
1787993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1788993a5519SMatthew G. Knepley             -monitor height,energy
1789c4762a1bSJed Brown 
17908d2c18caSMukkund Sunjii     test:
17918d2c18caSMukkund Sunjii       suffix: sw_hll
1792993a5519SMatthew G. Knepley       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1793993a5519SMatthew G. Knepley             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1794993a5519SMatthew G. Knepley             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1795993a5519SMatthew G. Knepley             -monitor height,energy
1796993a5519SMatthew G. Knepley 
1797de8de29fSMatthew G. Knepley   # 2D Euler
1798de8de29fSMatthew G. Knepley   testset:
1799de8de29fSMatthew G. Knepley     args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1800de8de29fSMatthew G. Knepley           -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy
1801de8de29fSMatthew G. Knepley 
1802de8de29fSMatthew G. Knepley     test:
1803277e32beSMatthew G. Knepley       suffix: euler_0
1804277e32beSMatthew G. Knepley       requires: exodusii !complex
1805277e32beSMatthew G. Knepley       args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1806277e32beSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1807277e32beSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1808277e32beSMatthew G. Knepley 
1809277e32beSMatthew G. Knepley     test:
1810de8de29fSMatthew G. Knepley       suffix: euler_ceed
1811de8de29fSMatthew G. Knepley       requires: exodusii libceed
1812de8de29fSMatthew G. Knepley       args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1813de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1814de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1815de8de29fSMatthew G. Knepley 
1816993a5519SMatthew G. Knepley   testset:
1817993a5519SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
18188d2c18caSMukkund Sunjii 
1819c4762a1bSJed Brown     # 2D Advection: p4est
1820c4762a1bSJed Brown     test:
1821c4762a1bSJed Brown       suffix: p4est_advec_2d
1822c4762a1bSJed Brown       requires: p4est
182330602db0SMatthew 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
1824c4762a1bSJed Brown 
1825c4762a1bSJed Brown     # Advection in a box
1826c4762a1bSJed Brown     test:
1827c4762a1bSJed Brown       suffix: adv_2d_quad_0
1828c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1829c4762a1bSJed Brown 
1830c4762a1bSJed Brown     test:
1831c4762a1bSJed Brown       suffix: adv_2d_quad_1
1832c4762a1bSJed 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
1833c4762a1bSJed Brown       timeoutfactor: 3
1834c4762a1bSJed Brown 
1835c4762a1bSJed Brown     test:
1836c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_0
1837c4762a1bSJed Brown       requires: p4est
1838c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1839c4762a1bSJed Brown 
1840c4762a1bSJed Brown     test:
1841c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_1
1842c4762a1bSJed Brown       requires: p4est
1843c4762a1bSJed 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
1844c4762a1bSJed Brown       timeoutfactor: 3
1845c4762a1bSJed Brown 
184653df731dSPierre Jolivet     test: # broken for quad precision
1847c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_adapt_0
184853df731dSPierre Jolivet       requires: p4est !__float128
1849c4762a1bSJed 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
1850c4762a1bSJed Brown       timeoutfactor: 3
1851c4762a1bSJed Brown 
1852c4762a1bSJed Brown     test:
1853c4762a1bSJed Brown       suffix: adv_0
1854c4762a1bSJed Brown       requires: exodusii
185546ac1a18SMatthew 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
1856c4762a1bSJed Brown 
1857c4762a1bSJed Brown     test:
1858c4762a1bSJed Brown       suffix: shock_0
1859c4762a1bSJed Brown       requires: p4est !single !complex
186030602db0SMatthew G. Knepley       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
186130602db0SMatthew 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 \
186230602db0SMatthew 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 \
186330602db0SMatthew 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. \
186430602db0SMatthew G. Knepley       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
186530602db0SMatthew G. Knepley       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
186630602db0SMatthew G. Knepley       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
1867c4762a1bSJed Brown       timeoutfactor: 3
1868c4762a1bSJed Brown 
1869c4762a1bSJed Brown     # Test GLVis visualization of PetscFV fields
1870c4762a1bSJed Brown     test:
1871c4762a1bSJed Brown       suffix: glvis_adv_2d_tet
187230602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
187330602db0SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
187430602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1875c4762a1bSJed Brown 
1876c4762a1bSJed Brown     test:
1877c4762a1bSJed Brown       suffix: glvis_adv_2d_quad
187830602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
187930602db0SMatthew G. Knepley             -dm_refine 5 -dm_plex_separate_marker \
188030602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1881c4762a1bSJed Brown 
1882c4762a1bSJed Brown TEST*/
1883