xref: /petsc/src/ts/tutorials/ex11.c (revision 910b42cbfb1b6524ef4594118dd5e013aee7da5d) 
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"));
402da0802e2SStefano Zampini   PetscFunctionReturn(PETSC_SUCCESS);
403de8de29fSMatthew G. Knepley }
404de8de29fSMatthew G. Knepley #endif
405de8de29fSMatthew G. Knepley 
406de8de29fSMatthew G. Knepley static PetscErrorCode SetupCEED_SW(DM dm, Physics physics)
407de8de29fSMatthew G. Knepley {
408de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
409de8de29fSMatthew G. Knepley   Ceed                 ceed;
410de8de29fSMatthew G. Knepley   CeedQFunctionContext qfCtx;
411de8de29fSMatthew G. Knepley #endif
412de8de29fSMatthew G. Knepley 
413de8de29fSMatthew G. Knepley   PetscFunctionBegin;
414de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
415de8de29fSMatthew G. Knepley   PetscCall(DMGetCeed(dm, &ceed));
416de8de29fSMatthew G. Knepley   PetscCall(CreateQFunctionContext_SW(physics, ceed, &qfCtx));
417de8de29fSMatthew G. Knepley   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_SW_Rusanov_CEED, PhysicsRiemann_SW_Rusanov_CEED_loc, qfCtx));
418de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
419de8de29fSMatthew G. Knepley #endif
420de8de29fSMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
421de8de29fSMatthew G. Knepley }
422de8de29fSMatthew G. Knepley 
423d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
424d71ae5a4SJacob Faibussowitsch {
425c4762a1bSJed Brown   Physics_SW *sw;
4268d2c18caSMukkund Sunjii   char        sw_riemann[64] = "rusanov";
427c4762a1bSJed Brown 
428c4762a1bSJed Brown   PetscFunctionBeginUser;
429c4762a1bSJed Brown   phys->field_desc = PhysicsFields_SW;
4309566063dSJacob Faibussowitsch   PetscCall(PetscNew(&sw));
431c4762a1bSJed Brown   phys->data     = sw;
432c4762a1bSJed Brown   mod->setupbc   = SetUpBC_SW;
433de8de29fSMatthew G. Knepley   mod->setupCEED = SetupCEED_SW;
434c4762a1bSJed Brown 
4353ba16761SJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov));
4363ba16761SJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL));
437de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
438de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov_ceed", PhysicsRiemann_SW_Rusanov_CEED));
439de8de29fSMatthew G. Knepley #endif
4408d2c18caSMukkund Sunjii 
441d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
442c4762a1bSJed Brown   {
4438d2c18caSMukkund Sunjii     void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
444c4762a1bSJed Brown     sw->gravity = 1.0;
4459566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
4469566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL));
4479566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW));
4488d2c18caSMukkund Sunjii     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_SW;
449c4762a1bSJed Brown   }
450d0609cedSBarry Smith   PetscOptionsHeadEnd();
451c4762a1bSJed Brown   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */
452c4762a1bSJed Brown 
4539566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
4549566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
4559566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
4569566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
4573ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
458c4762a1bSJed Brown }
459c4762a1bSJed Brown 
460c4762a1bSJed Brown /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
461c4762a1bSJed Brown /* An initial-value and self-similar solutions of the compressible Euler equations */
462c4762a1bSJed Brown /* Ravi Samtaney and D. I. Pullin */
463c4762a1bSJed Brown /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
4649371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_Euler[] = {
4659371c9d4SSatish Balay   {"Density",  1  },
4669371c9d4SSatish Balay   {"Momentum", DIM},
4679371c9d4SSatish Balay   {"Energy",   1  },
4689371c9d4SSatish Balay   {NULL,       0  }
4699371c9d4SSatish Balay };
470c4762a1bSJed Brown 
471c4762a1bSJed Brown /* initial condition */
472c4762a1bSJed Brown int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
473d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
474d71ae5a4SJacob Faibussowitsch {
475c4762a1bSJed Brown   PetscInt       i;
476c4762a1bSJed Brown   Physics        phys = (Physics)ctx;
477c4762a1bSJed Brown   Physics_Euler *eu   = (Physics_Euler *)phys->data;
478c4762a1bSJed Brown   EulerNode     *uu   = (EulerNode *)u;
479de8de29fSMatthew G. Knepley   PetscReal      p0, gamma, c = 0.;
4804d86920dSPierre Jolivet 
481c4762a1bSJed Brown   PetscFunctionBeginUser;
48254c59aa7SJacob Faibussowitsch   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
483c4762a1bSJed Brown 
484c4762a1bSJed Brown   for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
485c4762a1bSJed Brown   /* set E and rho */
486de8de29fSMatthew G. Knepley   gamma = eu->gamma;
487c4762a1bSJed Brown 
488c4762a1bSJed Brown   if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) {
489c4762a1bSJed Brown     /******************* Euler Density Shock ********************/
490c4762a1bSJed Brown     /* On initial-value and self-similar solutions of the compressible Euler equations */
491c4762a1bSJed Brown     /* Ravi Samtaney and D. I. Pullin */
492c4762a1bSJed Brown     /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
493c4762a1bSJed Brown     /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity,  */
494c4762a1bSJed Brown     p0 = 1.;
495de8de29fSMatthew G. Knepley     if (x[0] < 0.0 + x[1] * eu->itana) {
496c4762a1bSJed Brown       if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */
497c4762a1bSJed Brown         PetscReal amach, rho, press, gas1, p1;
498de8de29fSMatthew G. Knepley         amach     = eu->amach;
499c4762a1bSJed Brown         rho       = 1.;
500c4762a1bSJed Brown         press     = p0;
501c4762a1bSJed Brown         p1        = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0));
502c4762a1bSJed Brown         gas1      = (gamma - 1.0) / (gamma + 1.0);
503c4762a1bSJed Brown         uu->r     = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0);
504c4762a1bSJed Brown         uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach);
505c4762a1bSJed Brown         uu->E     = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0];
5069371c9d4SSatish Balay       } else {      /* left of discontinuity (0) */
507c4762a1bSJed Brown         uu->r = 1.; /* rho = 1 */
508c4762a1bSJed Brown         uu->E = p0 / (gamma - 1.0);
509c4762a1bSJed Brown       }
5109371c9d4SSatish Balay     } else { /* right of discontinuity (2) */
511de8de29fSMatthew G. Knepley       uu->r = eu->rhoR;
512c4762a1bSJed Brown       uu->E = p0 / (gamma - 1.0);
513c4762a1bSJed Brown     }
5149371c9d4SSatish Balay   } else if (eu->type == EULER_SHOCK_TUBE) {
515c4762a1bSJed Brown     /* For (x<x0) set (rho,u,p)=(8,0,10) and for (x>x0) set (rho,u,p)=(1,0,1). Choose x0 to the midpoint of the domain in the x-direction. */
516c4762a1bSJed Brown     if (x[0] < 0.0) {
517c4762a1bSJed Brown       uu->r = 8.;
518c4762a1bSJed Brown       uu->E = 10. / (gamma - 1.);
5199371c9d4SSatish Balay     } else {
520c4762a1bSJed Brown       uu->r = 1.;
521c4762a1bSJed Brown       uu->E = 1. / (gamma - 1.);
522c4762a1bSJed Brown     }
5239371c9d4SSatish Balay   } else if (eu->type == EULER_LINEAR_WAVE) {
524c4762a1bSJed Brown     initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
5259371c9d4SSatish Balay   } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type);
526c4762a1bSJed Brown 
527c4762a1bSJed Brown   /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */
528de8de29fSMatthew G. Knepley   PetscCall(SpeedOfSound_PG(gamma, uu, &c));
529c4762a1bSJed Brown   c = (uu->ru[0] / uu->r) + c;
530c4762a1bSJed Brown   if (c > phys->maxspeed) phys->maxspeed = c;
5313ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
532c4762a1bSJed Brown }
533c4762a1bSJed Brown 
534c4762a1bSJed Brown /* PetscReal* => EulerNode* conversion */
535d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
536d71ae5a4SJacob Faibussowitsch {
537c4762a1bSJed Brown   PetscInt         i;
538c4762a1bSJed Brown   const EulerNode *xI   = (const EulerNode *)a_xI;
539c4762a1bSJed Brown   EulerNode       *xG   = (EulerNode *)a_xG;
540c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
541c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
5424d86920dSPierre Jolivet 
543c4762a1bSJed Brown   PetscFunctionBeginUser;
544c4762a1bSJed Brown   xG->r = xI->r;                                     /* ghost cell density - same */
545c4762a1bSJed Brown   xG->E = xI->E;                                     /* ghost cell energy - same */
546c4762a1bSJed Brown   if (n[1] != 0.) {                                  /* top and bottom */
547c4762a1bSJed Brown     xG->ru[0] = xI->ru[0];                           /* copy tang to wall */
548c4762a1bSJed Brown     xG->ru[1] = -xI->ru[1];                          /* reflect perp to t/b wall */
5499371c9d4SSatish Balay   } else {                                           /* sides */
550c4762a1bSJed Brown     for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
551c4762a1bSJed Brown   }
552c4762a1bSJed Brown   if (eu->type == EULER_LINEAR_WAVE) { /* debug */
553c4762a1bSJed Brown #if 0
55463a3b9bcSJacob Faibussowitsch     PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]);
555c4762a1bSJed Brown #endif
556c4762a1bSJed Brown   }
5573ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
558c4762a1bSJed Brown }
559c4762a1bSJed Brown 
560d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
561d71ae5a4SJacob Faibussowitsch {
562c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
563c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
564c4762a1bSJed Brown   const EulerNode *x    = (const EulerNode *)xx;
565c4762a1bSJed Brown   PetscReal        p;
566c4762a1bSJed Brown 
567c4762a1bSJed Brown   PetscFunctionBeginUser;
568c4762a1bSJed Brown   f[eu->monitor.Density]  = x->r;
569c4762a1bSJed Brown   f[eu->monitor.Momentum] = NormDIM(x->ru);
570c4762a1bSJed Brown   f[eu->monitor.Energy]   = x->E;
571c4762a1bSJed Brown   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
572de8de29fSMatthew G. Knepley   PetscCall(Pressure_PG(eu->gamma, x, &p));
573c4762a1bSJed Brown   f[eu->monitor.Pressure] = p;
5743ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
575c4762a1bSJed Brown }
576c4762a1bSJed Brown 
577d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys)
578d71ae5a4SJacob Faibussowitsch {
579c4762a1bSJed Brown   Physics_Euler *eu = (Physics_Euler *)phys->data;
58045480ffeSMatthew G. Knepley   DMLabel        label;
58145480ffeSMatthew G. Knepley 
582362febeeSStefano Zampini   PetscFunctionBeginUser;
5839566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
584c4762a1bSJed Brown   if (eu->type == EULER_LINEAR_WAVE) {
585c4762a1bSJed Brown     const PetscInt wallids[] = {100, 101};
586dd39110bSPierre Jolivet     PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
5879371c9d4SSatish Balay   } else {
588c4762a1bSJed Brown     const PetscInt wallids[] = {100, 101, 200, 300};
589dd39110bSPierre Jolivet     PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
590c4762a1bSJed Brown   }
5913ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
592c4762a1bSJed Brown }
593c4762a1bSJed Brown 
594de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
595de8de29fSMatthew G. Knepley static PetscErrorCode CreateQFunctionContext_Euler(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
596de8de29fSMatthew G. Knepley {
597de8de29fSMatthew G. Knepley   Physics_Euler *in = (Physics_Euler *)phys->data;
598de8de29fSMatthew G. Knepley   Physics_Euler *eu;
599de8de29fSMatthew G. Knepley 
600de8de29fSMatthew G. Knepley   PetscFunctionBeginUser;
601de8de29fSMatthew G. Knepley   PetscCall(PetscCalloc1(1, &eu));
602de8de29fSMatthew G. Knepley 
603de8de29fSMatthew G. Knepley   eu->gamma = in->gamma;
604de8de29fSMatthew G. Knepley 
605de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
606de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*eu), eu));
607de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
608de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gamma", offsetof(Physics_Euler, gamma), 1, "Heat capacity ratio"));
609da0802e2SStefano Zampini   PetscFunctionReturn(PETSC_SUCCESS);
610de8de29fSMatthew G. Knepley }
611de8de29fSMatthew G. Knepley #endif
612de8de29fSMatthew G. Knepley 
613de8de29fSMatthew G. Knepley static PetscErrorCode SetupCEED_Euler(DM dm, Physics physics)
614de8de29fSMatthew G. Knepley {
615de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
616de8de29fSMatthew G. Knepley   Ceed                 ceed;
617de8de29fSMatthew G. Knepley   CeedQFunctionContext qfCtx;
618de8de29fSMatthew G. Knepley #endif
619de8de29fSMatthew G. Knepley 
620de8de29fSMatthew G. Knepley   PetscFunctionBegin;
621de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
622de8de29fSMatthew G. Knepley   PetscCall(DMGetCeed(dm, &ceed));
623de8de29fSMatthew G. Knepley   PetscCall(CreateQFunctionContext_Euler(physics, ceed, &qfCtx));
624de8de29fSMatthew G. Knepley   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_Euler_Godunov_CEED, PhysicsRiemann_Euler_Godunov_CEED_loc, qfCtx));
625de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
626de8de29fSMatthew G. Knepley #endif
627de8de29fSMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
628de8de29fSMatthew G. Knepley }
629de8de29fSMatthew G. Knepley 
630d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
631d71ae5a4SJacob Faibussowitsch {
632c4762a1bSJed Brown   Physics_Euler *eu;
633c4762a1bSJed Brown 
634c4762a1bSJed Brown   PetscFunctionBeginUser;
635c4762a1bSJed Brown   phys->field_desc = PhysicsFields_Euler;
636c4762a1bSJed Brown   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Euler_Godunov;
6379566063dSJacob Faibussowitsch   PetscCall(PetscNew(&eu));
638c4762a1bSJed Brown   phys->data     = eu;
639c4762a1bSJed Brown   mod->setupbc   = SetUpBC_Euler;
640de8de29fSMatthew G. Knepley   mod->setupCEED = SetupCEED_Euler;
641de8de29fSMatthew G. Knepley 
642de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov", PhysicsRiemann_Euler_Godunov));
643de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
644de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov_ceed", PhysicsRiemann_Euler_Godunov_CEED));
645de8de29fSMatthew G. Knepley #endif
646de8de29fSMatthew G. Knepley 
647d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
648c4762a1bSJed Brown   {
649de8de29fSMatthew G. Knepley     void (*PhysicsRiemann_Euler)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
650c4762a1bSJed Brown     PetscReal alpha;
651c4762a1bSJed Brown     char      type[64]       = "linear_wave";
652de8de29fSMatthew G. Knepley     char      eu_riemann[64] = "godunov";
653c4762a1bSJed Brown     PetscBool is;
654de8de29fSMatthew G. Knepley     eu->gamma = 1.4;
655de8de29fSMatthew G. Knepley     eu->amach = 2.02;
656de8de29fSMatthew G. Knepley     eu->rhoR  = 3.0;
657de8de29fSMatthew G. Knepley     eu->itana = 0.57735026918963; /* angle of Euler self similar (SS) shock */
658de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsFList("-eu_riemann", "Riemann solver", "", PhysicsRiemannList_Euler, eu_riemann, eu_riemann, sizeof eu_riemann, NULL));
659de8de29fSMatthew G. Knepley     PetscCall(PetscFunctionListFind(PhysicsRiemannList_Euler, eu_riemann, &PhysicsRiemann_Euler));
660de8de29fSMatthew G. Knepley     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Euler;
661de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->gamma, &eu->gamma, NULL));
662de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->amach, &eu->amach, NULL));
663de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->rhoR, &eu->rhoR, NULL));
664c4762a1bSJed Brown     alpha = 60.;
66552ce0ab5SPierre Jolivet     PetscCall(PetscOptionsRangeReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL, 0.0, 90.0));
666de8de29fSMatthew G. Knepley     eu->itana = 1. / PetscTanReal(alpha * PETSC_PI / 180.0);
6679566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL));
6689566063dSJacob Faibussowitsch     PetscCall(PetscStrcmp(type, "linear_wave", &is));
669c4762a1bSJed Brown     if (is) {
67030602db0SMatthew G. Knepley       /* Remember this should be periodic */
671c4762a1bSJed Brown       eu->type = EULER_LINEAR_WAVE;
6729566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave"));
6739371c9d4SSatish Balay     } else {
67454c59aa7SJacob Faibussowitsch       PetscCheck(DIM == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "DIM must be 2 unless linear wave test %s", type);
6759566063dSJacob Faibussowitsch       PetscCall(PetscStrcmp(type, "iv_shock", &is));
676c4762a1bSJed Brown       if (is) {
677c4762a1bSJed Brown         eu->type = EULER_IV_SHOCK;
6789566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock"));
6799371c9d4SSatish Balay       } else {
6809566063dSJacob Faibussowitsch         PetscCall(PetscStrcmp(type, "ss_shock", &is));
681c4762a1bSJed Brown         if (is) {
682c4762a1bSJed Brown           eu->type = EULER_SS_SHOCK;
6839566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock"));
6849371c9d4SSatish Balay         } else {
6859566063dSJacob Faibussowitsch           PetscCall(PetscStrcmp(type, "shock_tube", &is));
686c4762a1bSJed Brown           if (is) eu->type = EULER_SHOCK_TUBE;
68798921bdaSJacob Faibussowitsch           else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type);
6889566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube"));
689c4762a1bSJed Brown         }
690c4762a1bSJed Brown       }
691c4762a1bSJed Brown     }
692c4762a1bSJed Brown   }
693d0609cedSBarry Smith   PetscOptionsHeadEnd();
694c4762a1bSJed Brown   phys->maxspeed = 0.; /* will get set in solution */
6959566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
6969566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
6979566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
6989566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
6999566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
7009566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
7013ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
702c4762a1bSJed Brown }
703c4762a1bSJed Brown 
704d71ae5a4SJacob Faibussowitsch static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx)
705d71ae5a4SJacob Faibussowitsch {
706c4762a1bSJed Brown   PetscReal err = 0.;
707c4762a1bSJed Brown   PetscInt  i, j;
708c4762a1bSJed Brown 
709c4762a1bSJed Brown   PetscFunctionBeginUser;
710c4762a1bSJed Brown   for (i = 0; i < numComps; i++) {
711ad540459SPierre Jolivet     for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
712c4762a1bSJed Brown   }
713c4762a1bSJed Brown   *error = volume * err;
7143ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
715c4762a1bSJed Brown }
716c4762a1bSJed Brown 
717d71ae5a4SJacob Faibussowitsch PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
718d71ae5a4SJacob Faibussowitsch {
719c4762a1bSJed Brown   PetscSF      sfPoint;
720c4762a1bSJed Brown   PetscSection coordSection;
721c4762a1bSJed Brown   Vec          coordinates;
722c4762a1bSJed Brown   PetscSection sectionCell;
723c4762a1bSJed Brown   PetscScalar *part;
724c4762a1bSJed Brown   PetscInt     cStart, cEnd, c;
725c4762a1bSJed Brown   PetscMPIInt  rank;
726c4762a1bSJed Brown 
727c4762a1bSJed Brown   PetscFunctionBeginUser;
7289566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7299566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7309566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, dmCell));
7319566063dSJacob Faibussowitsch   PetscCall(DMGetPointSF(dm, &sfPoint));
7329566063dSJacob Faibussowitsch   PetscCall(DMSetPointSF(*dmCell, sfPoint));
7339566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
7349566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
7359566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
7369566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
7379566063dSJacob Faibussowitsch   PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
7389566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
73948a46eb9SPierre Jolivet   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
7409566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionCell));
7419566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(*dmCell, sectionCell));
7429566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionCell));
7439566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(*dmCell, partition));
7449566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
7459566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*partition, &part));
746c4762a1bSJed Brown   for (c = cStart; c < cEnd; ++c) {
747c4762a1bSJed Brown     PetscScalar *p;
748c4762a1bSJed Brown 
7499566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
750c4762a1bSJed Brown     p[0] = rank;
751c4762a1bSJed Brown   }
7529566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*partition, &part));
7533ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
754c4762a1bSJed Brown }
755c4762a1bSJed Brown 
756d71ae5a4SJacob Faibussowitsch PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
757d71ae5a4SJacob Faibussowitsch {
7583e9753d6SMatthew G. Knepley   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
759c4762a1bSJed Brown   PetscSection       coordSection;
760c4762a1bSJed Brown   Vec                coordinates, facegeom, cellgeom;
761c4762a1bSJed Brown   PetscSection       sectionMass;
762c4762a1bSJed Brown   PetscScalar       *m;
763c4762a1bSJed Brown   const PetscScalar *fgeom, *cgeom, *coords;
764c4762a1bSJed Brown   PetscInt           vStart, vEnd, v;
765c4762a1bSJed Brown 
766c4762a1bSJed Brown   PetscFunctionBeginUser;
7679566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
7689566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7699566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7709566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, &dmMass));
7719566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
7729566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
7739566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
7749566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
7759566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
776c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
777c4762a1bSJed Brown     PetscInt numFaces;
778c4762a1bSJed Brown 
7799566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
7809566063dSJacob Faibussowitsch     PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
781c4762a1bSJed Brown   }
7829566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionMass));
7839566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(dmMass, sectionMass));
7849566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionMass));
7859566063dSJacob Faibussowitsch   PetscCall(DMGetLocalVector(dmMass, massMatrix));
7869566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*massMatrix, &m));
7879566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
7889566063dSJacob Faibussowitsch   PetscCall(VecGetDM(facegeom, &dmFace));
7899566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(facegeom, &fgeom));
7909566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellgeom, &dmCell));
7919566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
7929566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
7939566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(coordinates, &coords));
794c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
795c4762a1bSJed Brown     const PetscInt  *faces;
796c4762a1bSJed Brown     PetscFVFaceGeom *fgA, *fgB, *cg;
797c4762a1bSJed Brown     PetscScalar     *vertex;
798c4762a1bSJed Brown     PetscInt         numFaces, sides[2], f, g;
799c4762a1bSJed Brown 
8009566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
8019566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
8029566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupport(dmMass, v, &faces));
803c4762a1bSJed Brown     for (f = 0; f < numFaces; ++f) {
804c4762a1bSJed Brown       sides[0] = faces[f];
8059566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
806c4762a1bSJed Brown       for (g = 0; g < numFaces; ++g) {
807c4762a1bSJed Brown         const PetscInt *cells = NULL;
808c4762a1bSJed Brown         PetscReal       area  = 0.0;
809c4762a1bSJed Brown         PetscInt        numCells;
810c4762a1bSJed Brown 
811c4762a1bSJed Brown         sides[1] = faces[g];
8129566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
8139566063dSJacob Faibussowitsch         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
81454c59aa7SJacob Faibussowitsch         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
8159566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
816c4762a1bSJed Brown         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
817c4762a1bSJed Brown         area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
818c4762a1bSJed Brown         m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
8199566063dSJacob Faibussowitsch         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
820c4762a1bSJed Brown       }
821c4762a1bSJed Brown     }
822c4762a1bSJed Brown   }
8239566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
8249566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
8259566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(coordinates, &coords));
8269566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*massMatrix, &m));
8279566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmMass));
8289566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
8293ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
830c4762a1bSJed Brown }
831c4762a1bSJed Brown 
832c4762a1bSJed Brown /* Behavior will be different for multi-physics or when using non-default boundary conditions */
833d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
834d71ae5a4SJacob Faibussowitsch {
835c4762a1bSJed Brown   PetscFunctionBeginUser;
836c4762a1bSJed Brown   mod->solution    = func;
837c4762a1bSJed Brown   mod->solutionctx = ctx;
8383ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
839c4762a1bSJed Brown }
840c4762a1bSJed Brown 
841d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
842d71ae5a4SJacob Faibussowitsch {
843c4762a1bSJed Brown   FunctionalLink link, *ptr;
844c4762a1bSJed Brown   PetscInt       lastoffset = -1;
845c4762a1bSJed Brown 
846c4762a1bSJed Brown   PetscFunctionBeginUser;
847c4762a1bSJed Brown   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
8489566063dSJacob Faibussowitsch   PetscCall(PetscNew(&link));
8499566063dSJacob Faibussowitsch   PetscCall(PetscStrallocpy(name, &link->name));
850c4762a1bSJed Brown   link->offset = lastoffset + 1;
851c4762a1bSJed Brown   link->func   = func;
852c4762a1bSJed Brown   link->ctx    = ctx;
853c4762a1bSJed Brown   link->next   = NULL;
854c4762a1bSJed Brown   *ptr         = link;
855c4762a1bSJed Brown   *offset      = link->offset;
8563ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
857c4762a1bSJed Brown }
858c4762a1bSJed Brown 
859d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject)
860d71ae5a4SJacob Faibussowitsch {
861c4762a1bSJed Brown   PetscInt       i, j;
862c4762a1bSJed Brown   FunctionalLink link;
863c4762a1bSJed Brown   char          *names[256];
864c4762a1bSJed Brown 
865c4762a1bSJed Brown   PetscFunctionBeginUser;
866dd39110bSPierre Jolivet   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
8679566063dSJacob Faibussowitsch   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
868c4762a1bSJed Brown   /* Create list of functionals that will be computed somehow */
8699566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
870c4762a1bSJed Brown   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
8719566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
872c4762a1bSJed Brown   mod->numCall = 0;
873c4762a1bSJed Brown   for (i = 0; i < mod->numMonitored; i++) {
874c4762a1bSJed Brown     for (link = mod->functionalRegistry; link; link = link->next) {
875c4762a1bSJed Brown       PetscBool match;
8769566063dSJacob Faibussowitsch       PetscCall(PetscStrcasecmp(names[i], link->name, &match));
877c4762a1bSJed Brown       if (match) break;
878c4762a1bSJed Brown     }
87954c59aa7SJacob Faibussowitsch     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
880c4762a1bSJed Brown     mod->functionalMonitored[i] = link;
881c4762a1bSJed Brown     for (j = 0; j < i; j++) {
882c4762a1bSJed Brown       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
883c4762a1bSJed Brown     }
884c4762a1bSJed Brown     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
885c4762a1bSJed Brown   next_name:
8869566063dSJacob Faibussowitsch     PetscCall(PetscFree(names[i]));
887c4762a1bSJed Brown   }
888c4762a1bSJed Brown 
889c4762a1bSJed Brown   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
890c4762a1bSJed Brown   mod->maxComputed = -1;
891c4762a1bSJed Brown   for (link = mod->functionalRegistry; link; link = link->next) {
892c4762a1bSJed Brown     for (i = 0; i < mod->numCall; i++) {
893c4762a1bSJed Brown       FunctionalLink call = mod->functionalCall[i];
894ad540459SPierre Jolivet       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
895c4762a1bSJed Brown     }
896c4762a1bSJed Brown   }
8973ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
898c4762a1bSJed Brown }
899c4762a1bSJed Brown 
900d71ae5a4SJacob Faibussowitsch static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
901d71ae5a4SJacob Faibussowitsch {
902c4762a1bSJed Brown   FunctionalLink l, next;
903c4762a1bSJed Brown 
904c4762a1bSJed Brown   PetscFunctionBeginUser;
9053ba16761SJacob Faibussowitsch   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
906c4762a1bSJed Brown   l     = *link;
907c4762a1bSJed Brown   *link = NULL;
908c4762a1bSJed Brown   for (; l; l = next) {
909c4762a1bSJed Brown     next = l->next;
9109566063dSJacob Faibussowitsch     PetscCall(PetscFree(l->name));
9119566063dSJacob Faibussowitsch     PetscCall(PetscFree(l));
912c4762a1bSJed Brown   }
9133ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
914c4762a1bSJed Brown }
915c4762a1bSJed Brown 
916c4762a1bSJed Brown /* put the solution callback into a functional callback */
917d71ae5a4SJacob Faibussowitsch static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
918d71ae5a4SJacob Faibussowitsch {
919c4762a1bSJed Brown   Model mod;
9204d86920dSPierre Jolivet 
9217510d9b0SBarry Smith   PetscFunctionBeginUser;
922c4762a1bSJed Brown   mod = (Model)modctx;
9239566063dSJacob Faibussowitsch   PetscCall((*mod->solution)(mod, time, x, u, mod->solutionctx));
9243ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
925c4762a1bSJed Brown }
926c4762a1bSJed Brown 
927d71ae5a4SJacob Faibussowitsch PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
928d71ae5a4SJacob Faibussowitsch {
929c4762a1bSJed Brown   PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
930c4762a1bSJed Brown   void *ctx[1];
931c4762a1bSJed Brown   Model mod = user->model;
932c4762a1bSJed Brown 
933c4762a1bSJed Brown   PetscFunctionBeginUser;
934c4762a1bSJed Brown   func[0] = SolutionFunctional;
935c4762a1bSJed Brown   ctx[0]  = (void *)mod;
9369566063dSJacob Faibussowitsch   PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X));
9373ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
938c4762a1bSJed Brown }
939c4762a1bSJed Brown 
940d71ae5a4SJacob Faibussowitsch static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
941d71ae5a4SJacob Faibussowitsch {
942c4762a1bSJed Brown   PetscFunctionBeginUser;
9439566063dSJacob Faibussowitsch   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
9449566063dSJacob Faibussowitsch   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
9459566063dSJacob Faibussowitsch   PetscCall(PetscViewerFileSetName(*viewer, filename));
9463ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
947c4762a1bSJed Brown }
948c4762a1bSJed Brown 
949d71ae5a4SJacob Faibussowitsch static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
950d71ae5a4SJacob Faibussowitsch {
951c4762a1bSJed Brown   User        user = (User)ctx;
9523e9753d6SMatthew G. Knepley   DM          dm, plex;
953c4762a1bSJed Brown   PetscViewer viewer;
954c4762a1bSJed Brown   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
955c4762a1bSJed Brown   PetscReal   xnorm;
956ecc87898SStefano Zampini   PetscBool   rollback;
957c4762a1bSJed Brown 
958c4762a1bSJed Brown   PetscFunctionBeginUser;
959ecc87898SStefano Zampini   PetscCall(TSGetStepRollBack(ts, &rollback));
960ecc87898SStefano Zampini   if (rollback) PetscFunctionReturn(PETSC_SUCCESS);
9619566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "u"));
9629566063dSJacob Faibussowitsch   PetscCall(VecGetDM(X, &dm));
9639566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));
964c4762a1bSJed Brown 
965ad540459SPierre Jolivet   if (stepnum >= 0) stepnum += user->monitorStepOffset;
966c4762a1bSJed Brown   if (stepnum >= 0) { /* No summary for final time */
967c4762a1bSJed Brown     Model              mod = user->model;
9683e9753d6SMatthew G. Knepley     Vec                cellgeom;
969c4762a1bSJed Brown     PetscInt           c, cStart, cEnd, fcount, i;
970c4762a1bSJed Brown     size_t             ftableused, ftablealloc;
971c4762a1bSJed Brown     const PetscScalar *cgeom, *x;
972c4762a1bSJed Brown     DM                 dmCell;
973c4762a1bSJed Brown     DMLabel            vtkLabel;
974c4762a1bSJed Brown     PetscReal         *fmin, *fmax, *fintegral, *ftmp;
9753e9753d6SMatthew G. Knepley 
9769566063dSJacob Faibussowitsch     PetscCall(DMConvert(dm, DMPLEX, &plex));
9779566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
978c4762a1bSJed Brown     fcount = mod->maxComputed + 1;
9799566063dSJacob Faibussowitsch     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
980c4762a1bSJed Brown     for (i = 0; i < fcount; i++) {
981c4762a1bSJed Brown       fmin[i]      = PETSC_MAX_REAL;
982c4762a1bSJed Brown       fmax[i]      = PETSC_MIN_REAL;
983c4762a1bSJed Brown       fintegral[i] = 0;
984c4762a1bSJed Brown     }
9859566063dSJacob Faibussowitsch     PetscCall(VecGetDM(cellgeom, &dmCell));
9869566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd));
9879566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
9889566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(X, &x));
9899566063dSJacob Faibussowitsch     PetscCall(DMGetLabel(dm, "vtk", &vtkLabel));
990c4762a1bSJed Brown     for (c = cStart; c < cEnd; ++c) {
991c4762a1bSJed Brown       PetscFVCellGeom   *cg;
992c4762a1bSJed Brown       const PetscScalar *cx     = NULL;
993c4762a1bSJed Brown       PetscInt           vtkVal = 0;
994c4762a1bSJed Brown 
995c4762a1bSJed Brown       /* not that these two routines as currently implemented work for any dm with a
996c4762a1bSJed Brown        * localSection/globalSection */
9979566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
9989566063dSJacob Faibussowitsch       PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
9999566063dSJacob Faibussowitsch       if (vtkLabel) PetscCall(DMLabelGetValue(vtkLabel, c, &vtkVal));
1000c4762a1bSJed Brown       if (!vtkVal || !cx) continue; /* ghost, or not a global cell */
1001c4762a1bSJed Brown       for (i = 0; i < mod->numCall; i++) {
1002c4762a1bSJed Brown         FunctionalLink flink = mod->functionalCall[i];
10039566063dSJacob Faibussowitsch         PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
1004c4762a1bSJed Brown       }
1005c4762a1bSJed Brown       for (i = 0; i < fcount; i++) {
1006c4762a1bSJed Brown         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1007c4762a1bSJed Brown         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1008c4762a1bSJed Brown         fintegral[i] += cg->volume * ftmp[i];
1009c4762a1bSJed Brown       }
1010c4762a1bSJed Brown     }
10119566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
10129566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(X, &x));
10139566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&plex));
1014712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1015712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1016712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
1017c4762a1bSJed Brown 
1018c4762a1bSJed Brown     ftablealloc = fcount * 100;
1019c4762a1bSJed Brown     ftableused  = 0;
10209566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(ftablealloc, &ftable));
1021c4762a1bSJed Brown     for (i = 0; i < mod->numMonitored; i++) {
1022c4762a1bSJed Brown       size_t         countused;
1023c4762a1bSJed Brown       char           buffer[256], *p;
1024c4762a1bSJed Brown       FunctionalLink flink = mod->functionalMonitored[i];
1025c4762a1bSJed Brown       PetscInt       id    = flink->offset;
1026c4762a1bSJed Brown       if (i % 3) {
10279566063dSJacob Faibussowitsch         PetscCall(PetscArraycpy(buffer, "  ", 2));
1028c4762a1bSJed Brown         p = buffer + 2;
1029c4762a1bSJed Brown       } else if (i) {
1030c4762a1bSJed Brown         char newline[] = "\n";
10319566063dSJacob Faibussowitsch         PetscCall(PetscMemcpy(buffer, newline, sizeof(newline) - 1));
1032c4762a1bSJed Brown         p = buffer + sizeof(newline) - 1;
1033c4762a1bSJed Brown       } else {
1034c4762a1bSJed Brown         p = buffer;
1035c4762a1bSJed Brown       }
10369566063dSJacob 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]));
1037c4762a1bSJed Brown       countused--;
1038c4762a1bSJed Brown       countused += p - buffer;
1039c4762a1bSJed Brown       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1040c4762a1bSJed Brown         char *ftablenew;
1041c4762a1bSJed Brown         ftablealloc = 2 * ftablealloc + countused;
10429566063dSJacob Faibussowitsch         PetscCall(PetscMalloc(ftablealloc, &ftablenew));
10439566063dSJacob Faibussowitsch         PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
10449566063dSJacob Faibussowitsch         PetscCall(PetscFree(ftable));
1045c4762a1bSJed Brown         ftable = ftablenew;
1046c4762a1bSJed Brown       }
10479566063dSJacob Faibussowitsch       PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
1048c4762a1bSJed Brown       ftableused += countused;
1049c4762a1bSJed Brown       ftable[ftableused] = 0;
1050c4762a1bSJed Brown     }
10519566063dSJacob Faibussowitsch     PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));
1052c4762a1bSJed Brown 
105363a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
10549566063dSJacob Faibussowitsch     PetscCall(PetscFree(ftable));
1055c4762a1bSJed Brown   }
10563ba16761SJacob Faibussowitsch   if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
1057c4762a1bSJed Brown   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1058c4762a1bSJed Brown     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
10599566063dSJacob Faibussowitsch       PetscCall(TSGetStepNumber(ts, &stepnum));
1060c4762a1bSJed Brown     }
106163a3b9bcSJacob Faibussowitsch     PetscCall(PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum));
10629566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, filename, &viewer));
10639566063dSJacob Faibussowitsch     PetscCall(VecView(X, viewer));
10649566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
1065c4762a1bSJed Brown   }
10663ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1067c4762a1bSJed Brown }
1068c4762a1bSJed Brown 
1069d71ae5a4SJacob Faibussowitsch static PetscErrorCode initializeTS(DM dm, User user, TS *ts)
1070d71ae5a4SJacob Faibussowitsch {
1071de8de29fSMatthew G. Knepley   PetscBool useCeed;
1072de8de29fSMatthew G. Knepley 
10737510d9b0SBarry Smith   PetscFunctionBeginUser;
10749566063dSJacob Faibussowitsch   PetscCall(TSCreate(PetscObjectComm((PetscObject)dm), ts));
10759566063dSJacob Faibussowitsch   PetscCall(TSSetType(*ts, TSSSP));
10769566063dSJacob Faibussowitsch   PetscCall(TSSetDM(*ts, dm));
10771baa6e33SBarry Smith   if (user->vtkmon) PetscCall(TSMonitorSet(*ts, MonitorVTK, user, NULL));
1078de8de29fSMatthew G. Knepley   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
10799566063dSJacob Faibussowitsch   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user));
1080de8de29fSMatthew G. Knepley   if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user));
1081de8de29fSMatthew G. Knepley   else PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
10829566063dSJacob Faibussowitsch   PetscCall(TSSetMaxTime(*ts, 2.0));
10839566063dSJacob Faibussowitsch   PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
10843ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1085c4762a1bSJed Brown }
1086c4762a1bSJed Brown 
108752115386SStefano Zampini typedef struct {
108852115386SStefano Zampini   PetscFV      fvm;
108952115386SStefano Zampini   VecTagger    refineTag;
109052115386SStefano Zampini   VecTagger    coarsenTag;
109152115386SStefano Zampini   DM           adaptedDM;
109252115386SStefano Zampini   User         user;
109352115386SStefano Zampini   PetscReal    cfl;
109452115386SStefano Zampini   PetscLimiter limiter;
109552115386SStefano Zampini   PetscLimiter noneLimiter;
109652115386SStefano Zampini } TransferCtx;
109752115386SStefano Zampini 
109852115386SStefano Zampini static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
1099d71ae5a4SJacob Faibussowitsch {
110052115386SStefano Zampini   TransferCtx       *tctx       = (TransferCtx *)ctx;
110152115386SStefano Zampini   PetscFV            fvm        = tctx->fvm;
110252115386SStefano Zampini   VecTagger          refineTag  = tctx->refineTag;
110352115386SStefano Zampini   VecTagger          coarsenTag = tctx->coarsenTag;
110452115386SStefano Zampini   User               user       = tctx->user;
1105c4762a1bSJed Brown   DM                 dm, gradDM, plex, cellDM, adaptedDM = NULL;
1106c4762a1bSJed Brown   Vec                cellGeom, faceGeom;
1107ecc87898SStefano Zampini   PetscBool          computeGradient;
1108c4762a1bSJed Brown   Vec                grad, locGrad, locX, errVec;
1109c4762a1bSJed Brown   PetscInt           cStart, cEnd, c, dim, nRefine, nCoarsen;
111052115386SStefano Zampini   PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1111c4762a1bSJed Brown   PetscScalar       *errArray;
1112c4762a1bSJed Brown   const PetscScalar *pointVals;
1113c4762a1bSJed Brown   const PetscScalar *pointGrads;
1114c4762a1bSJed Brown   const PetscScalar *pointGeom;
1115c4762a1bSJed Brown   DMLabel            adaptLabel = NULL;
1116c4762a1bSJed Brown   IS                 refineIS, coarsenIS;
1117c4762a1bSJed Brown 
11187510d9b0SBarry Smith   PetscFunctionBeginUser;
111952115386SStefano Zampini   *resize = PETSC_FALSE;
11209566063dSJacob Faibussowitsch   PetscCall(VecGetDM(sol, &dm));
11219566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
112252115386SStefano Zampini   PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
11239566063dSJacob Faibussowitsch   PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
11249566063dSJacob Faibussowitsch   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
11259566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
11269566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
11279566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(plex, &locX));
11289566063dSJacob Faibussowitsch   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
11299566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
11309566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
11319566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(gradDM, &grad));
11329566063dSJacob Faibussowitsch   PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
11339566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(gradDM, &locGrad));
11349566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
11359566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
11369566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&grad));
11379566063dSJacob Faibussowitsch   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
11389566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(locGrad, &pointGrads));
11399566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
11409566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(locX, &pointVals));
11419566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellGeom, &cellDM));
11429566063dSJacob Faibussowitsch   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
114377433607SBarry Smith   PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
11449566063dSJacob Faibussowitsch   PetscCall(VecSetUp(errVec));
11459566063dSJacob Faibussowitsch   PetscCall(VecGetArray(errVec, &errArray));
1146c4762a1bSJed Brown   for (c = cStart; c < cEnd; c++) {
1147c4762a1bSJed Brown     PetscReal        errInd = 0.;
1148c4762a1bSJed Brown     PetscScalar     *pointGrad;
1149c4762a1bSJed Brown     PetscScalar     *pointVal;
1150c4762a1bSJed Brown     PetscFVCellGeom *cg;
1151c4762a1bSJed Brown 
11529566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
11539566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
11549566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));
1155c4762a1bSJed Brown 
1156ecc87898SStefano Zampini     PetscCall((*user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1157c4762a1bSJed Brown     errArray[c - cStart] = errInd;
1158c4762a1bSJed Brown     minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
1159c4762a1bSJed Brown     minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
1160c4762a1bSJed Brown   }
11619566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(errVec, &errArray));
11629566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(locX, &pointVals));
11639566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
11649566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
11659566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&locGrad));
11669566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&locX));
11679566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1168c4762a1bSJed Brown 
11699566063dSJacob Faibussowitsch   PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
11709566063dSJacob Faibussowitsch   PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
11719566063dSJacob Faibussowitsch   PetscCall(ISGetSize(refineIS, &nRefine));
11729566063dSJacob Faibussowitsch   PetscCall(ISGetSize(coarsenIS, &nCoarsen));
11739566063dSJacob Faibussowitsch   if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
11749566063dSJacob Faibussowitsch   if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
11759566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&coarsenIS));
11769566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&refineIS));
11779566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&errVec));
1178c4762a1bSJed Brown 
11799566063dSJacob Faibussowitsch   PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
118052115386SStefano Zampini   PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1181c4762a1bSJed Brown   minMaxInd[1] = -minMaxInd[1];
1182712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
118352115386SStefano Zampini   PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1184c4762a1bSJed Brown   if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
11859566063dSJacob Faibussowitsch     PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1186c4762a1bSJed Brown   }
11879566063dSJacob Faibussowitsch   PetscCall(DMLabelDestroy(&adaptLabel));
1188c4762a1bSJed Brown   if (adaptedDM) {
118952115386SStefano Zampini     tctx->adaptedDM = adaptedDM;
119052115386SStefano 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));
119152115386SStefano Zampini     *resize = PETSC_TRUE;
1192c4762a1bSJed Brown   } else {
119352115386SStefano Zampini     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1194c4762a1bSJed Brown   }
11953ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1196c4762a1bSJed Brown }
1197c4762a1bSJed Brown 
119852115386SStefano Zampini static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
119952115386SStefano Zampini {
120052115386SStefano Zampini   TransferCtx *tctx = (TransferCtx *)ctx;
120152115386SStefano Zampini   DM           dm;
120252115386SStefano Zampini   PetscReal    time;
120352115386SStefano Zampini 
120452115386SStefano Zampini   PetscFunctionBeginUser;
120552115386SStefano Zampini   PetscCall(TSGetDM(ts, &dm));
120652115386SStefano Zampini   PetscCall(TSGetTime(ts, &time));
120752115386SStefano Zampini   PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
120852115386SStefano Zampini   for (PetscInt i = 0; i < nv; i++) {
120952115386SStefano Zampini     PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
121052115386SStefano Zampini     PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
121152115386SStefano Zampini   }
121252115386SStefano Zampini   PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */
121352115386SStefano Zampini 
121452115386SStefano Zampini   Model     mod  = tctx->user->model;
121552115386SStefano Zampini   Physics   phys = mod->physics;
121652115386SStefano Zampini   PetscReal minRadius;
121752115386SStefano Zampini 
121852115386SStefano Zampini   PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius));
121952115386SStefano Zampini   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
122052115386SStefano Zampini   PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed");
122152115386SStefano Zampini 
122252115386SStefano Zampini   PetscReal dt = tctx->cfl * minRadius / mod->maxspeed;
122352115386SStefano Zampini   PetscCall(TSSetTimeStep(ts, dt));
122452115386SStefano Zampini 
122552115386SStefano Zampini   PetscCall(TSSetDM(ts, tctx->adaptedDM));
122652115386SStefano Zampini   PetscCall(DMDestroy(&tctx->adaptedDM));
122752115386SStefano Zampini   PetscFunctionReturn(PETSC_SUCCESS);
122852115386SStefano Zampini }
122952115386SStefano Zampini 
1230d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
1231d71ae5a4SJacob Faibussowitsch {
1232c4762a1bSJed Brown   MPI_Comm          comm;
1233c4762a1bSJed Brown   PetscDS           prob;
1234c4762a1bSJed Brown   PetscFV           fvm;
1235c4762a1bSJed Brown   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1236c4762a1bSJed Brown   User              user;
1237c4762a1bSJed Brown   Model             mod;
1238c4762a1bSJed Brown   Physics           phys;
12393e9753d6SMatthew G. Knepley   DM                dm, plex;
1240c4762a1bSJed Brown   PetscReal         ftime, cfl, dt, minRadius;
1241c4762a1bSJed Brown   PetscInt          dim, nsteps;
1242c4762a1bSJed Brown   TS                ts;
1243c4762a1bSJed Brown   TSConvergedReason reason;
1244c4762a1bSJed Brown   Vec               X;
1245c4762a1bSJed Brown   PetscViewer       viewer;
1246b5a892a1SMatthew G. Knepley   PetscBool         vtkCellGeom, useAMR;
124730602db0SMatthew G. Knepley   PetscInt          adaptInterval;
1248c4762a1bSJed Brown   char              physname[256] = "advect";
1249c4762a1bSJed Brown   VecTagger         refineTag = NULL, coarsenTag = NULL;
125052115386SStefano Zampini   TransferCtx       tctx;
1251c4762a1bSJed Brown 
1252327415f7SBarry Smith   PetscFunctionBeginUser;
12539566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1254c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
1255c4762a1bSJed Brown 
12569566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user));
12579566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model));
12589566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model->physics));
1259c4762a1bSJed Brown   mod           = user->model;
1260c4762a1bSJed Brown   phys          = mod->physics;
1261c4762a1bSJed Brown   mod->comm     = comm;
1262c4762a1bSJed Brown   useAMR        = PETSC_FALSE;
1263c4762a1bSJed Brown   adaptInterval = 1;
1264c4762a1bSJed Brown 
1265c4762a1bSJed Brown   /* Register physical models to be available on the command line */
12669566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
12679566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
12689566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1269c4762a1bSJed Brown 
1270d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1271c4762a1bSJed Brown   {
1272c4762a1bSJed Brown     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
12739566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1274c4762a1bSJed Brown     user->vtkInterval = 1;
12759566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1276c4762a1bSJed Brown     user->vtkmon = PETSC_TRUE;
12779566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1278c4762a1bSJed Brown     vtkCellGeom = PETSC_FALSE;
1279c6a7a370SJeremy L Thompson     PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
12809566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
12819566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
12829566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
12839566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1284c4762a1bSJed Brown   }
1285d0609cedSBarry Smith   PetscOptionsEnd();
1286c4762a1bSJed Brown 
1287c4762a1bSJed Brown   if (useAMR) {
1288c4762a1bSJed Brown     VecTaggerBox refineBox, coarsenBox;
1289c4762a1bSJed Brown 
1290c4762a1bSJed Brown     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1291c4762a1bSJed Brown     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;
1292c4762a1bSJed Brown 
12939566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &refineTag));
12949566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
12959566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
12969566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
12979566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(refineTag));
12989566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(refineTag));
12999566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
1300c4762a1bSJed Brown 
13019566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &coarsenTag));
13029566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
13039566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
13049566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
13059566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(coarsenTag));
13069566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(coarsenTag));
13079566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1308c4762a1bSJed Brown   }
1309c4762a1bSJed Brown 
1310d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1311c4762a1bSJed Brown   {
1312c4762a1bSJed Brown     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
13139566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
13149566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
13159566063dSJacob Faibussowitsch     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
13169566063dSJacob Faibussowitsch     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1317c4762a1bSJed Brown     /* Count number of fields and dofs */
1318c4762a1bSJed Brown     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
131954c59aa7SJacob Faibussowitsch     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
13209566063dSJacob Faibussowitsch     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1321c4762a1bSJed Brown   }
1322d0609cedSBarry Smith   PetscOptionsEnd();
1323c4762a1bSJed Brown 
1324c4762a1bSJed Brown   /* Create mesh */
1325c4762a1bSJed Brown   {
132630602db0SMatthew G. Knepley     PetscInt i;
132730602db0SMatthew G. Knepley 
13289566063dSJacob Faibussowitsch     PetscCall(DMCreate(comm, &dm));
13299566063dSJacob Faibussowitsch     PetscCall(DMSetType(dm, DMPLEX));
13309566063dSJacob Faibussowitsch     PetscCall(DMSetFromOptions(dm));
13319371c9d4SSatish Balay     for (i = 0; i < DIM; i++) {
13329371c9d4SSatish Balay       mod->bounds[2 * i]     = 0.;
13339371c9d4SSatish Balay       mod->bounds[2 * i + 1] = 1.;
13349371c9d4SSatish Balay     };
1335c4762a1bSJed Brown     dim = DIM;
133630602db0SMatthew G. Knepley     { /* a null name means just do a hex box */
133730602db0SMatthew G. Knepley       PetscInt  cells[3] = {1, 1, 1}, n = 3;
133830602db0SMatthew G. Knepley       PetscBool flg2, skew              = PETSC_FALSE;
1339c4762a1bSJed Brown       PetscInt  nret2 = 2 * DIM;
1340d0609cedSBarry Smith       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
13419566063dSJacob 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));
13429566063dSJacob Faibussowitsch       PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
13439566063dSJacob Faibussowitsch       PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1344d0609cedSBarry Smith       PetscOptionsEnd();
134530602db0SMatthew G. Knepley       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1346c4762a1bSJed Brown       if (flg2) {
1347c4762a1bSJed Brown         PetscInt     dimEmbed, i;
1348c4762a1bSJed Brown         PetscInt     nCoords;
1349c4762a1bSJed Brown         PetscScalar *coords;
1350c4762a1bSJed Brown         Vec          coordinates;
1351c4762a1bSJed Brown 
13529566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
13539566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
13549566063dSJacob Faibussowitsch         PetscCall(VecGetLocalSize(coordinates, &nCoords));
135554c59aa7SJacob Faibussowitsch         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
13569566063dSJacob Faibussowitsch         PetscCall(VecGetArray(coordinates, &coords));
1357c4762a1bSJed Brown         for (i = 0; i < nCoords; i += dimEmbed) {
1358c4762a1bSJed Brown           PetscInt j;
1359c4762a1bSJed Brown 
1360c4762a1bSJed Brown           PetscScalar *coord = &coords[i];
1361c4762a1bSJed Brown           for (j = 0; j < dimEmbed; j++) {
1362c4762a1bSJed Brown             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1363c4762a1bSJed Brown             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1364c4762a1bSJed Brown               if (cells[0] == 2 && i == 8) {
1365c4762a1bSJed Brown                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
136630602db0SMatthew G. Knepley               } else if (cells[0] == 3) {
1367c4762a1bSJed Brown                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1368c4762a1bSJed Brown                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1369c4762a1bSJed Brown                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1370c4762a1bSJed Brown               }
1371c4762a1bSJed Brown             }
1372c4762a1bSJed Brown           }
1373c4762a1bSJed Brown         }
13749566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(coordinates, &coords));
13759566063dSJacob Faibussowitsch         PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1376c4762a1bSJed Brown       }
1377c4762a1bSJed Brown     }
1378c4762a1bSJed Brown   }
13799566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
13809566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
1381c4762a1bSJed Brown 
1382c4762a1bSJed Brown   /* set up BCs, functions, tags */
13839566063dSJacob Faibussowitsch   PetscCall(DMCreateLabel(dm, "Face Sets"));
1384c4762a1bSJed Brown   mod->errorIndicator = ErrorIndicator_Simple;
1385c4762a1bSJed Brown 
1386c4762a1bSJed Brown   {
1387c4762a1bSJed Brown     DM gdm;
1388c4762a1bSJed Brown 
13899566063dSJacob Faibussowitsch     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
13909566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dm));
1391c4762a1bSJed Brown     dm = gdm;
13929566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1393c4762a1bSJed Brown   }
1394c4762a1bSJed Brown 
13959566063dSJacob Faibussowitsch   PetscCall(PetscFVCreate(comm, &fvm));
13969566063dSJacob Faibussowitsch   PetscCall(PetscFVSetFromOptions(fvm));
13979566063dSJacob Faibussowitsch   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
13989566063dSJacob Faibussowitsch   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
13999566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1400c4762a1bSJed Brown   {
1401c4762a1bSJed Brown     PetscInt f, dof;
1402c4762a1bSJed Brown     for (f = 0, dof = 0; f < phys->nfields; f++) {
1403c4762a1bSJed Brown       PetscInt newDof = phys->field_desc[f].dof;
1404c4762a1bSJed Brown 
1405c4762a1bSJed Brown       if (newDof == 1) {
14069566063dSJacob Faibussowitsch         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
14079371c9d4SSatish Balay       } else {
1408c4762a1bSJed Brown         PetscInt j;
1409c4762a1bSJed Brown 
1410c4762a1bSJed Brown         for (j = 0; j < newDof; j++) {
1411c4762a1bSJed Brown           char compName[256] = "Unknown";
1412c4762a1bSJed Brown 
141363a3b9bcSJacob Faibussowitsch           PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
14149566063dSJacob Faibussowitsch           PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1415c4762a1bSJed Brown         }
1416c4762a1bSJed Brown       }
1417c4762a1bSJed Brown       dof += newDof;
1418c4762a1bSJed Brown     }
1419c4762a1bSJed Brown   }
1420c4762a1bSJed Brown   /* FV is now structured with one field having all physics as components */
14219566063dSJacob Faibussowitsch   PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
14229566063dSJacob Faibussowitsch   PetscCall(DMCreateDS(dm));
14239566063dSJacob Faibussowitsch   PetscCall(DMGetDS(dm, &prob));
14249566063dSJacob Faibussowitsch   PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
14259566063dSJacob Faibussowitsch   PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
14269566063dSJacob Faibussowitsch   PetscCall((*mod->setupbc)(dm, prob, phys));
14279566063dSJacob Faibussowitsch   PetscCall(PetscDSSetFromOptions(prob));
1428c4762a1bSJed Brown   {
1429c4762a1bSJed Brown     char      convType[256];
1430c4762a1bSJed Brown     PetscBool flg;
1431c4762a1bSJed Brown 
1432d0609cedSBarry Smith     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
14339566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1434d0609cedSBarry Smith     PetscOptionsEnd();
1435c4762a1bSJed Brown     if (flg) {
1436c4762a1bSJed Brown       DM dmConv;
1437c4762a1bSJed Brown 
14389566063dSJacob Faibussowitsch       PetscCall(DMConvert(dm, convType, &dmConv));
1439c4762a1bSJed Brown       if (dmConv) {
14409566063dSJacob Faibussowitsch         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
14419566063dSJacob Faibussowitsch         PetscCall(DMDestroy(&dm));
1442c4762a1bSJed Brown         dm = dmConv;
14439566063dSJacob Faibussowitsch         PetscCall(DMSetFromOptions(dm));
1444c4762a1bSJed Brown       }
1445c4762a1bSJed Brown     }
1446c4762a1bSJed Brown   }
1447de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
1448de8de29fSMatthew G. Knepley   {
1449de8de29fSMatthew G. Knepley     PetscBool useCeed;
1450de8de29fSMatthew G. Knepley     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1451de8de29fSMatthew G. Knepley     if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1452de8de29fSMatthew G. Knepley   }
1453de8de29fSMatthew G. Knepley #endif
1454c4762a1bSJed Brown 
14559566063dSJacob Faibussowitsch   PetscCall(initializeTS(dm, user, &ts));
1456c4762a1bSJed Brown 
14579566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(dm, &X));
14589566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
14599566063dSJacob Faibussowitsch   PetscCall(SetInitialCondition(dm, X, user));
1460c4762a1bSJed Brown   if (useAMR) {
1461c4762a1bSJed Brown     PetscInt adaptIter;
1462c4762a1bSJed Brown 
1463c4762a1bSJed Brown     /* use no limiting when reconstructing gradients for adaptivity */
14649566063dSJacob Faibussowitsch     PetscCall(PetscFVGetLimiter(fvm, &limiter));
14659566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)limiter));
14669566063dSJacob Faibussowitsch     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
14679566063dSJacob Faibussowitsch     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
1468c4762a1bSJed Brown 
146952115386SStefano Zampini     /* Refinement context */
147052115386SStefano Zampini     tctx.fvm         = fvm;
147152115386SStefano Zampini     tctx.refineTag   = refineTag;
147252115386SStefano Zampini     tctx.coarsenTag  = coarsenTag;
147352115386SStefano Zampini     tctx.adaptedDM   = NULL;
147452115386SStefano Zampini     tctx.user        = user;
147552115386SStefano Zampini     tctx.noneLimiter = noneLimiter;
147652115386SStefano Zampini     tctx.limiter     = limiter;
147752115386SStefano Zampini     tctx.cfl         = cfl;
147852115386SStefano Zampini 
147952115386SStefano Zampini     /* Do some initial refinement steps */
1480c4762a1bSJed Brown     for (adaptIter = 0;; ++adaptIter) {
1481c4762a1bSJed Brown       PetscLogDouble bytes;
148252115386SStefano Zampini       PetscBool      resize;
1483c4762a1bSJed Brown 
14849566063dSJacob Faibussowitsch       PetscCall(PetscMemoryGetCurrentUsage(&bytes));
148563a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
14869566063dSJacob Faibussowitsch       PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
14879566063dSJacob Faibussowitsch       PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));
1488c4762a1bSJed Brown 
148952115386SStefano Zampini       PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
149052115386SStefano Zampini       if (!resize) break;
14919566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dm));
14929566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&X));
149352115386SStefano Zampini       dm             = tctx.adaptedDM;
149452115386SStefano Zampini       tctx.adaptedDM = NULL;
149552115386SStefano Zampini       PetscCall(TSSetDM(ts, dm));
14969566063dSJacob Faibussowitsch       PetscCall(DMCreateGlobalVector(dm, &X));
14979566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
14989566063dSJacob Faibussowitsch       PetscCall(SetInitialCondition(dm, X, user));
1499c4762a1bSJed Brown     }
1500c4762a1bSJed Brown   }
1501c4762a1bSJed Brown 
15029566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
1503c4762a1bSJed Brown   if (vtkCellGeom) {
1504c4762a1bSJed Brown     DM  dmCell;
1505c4762a1bSJed Brown     Vec cellgeom, partition;
1506c4762a1bSJed Brown 
15079566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
15089566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
15099566063dSJacob Faibussowitsch     PetscCall(VecView(cellgeom, viewer));
15109566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15119566063dSJacob Faibussowitsch     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
15129566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
15139566063dSJacob Faibussowitsch     PetscCall(VecView(partition, viewer));
15149566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15159566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&partition));
15169566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dmCell));
1517c4762a1bSJed Brown   }
1518c4762a1bSJed Brown   /* collect max maxspeed from all processes -- todo */
15199566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
15209566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1521712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
152254c59aa7SJacob Faibussowitsch   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1523c4762a1bSJed Brown   dt = cfl * minRadius / mod->maxspeed;
15249566063dSJacob Faibussowitsch   PetscCall(TSSetTimeStep(ts, dt));
15259566063dSJacob Faibussowitsch   PetscCall(TSSetFromOptions(ts));
1526c4762a1bSJed Brown 
152752115386SStefano Zampini   /* When using adaptive mesh refinement
152852115386SStefano Zampini      specify callbacks to refine the solution
1529ecc87898SStefano Zampini      and interpolate data from old to new mesh
1530ecc87898SStefano Zampini      When mesh adaption is requested, the step will be restarted
1531ecc87898SStefano Zampini   */
1532ecc87898SStefano Zampini   if (useAMR) PetscCall(TSSetResize(ts, PETSC_TRUE, adaptToleranceFVMSetUp, Transfer, &tctx));
153352115386SStefano Zampini   PetscCall(TSSetSolution(ts, X));
15349566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&X));
153552115386SStefano Zampini   PetscCall(TSSolve(ts, NULL));
15369566063dSJacob Faibussowitsch   PetscCall(TSGetSolveTime(ts, &ftime));
15379566063dSJacob Faibussowitsch   PetscCall(TSGetStepNumber(ts, &nsteps));
15389566063dSJacob Faibussowitsch   PetscCall(TSGetConvergedReason(ts, &reason));
153963a3b9bcSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
15409566063dSJacob Faibussowitsch   PetscCall(TSDestroy(&ts));
1541c4762a1bSJed Brown 
15429566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&refineTag));
15439566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&coarsenTag));
15449566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsList));
15459566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1546de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
15479566063dSJacob Faibussowitsch   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
15489566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalMonitored));
15499566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalCall));
15509566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics->data));
15519566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics));
15529566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model));
15539566063dSJacob Faibussowitsch   PetscCall(PetscFree(user));
15549566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&limiter));
15559566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&noneLimiter));
15569566063dSJacob Faibussowitsch   PetscCall(PetscFVDestroy(&fvm));
15579566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dm));
15589566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
1559b122ec5aSJacob Faibussowitsch   return 0;
1560c4762a1bSJed Brown }
1561c4762a1bSJed Brown 
1562c4762a1bSJed Brown /* Subroutine to set up the initial conditions for the */
1563c4762a1bSJed Brown /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1564c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1565d71ae5a4SJacob Faibussowitsch int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1566d71ae5a4SJacob Faibussowitsch {
1567c4762a1bSJed Brown   int j, k;
1568c4762a1bSJed Brown   /*      Wc=matmul(lv,Ueq) 3 vars */
1569c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1570c4762a1bSJed Brown     wc[k] = 0.;
1571ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1572c4762a1bSJed Brown   }
1573c4762a1bSJed Brown   return 0;
1574c4762a1bSJed Brown }
1575c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1576d71ae5a4SJacob Faibussowitsch int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1577d71ae5a4SJacob Faibussowitsch {
1578c4762a1bSJed Brown   int k, j;
1579c4762a1bSJed Brown   /*      V=matmul(rv,WC) 3 vars */
1580c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1581c4762a1bSJed Brown     v[k] = 0.;
1582ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1583c4762a1bSJed Brown   }
1584c4762a1bSJed Brown   return 0;
1585c4762a1bSJed Brown }
1586c4762a1bSJed Brown /* ---------------------------------------------------------------------- */
1587d71ae5a4SJacob Faibussowitsch int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1588d71ae5a4SJacob Faibussowitsch {
1589c4762a1bSJed Brown   int       j, k;
1590c4762a1bSJed Brown   PetscReal rho, csnd, p0;
1591c4762a1bSJed Brown   /* PetscScalar u; */
1592c4762a1bSJed Brown 
15939371c9d4SSatish Balay   for (k = 0; k < 3; ++k)
15949371c9d4SSatish Balay     for (j = 0; j < 3; ++j) {
15959371c9d4SSatish Balay       lv[k][j] = 0.;
15969371c9d4SSatish Balay       rv[k][j] = 0.;
15979371c9d4SSatish Balay     }
1598c4762a1bSJed Brown   rho = ueq[0];
1599c4762a1bSJed Brown   /* u = ueq[1]; */
1600c4762a1bSJed Brown   p0       = ueq[2];
1601c4762a1bSJed Brown   csnd     = PetscSqrtReal(gamma * p0 / rho);
1602c4762a1bSJed Brown   lv[0][1] = rho * .5;
1603c4762a1bSJed Brown   lv[0][2] = -.5 / csnd;
1604c4762a1bSJed Brown   lv[1][0] = csnd;
1605c4762a1bSJed Brown   lv[1][2] = -1. / csnd;
1606c4762a1bSJed Brown   lv[2][1] = rho * .5;
1607c4762a1bSJed Brown   lv[2][2] = .5 / csnd;
1608c4762a1bSJed Brown   rv[0][0] = -1. / csnd;
1609c4762a1bSJed Brown   rv[1][0] = 1. / rho;
1610c4762a1bSJed Brown   rv[2][0] = -csnd;
1611c4762a1bSJed Brown   rv[0][1] = 1. / csnd;
1612c4762a1bSJed Brown   rv[0][2] = 1. / csnd;
1613c4762a1bSJed Brown   rv[1][2] = 1. / rho;
1614c4762a1bSJed Brown   rv[2][2] = csnd;
1615c4762a1bSJed Brown   return 0;
1616c4762a1bSJed Brown }
1617c4762a1bSJed Brown 
1618d71ae5a4SJacob Faibussowitsch int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1619d71ae5a4SJacob Faibussowitsch {
1620c4762a1bSJed Brown   PetscReal p0, u0, wcp[3], wc[3];
1621c4762a1bSJed Brown   PetscReal lv[3][3];
1622c4762a1bSJed Brown   PetscReal vp[3];
1623c4762a1bSJed Brown   PetscReal rv[3][3];
1624c4762a1bSJed Brown   PetscReal eps, ueq[3], rho0, twopi;
1625c4762a1bSJed Brown 
1626c4762a1bSJed Brown   /* Function Body */
1627c4762a1bSJed Brown   twopi  = 2. * PETSC_PI;
1628c4762a1bSJed Brown   eps    = 1e-4;    /* perturbation */
1629c4762a1bSJed Brown   rho0   = 1e3;     /* density of water */
1630c4762a1bSJed Brown   p0     = 101325.; /* init pressure of 1 atm (?) */
1631c4762a1bSJed Brown   u0     = 0.;
1632c4762a1bSJed Brown   ueq[0] = rho0;
1633c4762a1bSJed Brown   ueq[1] = u0;
1634c4762a1bSJed Brown   ueq[2] = p0;
1635c4762a1bSJed Brown   /* Project initial state to characteristic variables */
1636c4762a1bSJed Brown   eigenvectors(rv, lv, ueq, gamma);
1637c4762a1bSJed Brown   projecteqstate(wc, ueq, lv);
1638c4762a1bSJed Brown   wcp[0] = wc[0];
1639c4762a1bSJed Brown   wcp[1] = wc[1];
1640c4762a1bSJed Brown   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1641c4762a1bSJed Brown   projecttoprim(vp, wcp, rv);
1642c4762a1bSJed Brown   ux->r     = vp[0];         /* density */
1643c4762a1bSJed Brown   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1644c4762a1bSJed Brown   ux->ru[1] = 0.;
1645c4762a1bSJed Brown #if defined DIM > 2
1646c4762a1bSJed Brown   if (dim > 2) ux->ru[2] = 0.;
1647c4762a1bSJed Brown #endif
1648c4762a1bSJed Brown   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1649c4762a1bSJed Brown   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1650c4762a1bSJed Brown   return 0;
1651c4762a1bSJed Brown }
1652c4762a1bSJed Brown 
1653c4762a1bSJed Brown /*TEST
1654c4762a1bSJed Brown 
165530602db0SMatthew G. Knepley   testset:
165630602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0
165730602db0SMatthew G. Knepley 
165830602db0SMatthew G. Knepley     test:
165930602db0SMatthew G. Knepley       suffix: adv_2d_tri_0
166030602db0SMatthew G. Knepley       requires: triangle
166130602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166230602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -simplex -dm_refine 3 -dm_plex_faces 1,1 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
166330602db0SMatthew G. Knepley 
166430602db0SMatthew G. Knepley     test:
166530602db0SMatthew G. Knepley       suffix: adv_2d_tri_1
166630602db0SMatthew G. Knepley       requires: triangle
166730602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166830602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -simplex -dm_refine 5 -dm_plex_faces 1,1 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
166930602db0SMatthew G. Knepley 
167030602db0SMatthew G. Knepley     test:
167130602db0SMatthew G. Knepley       suffix: tut_1
167230602db0SMatthew G. Knepley       requires: exodusii
167330602db0SMatthew G. Knepley       nsize: 1
167446ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -
167530602db0SMatthew G. Knepley 
167630602db0SMatthew G. Knepley     test:
167730602db0SMatthew G. Knepley       suffix: tut_2
167830602db0SMatthew G. Knepley       requires: exodusii
167930602db0SMatthew G. Knepley       nsize: 1
168046ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw
168130602db0SMatthew G. Knepley 
168230602db0SMatthew G. Knepley     test:
168330602db0SMatthew G. Knepley       suffix: tut_3
168430602db0SMatthew G. Knepley       requires: exodusii
168530602db0SMatthew G. Knepley       nsize: 4
1686e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -monitor Error -advect_sol_type bump -petscfv_type leastsquares -petsclimiter_type sin
168730602db0SMatthew G. Knepley 
168830602db0SMatthew G. Knepley     test:
168930602db0SMatthew G. Knepley       suffix: tut_4
1690*910b42cbSStefano Zampini       requires: exodusii !single
169130602db0SMatthew G. Knepley       nsize: 4
1692e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -physics sw -monitor Height,Energy -petscfv_type leastsquares -petsclimiter_type minmod
169330602db0SMatthew G. Knepley 
169430602db0SMatthew G. Knepley   testset:
169530602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
169630602db0SMatthew G. Knepley 
1697c4762a1bSJed Brown     # 2D Advection 0-10
1698c4762a1bSJed Brown     test:
1699c4762a1bSJed Brown       suffix: 0
1700c4762a1bSJed Brown       requires: exodusii
170146ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1702c4762a1bSJed Brown 
1703c4762a1bSJed Brown     test:
1704c4762a1bSJed Brown       suffix: 1
1705c4762a1bSJed Brown       requires: exodusii
170630602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1707c4762a1bSJed Brown 
1708c4762a1bSJed Brown     test:
1709c4762a1bSJed Brown       suffix: 2
1710c4762a1bSJed Brown       requires: exodusii
1711c4762a1bSJed Brown       nsize: 2
171246ac1a18SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1713c4762a1bSJed Brown 
1714c4762a1bSJed Brown     test:
1715c4762a1bSJed Brown       suffix: 3
1716c4762a1bSJed Brown       requires: exodusii
1717c4762a1bSJed Brown       nsize: 2
1718e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1719c4762a1bSJed Brown 
1720c4762a1bSJed Brown     test:
1721c4762a1bSJed Brown       suffix: 4
1722c4762a1bSJed Brown       requires: exodusii
17236c2b77d5SStefano Zampini       nsize: 4
17246c2b77d5SStefano Zampini       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple
1725c4762a1bSJed Brown 
1726c4762a1bSJed Brown     test:
1727c4762a1bSJed Brown       suffix: 5
1728c4762a1bSJed Brown       requires: exodusii
172946ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1
1730c4762a1bSJed Brown 
1731c4762a1bSJed Brown     test:
1732c4762a1bSJed Brown       suffix: 7
1733c4762a1bSJed Brown       requires: exodusii
173430602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1735c4762a1bSJed Brown 
1736c4762a1bSJed Brown     test:
1737c4762a1bSJed Brown       suffix: 8
1738c4762a1bSJed Brown       requires: exodusii
1739c4762a1bSJed Brown       nsize: 2
1740e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1741c4762a1bSJed Brown 
1742c4762a1bSJed Brown     test:
1743c4762a1bSJed Brown       suffix: 9
1744c4762a1bSJed Brown       requires: exodusii
1745c4762a1bSJed Brown       nsize: 8
1746e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1747c4762a1bSJed Brown 
1748c4762a1bSJed Brown     test:
1749c4762a1bSJed Brown       suffix: 10
1750c4762a1bSJed Brown       requires: exodusii
175130602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo
1752c4762a1bSJed Brown 
1753c4762a1bSJed Brown   # 2D Shallow water
1754993a5519SMatthew G. Knepley   testset:
1755993a5519SMatthew G. Knepley     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0
1756993a5519SMatthew G. Knepley 
1757c4762a1bSJed Brown     test:
1758c4762a1bSJed Brown       suffix: sw_0
1759c4762a1bSJed Brown       requires: exodusii
1760993a5519SMatthew G. Knepley       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1761993a5519SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1762993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1763993a5519SMatthew G. Knepley             -monitor height,energy
1764993a5519SMatthew G. Knepley 
1765993a5519SMatthew G. Knepley     test:
1766de8de29fSMatthew G. Knepley       suffix: sw_ceed
1767de8de29fSMatthew G. Knepley       requires: exodusii libceed
1768de8de29fSMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1769de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1770de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1771de8de29fSMatthew G. Knepley             -monitor height,energy
1772de8de29fSMatthew G. Knepley 
1773de8de29fSMatthew G. Knepley     test:
1774a7d931c6SMatthew G. Knepley       suffix: sw_ceed_small
1775a7d931c6SMatthew G. Knepley       requires: exodusii libceed
1776a7d931c6SMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \
1777a7d931c6SMatthew G. Knepley             -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_lower 0,1 -dm_plex_box_upper 6.28,3 -dm_plex_box_faces 8,2 \
1778a7d931c6SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1779a7d931c6SMatthew G. Knepley             -monitor height,energy
1780a7d931c6SMatthew G. Knepley 
1781a7d931c6SMatthew G. Knepley     test:
1782993a5519SMatthew G. Knepley       suffix: sw_1
1783993a5519SMatthew G. Knepley       nsize: 2
1784993a5519SMatthew G. Knepley       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1785993a5519SMatthew G. Knepley             -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_faces 24,12 -dm_plex_box_lower 0,1 -dm_plex_box_upper 1,3 -dm_distribute_overlap 1 \
1786993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1787993a5519SMatthew G. Knepley             -monitor height,energy
1788c4762a1bSJed Brown 
17898d2c18caSMukkund Sunjii     test:
17908d2c18caSMukkund Sunjii       suffix: sw_hll
1791993a5519SMatthew G. Knepley       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1792993a5519SMatthew G. Knepley             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1793993a5519SMatthew G. Knepley             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1794993a5519SMatthew G. Knepley             -monitor height,energy
1795993a5519SMatthew G. Knepley 
1796de8de29fSMatthew G. Knepley   # 2D Euler
1797de8de29fSMatthew G. Knepley   testset:
1798de8de29fSMatthew G. Knepley     args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1799de8de29fSMatthew G. Knepley           -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy
1800de8de29fSMatthew G. Knepley 
1801de8de29fSMatthew G. Knepley     test:
1802277e32beSMatthew G. Knepley       suffix: euler_0
1803*910b42cbSStefano Zampini       requires: exodusii !complex !single
1804277e32beSMatthew G. Knepley       args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1805277e32beSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1806277e32beSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1807277e32beSMatthew G. Knepley 
1808277e32beSMatthew G. Knepley     test:
1809de8de29fSMatthew G. Knepley       suffix: euler_ceed
1810de8de29fSMatthew G. Knepley       requires: exodusii libceed
1811de8de29fSMatthew G. Knepley       args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1812de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1813de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1814de8de29fSMatthew G. Knepley 
1815993a5519SMatthew G. Knepley   testset:
1816993a5519SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
18178d2c18caSMukkund Sunjii 
1818c4762a1bSJed Brown     # 2D Advection: p4est
1819c4762a1bSJed Brown     test:
1820c4762a1bSJed Brown       suffix: p4est_advec_2d
1821c4762a1bSJed Brown       requires: p4est
182230602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_type p4est -dm_forest_minimum_refinement 1 -dm_forest_initial_refinement 2 -dm_p4est_refine_pattern hash -dm_forest_maximum_refinement 5
1823c4762a1bSJed Brown 
1824c4762a1bSJed Brown     # Advection in a box
1825c4762a1bSJed Brown     test:
1826c4762a1bSJed Brown       suffix: adv_2d_quad_0
1827c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1828c4762a1bSJed Brown 
1829c4762a1bSJed Brown     test:
1830c4762a1bSJed Brown       suffix: adv_2d_quad_1
1831c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
1832c4762a1bSJed Brown       timeoutfactor: 3
1833c4762a1bSJed Brown 
1834c4762a1bSJed Brown     test:
1835c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_0
1836c4762a1bSJed Brown       requires: p4est
1837c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1838c4762a1bSJed Brown 
1839c4762a1bSJed Brown     test:
1840c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_1
1841c4762a1bSJed Brown       requires: p4est
1842c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1
1843c4762a1bSJed Brown       timeoutfactor: 3
1844c4762a1bSJed Brown 
184553df731dSPierre Jolivet     test: # broken for quad precision
1846c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_adapt_0
184753df731dSPierre Jolivet       requires: p4est !__float128
1848c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_type p4est -dm_plex_separate_marker -grid_bounds -0.5,0.5,-0.5,0.5 -bc_inflow 1,2,4 -bc_outflow 3 -advect_sol_type bump -advect_bump_center 0.25,0 -advect_bump_radius 0.1 -ufv_use_amr -refine_vec_tagger_box 0.005,inf -coarsen_vec_tagger_box 0,1.e-5 -petscfv_type leastsquares -ts_max_time 0.01
1849c4762a1bSJed Brown       timeoutfactor: 3
1850c4762a1bSJed Brown 
1851c4762a1bSJed Brown     test:
1852c4762a1bSJed Brown       suffix: adv_0
1853c4762a1bSJed Brown       requires: exodusii
185446ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/blockcylinder-50.exo -bc_inflow 100,101,200 -bc_outflow 201
1855c4762a1bSJed Brown 
1856c4762a1bSJed Brown     test:
1857c4762a1bSJed Brown       suffix: shock_0
1858c4762a1bSJed Brown       requires: p4est !single !complex
185930602db0SMatthew G. Knepley       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
186030602db0SMatthew G. Knepley       -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 6 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
186130602db0SMatthew G. Knepley       -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
186230602db0SMatthew G. Knepley       -bc_wall 1,2,3,4 -physics euler -eu_type iv_shock -ufv_cfl 10 -eu_alpha 60. -eu_gamma 1.4 -eu_amach 2.02 -eu_rho2 3. \
186330602db0SMatthew G. Knepley       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
186430602db0SMatthew G. Knepley       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
186530602db0SMatthew G. Knepley       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
1866c4762a1bSJed Brown       timeoutfactor: 3
1867c4762a1bSJed Brown 
1868c4762a1bSJed Brown     # Test GLVis visualization of PetscFV fields
1869c4762a1bSJed Brown     test:
1870c4762a1bSJed Brown       suffix: glvis_adv_2d_tet
187130602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
187230602db0SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
187330602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1874c4762a1bSJed Brown 
1875c4762a1bSJed Brown     test:
1876c4762a1bSJed Brown       suffix: glvis_adv_2d_quad
187730602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
187830602db0SMatthew G. Knepley             -dm_refine 5 -dm_plex_separate_marker \
187930602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1880c4762a1bSJed Brown 
1881c4762a1bSJed Brown TEST*/
1882