xref: /petsc/src/ts/tutorials/ex11.c (revision 53df731d4a58a13d32dc697c3d2b74491914e802) 
1c4762a1bSJed Brown static char help[] = "Second Order TVD Finite Volume Example.\n";
2c4762a1bSJed Brown /*F
3c4762a1bSJed Brown 
4c4762a1bSJed Brown We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops
5c4762a1bSJed Brown over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values,
6c4762a1bSJed Brown \begin{equation}
7c4762a1bSJed Brown   f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
8c4762a1bSJed Brown \end{equation}
9c4762a1bSJed Brown and then update the cell values given the cell volume.
10c4762a1bSJed Brown \begin{eqnarray}
11c4762a1bSJed Brown     f^L_i &-=& \frac{f_i}{vol^L} \\
12c4762a1bSJed Brown     f^R_i &+=& \frac{f_i}{vol^R}
13c4762a1bSJed Brown \end{eqnarray}
14c4762a1bSJed Brown 
15c4762a1bSJed Brown As an example, we can consider the shallow water wave equation,
16c4762a1bSJed Brown \begin{eqnarray}
17c4762a1bSJed Brown      h_t + \nabla\cdot \left( uh                              \right) &=& 0 \\
18c4762a1bSJed Brown   (uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0
19c4762a1bSJed Brown \end{eqnarray}
20c4762a1bSJed Brown where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity.
21c4762a1bSJed Brown 
22c4762a1bSJed Brown A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
23c4762a1bSJed Brown \begin{eqnarray}
24c4762a1bSJed Brown   f^{L,R}_h    &=& uh^{L,R} \cdot \hat n \\
25c4762a1bSJed Brown   f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
26c4762a1bSJed Brown   c^{L,R}      &=& \sqrt{g h^{L,R}} \\
27c4762a1bSJed Brown   s            &=& \max\left( \left|\frac{uh^L \cdot \hat n}{h^L}\right| + c^L, \left|\frac{uh^R \cdot \hat n}{h^R}\right| + c^R \right) \\
28c4762a1bSJed Brown   f_i          &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
29c4762a1bSJed Brown \end{eqnarray}
30c4762a1bSJed Brown where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.
31c4762a1bSJed Brown 
32c4762a1bSJed Brown The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial
33c4762a1bSJed Brown over a neighborhood of the given element.
34c4762a1bSJed Brown 
35c4762a1bSJed Brown The mesh is read in from an ExodusII file, usually generated by Cubit.
3652115386SStefano Zampini 
3752115386SStefano Zampini The example also shows how to handle AMR in a time-dependent TS solver.
38c4762a1bSJed Brown F*/
39c4762a1bSJed Brown #include <petscdmplex.h>
40c4762a1bSJed Brown #include <petscdmforest.h>
41c4762a1bSJed Brown #include <petscds.h>
42c4762a1bSJed Brown #include <petscts.h>
43c4762a1bSJed Brown 
44de8de29fSMatthew G. Knepley #include "ex11.h"
45c4762a1bSJed Brown 
46de8de29fSMatthew G. Knepley static PetscFunctionList PhysicsList, PhysicsRiemannList_SW, PhysicsRiemannList_Euler;
47c4762a1bSJed Brown 
48c4762a1bSJed Brown /* 'User' implements a discretization of a continuous model. */
49c4762a1bSJed Brown typedef struct _n_User *User;
50c4762a1bSJed Brown typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
5145480ffeSMatthew G. Knepley typedef PetscErrorCode (*SetUpBCFunction)(DM, PetscDS, Physics);
52c4762a1bSJed Brown typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
53c4762a1bSJed Brown typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
54c4762a1bSJed Brown static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
55c4762a1bSJed Brown static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
56c4762a1bSJed Brown static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);
57c4762a1bSJed Brown 
58c4762a1bSJed Brown typedef struct _n_FunctionalLink *FunctionalLink;
59c4762a1bSJed Brown struct _n_FunctionalLink {
60c4762a1bSJed Brown   char              *name;
61c4762a1bSJed Brown   FunctionalFunction func;
62c4762a1bSJed Brown   void              *ctx;
63c4762a1bSJed Brown   PetscInt           offset;
64c4762a1bSJed Brown   FunctionalLink     next;
65c4762a1bSJed Brown };
66c4762a1bSJed Brown 
67c4762a1bSJed Brown struct _n_Model {
68da81f932SPierre Jolivet   MPI_Comm         comm; /* Does not do collective communication, but some error conditions can be collective */
69c4762a1bSJed Brown   Physics          physics;
70c4762a1bSJed Brown   FunctionalLink   functionalRegistry;
71c4762a1bSJed Brown   PetscInt         maxComputed;
72c4762a1bSJed Brown   PetscInt         numMonitored;
73c4762a1bSJed Brown   FunctionalLink  *functionalMonitored;
74c4762a1bSJed Brown   PetscInt         numCall;
75c4762a1bSJed Brown   FunctionalLink  *functionalCall;
76c4762a1bSJed Brown   SolutionFunction solution;
77c4762a1bSJed Brown   SetUpBCFunction  setupbc;
78c4762a1bSJed Brown   void            *solutionctx;
79c4762a1bSJed Brown   PetscReal        maxspeed; /* estimate of global maximum speed (for CFL calculation) */
80c4762a1bSJed Brown   PetscReal        bounds[2 * DIM];
81c4762a1bSJed Brown   PetscErrorCode (*errorIndicator)(PetscInt, PetscReal, PetscInt, const PetscScalar[], const PetscScalar[], PetscReal *, void *);
82c4762a1bSJed Brown   void *errorCtx;
83de8de29fSMatthew G. Knepley   PetscErrorCode (*setupCEED)(DM, Physics);
84c4762a1bSJed Brown };
85c4762a1bSJed Brown 
86c4762a1bSJed Brown struct _n_User {
87c4762a1bSJed Brown   PetscInt  vtkInterval;                        /* For monitor */
88c4762a1bSJed Brown   char      outputBasename[PETSC_MAX_PATH_LEN]; /* Basename for output files */
89c4762a1bSJed Brown   PetscInt  monitorStepOffset;
90c4762a1bSJed Brown   Model     model;
91c4762a1bSJed Brown   PetscBool vtkmon;
92c4762a1bSJed Brown };
93c4762a1bSJed Brown 
94de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
95de8de29fSMatthew G. Knepley // Free a plain data context that was allocated using PETSc; returning libCEED error codes
96de8de29fSMatthew G. Knepley static int FreeContextPetsc(void *data)
97d71ae5a4SJacob Faibussowitsch {
98de8de29fSMatthew G. Knepley   if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed");
99de8de29fSMatthew G. Knepley   return CEED_ERROR_SUCCESS;
100c4762a1bSJed Brown }
101de8de29fSMatthew G. Knepley #endif
102c4762a1bSJed Brown 
103c4762a1bSJed Brown /******************* Advect ********************/
1049371c9d4SSatish Balay typedef enum {
1059371c9d4SSatish Balay   ADVECT_SOL_TILTED,
1069371c9d4SSatish Balay   ADVECT_SOL_BUMP,
1079371c9d4SSatish Balay   ADVECT_SOL_BUMP_CAVITY
1089371c9d4SSatish Balay } AdvectSolType;
109c4762a1bSJed Brown static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "BUMP_CAVITY", "AdvectSolType", "ADVECT_SOL_", 0};
1109371c9d4SSatish Balay typedef enum {
1119371c9d4SSatish Balay   ADVECT_SOL_BUMP_CONE,
1129371c9d4SSatish Balay   ADVECT_SOL_BUMP_COS
1139371c9d4SSatish Balay } AdvectSolBumpType;
114c4762a1bSJed Brown static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};
115c4762a1bSJed Brown 
116c4762a1bSJed Brown typedef struct {
117c4762a1bSJed Brown   PetscReal wind[DIM];
118c4762a1bSJed Brown } Physics_Advect_Tilted;
119c4762a1bSJed Brown typedef struct {
120c4762a1bSJed Brown   PetscReal         center[DIM];
121c4762a1bSJed Brown   PetscReal         radius;
122c4762a1bSJed Brown   AdvectSolBumpType type;
123c4762a1bSJed Brown } Physics_Advect_Bump;
124c4762a1bSJed Brown 
125c4762a1bSJed Brown typedef struct {
126c4762a1bSJed Brown   PetscReal     inflowState;
127c4762a1bSJed Brown   AdvectSolType soltype;
1289371c9d4SSatish Balay   union
1299371c9d4SSatish Balay   {
130c4762a1bSJed Brown     Physics_Advect_Tilted tilted;
131c4762a1bSJed Brown     Physics_Advect_Bump   bump;
132c4762a1bSJed Brown   } sol;
133c4762a1bSJed Brown   struct {
134c4762a1bSJed Brown     PetscInt Solution;
135c4762a1bSJed Brown     PetscInt Error;
136c4762a1bSJed Brown   } functional;
137c4762a1bSJed Brown } Physics_Advect;
138c4762a1bSJed Brown 
1399371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_Advect[] = {
1409371c9d4SSatish Balay   {"U",  1},
1419371c9d4SSatish Balay   {NULL, 0}
1429371c9d4SSatish Balay };
143c4762a1bSJed Brown 
144d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
145d71ae5a4SJacob Faibussowitsch {
146c4762a1bSJed Brown   Physics         phys   = (Physics)ctx;
147c4762a1bSJed Brown   Physics_Advect *advect = (Physics_Advect *)phys->data;
148c4762a1bSJed Brown 
149c4762a1bSJed Brown   PetscFunctionBeginUser;
150c4762a1bSJed Brown   xG[0] = advect->inflowState;
1513ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
152c4762a1bSJed Brown }
153c4762a1bSJed Brown 
154d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
155d71ae5a4SJacob Faibussowitsch {
156c4762a1bSJed Brown   PetscFunctionBeginUser;
157c4762a1bSJed Brown   xG[0] = xI[0];
1583ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
159c4762a1bSJed Brown }
160c4762a1bSJed Brown 
161d71ae5a4SJacob Faibussowitsch static void PhysicsRiemann_Advect(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, Physics phys)
162d71ae5a4SJacob Faibussowitsch {
163c4762a1bSJed Brown   Physics_Advect *advect = (Physics_Advect *)phys->data;
164c4762a1bSJed Brown   PetscReal       wind[DIM], wn;
165c4762a1bSJed Brown 
166c4762a1bSJed Brown   switch (advect->soltype) {
167c4762a1bSJed Brown   case ADVECT_SOL_TILTED: {
168c4762a1bSJed Brown     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
169c4762a1bSJed Brown     wind[0]                       = tilted->wind[0];
170c4762a1bSJed Brown     wind[1]                       = tilted->wind[1];
171c4762a1bSJed Brown   } break;
172c4762a1bSJed Brown   case ADVECT_SOL_BUMP:
173c4762a1bSJed Brown     wind[0] = -qp[1];
174c4762a1bSJed Brown     wind[1] = qp[0];
175c4762a1bSJed Brown     break;
1769371c9d4SSatish Balay   case ADVECT_SOL_BUMP_CAVITY: {
177c4762a1bSJed Brown     PetscInt  i;
178c4762a1bSJed Brown     PetscReal comp2[3] = {0., 0., 0.}, rad2;
179c4762a1bSJed Brown 
180c4762a1bSJed Brown     rad2 = 0.;
181c4762a1bSJed Brown     for (i = 0; i < dim; i++) {
182c4762a1bSJed Brown       comp2[i] = qp[i] * qp[i];
183c4762a1bSJed Brown       rad2 += comp2[i];
184c4762a1bSJed Brown     }
185c4762a1bSJed Brown 
186c4762a1bSJed Brown     wind[0] = -qp[1];
187c4762a1bSJed Brown     wind[1] = qp[0];
188c4762a1bSJed Brown     if (rad2 > 1.) {
189c4762a1bSJed Brown       PetscInt  maxI     = 0;
190c4762a1bSJed Brown       PetscReal maxComp2 = comp2[0];
191c4762a1bSJed Brown 
192c4762a1bSJed Brown       for (i = 1; i < dim; i++) {
193c4762a1bSJed Brown         if (comp2[i] > maxComp2) {
194c4762a1bSJed Brown           maxI     = i;
195c4762a1bSJed Brown           maxComp2 = comp2[i];
196c4762a1bSJed Brown         }
197c4762a1bSJed Brown       }
198c4762a1bSJed Brown       wind[maxI] = 0.;
199c4762a1bSJed Brown     }
2009371c9d4SSatish Balay   } break;
2019371c9d4SSatish Balay   default: {
202c4762a1bSJed Brown     PetscInt i;
203c4762a1bSJed Brown     for (i = 0; i < DIM; ++i) wind[i] = 0.0;
204c4762a1bSJed Brown   }
20598921bdaSJacob Faibussowitsch     /* default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
206c4762a1bSJed Brown   }
207c4762a1bSJed Brown   wn      = Dot2Real(wind, n);
208c4762a1bSJed Brown   flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
209c4762a1bSJed Brown }
210c4762a1bSJed Brown 
211d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
212d71ae5a4SJacob Faibussowitsch {
213c4762a1bSJed Brown   Physics         phys   = (Physics)ctx;
214c4762a1bSJed Brown   Physics_Advect *advect = (Physics_Advect *)phys->data;
215c4762a1bSJed Brown 
216c4762a1bSJed Brown   PetscFunctionBeginUser;
217c4762a1bSJed Brown   switch (advect->soltype) {
218c4762a1bSJed Brown   case ADVECT_SOL_TILTED: {
219c4762a1bSJed Brown     PetscReal              x0[DIM];
220c4762a1bSJed Brown     Physics_Advect_Tilted *tilted = &advect->sol.tilted;
221c4762a1bSJed Brown     Waxpy2Real(-time, tilted->wind, x, x0);
222c4762a1bSJed Brown     if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
223c4762a1bSJed Brown     else u[0] = advect->inflowState;
224c4762a1bSJed Brown   } break;
225c4762a1bSJed Brown   case ADVECT_SOL_BUMP_CAVITY:
226c4762a1bSJed Brown   case ADVECT_SOL_BUMP: {
227c4762a1bSJed Brown     Physics_Advect_Bump *bump = &advect->sol.bump;
228c4762a1bSJed Brown     PetscReal            x0[DIM], v[DIM], r, cost, sint;
229c4762a1bSJed Brown     cost  = PetscCosReal(time);
230c4762a1bSJed Brown     sint  = PetscSinReal(time);
231c4762a1bSJed Brown     x0[0] = cost * x[0] + sint * x[1];
232c4762a1bSJed Brown     x0[1] = -sint * x[0] + cost * x[1];
233c4762a1bSJed Brown     Waxpy2Real(-1, bump->center, x0, v);
234c4762a1bSJed Brown     r = Norm2Real(v);
235c4762a1bSJed Brown     switch (bump->type) {
236d71ae5a4SJacob Faibussowitsch     case ADVECT_SOL_BUMP_CONE:
237d71ae5a4SJacob Faibussowitsch       u[0] = PetscMax(1 - r / bump->radius, 0);
238d71ae5a4SJacob Faibussowitsch       break;
239d71ae5a4SJacob Faibussowitsch     case ADVECT_SOL_BUMP_COS:
240d71ae5a4SJacob Faibussowitsch       u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
241d71ae5a4SJacob Faibussowitsch       break;
242c4762a1bSJed Brown     }
243c4762a1bSJed Brown   } break;
244d71ae5a4SJacob Faibussowitsch   default:
245d71ae5a4SJacob Faibussowitsch     SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
246c4762a1bSJed Brown   }
2473ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
248c4762a1bSJed Brown }
249c4762a1bSJed Brown 
250d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, void *ctx)
251d71ae5a4SJacob Faibussowitsch {
252c4762a1bSJed Brown   Physics         phys      = (Physics)ctx;
253c4762a1bSJed Brown   Physics_Advect *advect    = (Physics_Advect *)phys->data;
254391faea4SSatish Balay   PetscScalar     yexact[1] = {0.0};
255c4762a1bSJed Brown 
256c4762a1bSJed Brown   PetscFunctionBeginUser;
2579566063dSJacob Faibussowitsch   PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
258c4762a1bSJed Brown   f[advect->functional.Solution] = PetscRealPart(y[0]);
259c4762a1bSJed Brown   f[advect->functional.Error]    = PetscAbsScalar(y[0] - yexact[0]);
2603ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
261c4762a1bSJed Brown }
262c4762a1bSJed Brown 
263d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_Advect(DM dm, PetscDS prob, Physics phys)
264d71ae5a4SJacob Faibussowitsch {
265c4762a1bSJed Brown   const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
26645480ffeSMatthew G. Knepley   DMLabel        label;
267c4762a1bSJed Brown 
268c4762a1bSJed Brown   PetscFunctionBeginUser;
269c4762a1bSJed Brown   /* Register "canned" boundary conditions and defaults for where to apply. */
2709566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
271dd39110bSPierre Jolivet   PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Inflow, NULL, phys, NULL));
272dd39110bSPierre Jolivet   PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Advect_Outflow, NULL, phys, NULL));
2733ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
274c4762a1bSJed Brown }
275c4762a1bSJed Brown 
276d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_Advect(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
277d71ae5a4SJacob Faibussowitsch {
278c4762a1bSJed Brown   Physics_Advect *advect;
279c4762a1bSJed Brown 
280c4762a1bSJed Brown   PetscFunctionBeginUser;
281c4762a1bSJed Brown   phys->field_desc = PhysicsFields_Advect;
282c4762a1bSJed Brown   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Advect;
2839566063dSJacob Faibussowitsch   PetscCall(PetscNew(&advect));
284c4762a1bSJed Brown   phys->data   = advect;
285c4762a1bSJed Brown   mod->setupbc = SetUpBC_Advect;
286c4762a1bSJed Brown 
287d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
288c4762a1bSJed Brown   {
289c4762a1bSJed Brown     PetscInt two = 2, dof = 1;
290c4762a1bSJed Brown     advect->soltype = ADVECT_SOL_TILTED;
2919566063dSJacob Faibussowitsch     PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL));
292c4762a1bSJed Brown     switch (advect->soltype) {
293c4762a1bSJed Brown     case ADVECT_SOL_TILTED: {
294c4762a1bSJed Brown       Physics_Advect_Tilted *tilted = &advect->sol.tilted;
295c4762a1bSJed Brown       two                           = 2;
296c4762a1bSJed Brown       tilted->wind[0]               = 0.0;
297c4762a1bSJed Brown       tilted->wind[1]               = 1.0;
2989566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
299c4762a1bSJed Brown       advect->inflowState = -2.0;
3009566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL));
301c4762a1bSJed Brown       phys->maxspeed = Norm2Real(tilted->wind);
302c4762a1bSJed Brown     } break;
303c4762a1bSJed Brown     case ADVECT_SOL_BUMP_CAVITY:
304c4762a1bSJed Brown     case ADVECT_SOL_BUMP: {
305c4762a1bSJed Brown       Physics_Advect_Bump *bump = &advect->sol.bump;
306c4762a1bSJed Brown       two                       = 2;
307c4762a1bSJed Brown       bump->center[0]           = 2.;
308c4762a1bSJed Brown       bump->center[1]           = 0.;
3099566063dSJacob Faibussowitsch       PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
310c4762a1bSJed Brown       bump->radius = 0.9;
3119566063dSJacob Faibussowitsch       PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL));
312c4762a1bSJed Brown       bump->type = ADVECT_SOL_BUMP_CONE;
3139566063dSJacob Faibussowitsch       PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL));
314c4762a1bSJed Brown       phys->maxspeed = 3.; /* radius of mesh, kludge */
315c4762a1bSJed Brown     } break;
316c4762a1bSJed Brown     }
317c4762a1bSJed Brown   }
318d0609cedSBarry Smith   PetscOptionsHeadEnd();
319c4762a1bSJed Brown   /* Initial/transient solution with default boundary conditions */
3209566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys));
321c4762a1bSJed Brown   /* Register "canned" functionals */
3229566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Solution", &advect->functional.Solution, PhysicsFunctional_Advect, phys));
3239566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys));
3243ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
325c4762a1bSJed Brown }
326c4762a1bSJed Brown 
327c4762a1bSJed Brown /******************* Shallow Water ********************/
3289371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_SW[] = {
3299371c9d4SSatish Balay   {"Height",   1  },
3309371c9d4SSatish Balay   {"Momentum", DIM},
3319371c9d4SSatish Balay   {NULL,       0  }
3329371c9d4SSatish Balay };
333c4762a1bSJed Brown 
334d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
335d71ae5a4SJacob Faibussowitsch {
336c4762a1bSJed Brown   PetscFunctionBeginUser;
337c4762a1bSJed Brown   xG[0] = xI[0];
338c4762a1bSJed Brown   xG[1] = -xI[1];
339c4762a1bSJed Brown   xG[2] = -xI[2];
3403ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
341c4762a1bSJed Brown }
342c4762a1bSJed Brown 
343d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_SW(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
344d71ae5a4SJacob Faibussowitsch {
345c4762a1bSJed Brown   PetscReal dx[2], r, sigma;
346c4762a1bSJed Brown 
347c4762a1bSJed Brown   PetscFunctionBeginUser;
34854c59aa7SJacob Faibussowitsch   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
349c4762a1bSJed Brown   dx[0] = x[0] - 1.5;
350c4762a1bSJed Brown   dx[1] = x[1] - 1.0;
351c4762a1bSJed Brown   r     = Norm2Real(dx);
352c4762a1bSJed Brown   sigma = 0.5;
353c4762a1bSJed Brown   u[0]  = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma)));
354c4762a1bSJed Brown   u[1]  = 0.0;
355c4762a1bSJed Brown   u[2]  = 0.0;
3563ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
357c4762a1bSJed Brown }
358c4762a1bSJed Brown 
359d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
360d71ae5a4SJacob Faibussowitsch {
361c4762a1bSJed Brown   Physics       phys = (Physics)ctx;
362c4762a1bSJed Brown   Physics_SW   *sw   = (Physics_SW *)phys->data;
363c4762a1bSJed Brown   const SWNode *x    = (const SWNode *)xx;
364c4762a1bSJed Brown   PetscReal     u[2];
365c4762a1bSJed Brown   PetscReal     h;
366c4762a1bSJed Brown 
367c4762a1bSJed Brown   PetscFunctionBeginUser;
368c4762a1bSJed Brown   h = x->h;
369c4762a1bSJed Brown   Scale2Real(1. / x->h, x->uh, u);
370c4762a1bSJed Brown   f[sw->functional.Height] = h;
371c4762a1bSJed Brown   f[sw->functional.Speed]  = Norm2Real(u) + PetscSqrtReal(sw->gravity * h);
372c4762a1bSJed Brown   f[sw->functional.Energy] = 0.5 * (Dot2Real(x->uh, u) + sw->gravity * PetscSqr(h));
3733ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
374c4762a1bSJed Brown }
375c4762a1bSJed Brown 
376d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_SW(DM dm, PetscDS prob, Physics phys)
377d71ae5a4SJacob Faibussowitsch {
378c4762a1bSJed Brown   const PetscInt wallids[] = {100, 101, 200, 300};
37945480ffeSMatthew G. Knepley   DMLabel        label;
38045480ffeSMatthew G. Knepley 
381c4762a1bSJed Brown   PetscFunctionBeginUser;
3829566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
383dd39110bSPierre Jolivet   PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_SW_Wall, NULL, phys, NULL));
3843ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
385c4762a1bSJed Brown }
386c4762a1bSJed Brown 
387de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
388de8de29fSMatthew G. Knepley static PetscErrorCode CreateQFunctionContext_SW(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
389de8de29fSMatthew G. Knepley {
390de8de29fSMatthew G. Knepley   Physics_SW *in = (Physics_SW *)phys->data;
391de8de29fSMatthew G. Knepley   Physics_SW *sw;
392de8de29fSMatthew G. Knepley 
393de8de29fSMatthew G. Knepley   PetscFunctionBeginUser;
394de8de29fSMatthew G. Knepley   PetscCall(PetscCalloc1(1, &sw));
395de8de29fSMatthew G. Knepley 
396de8de29fSMatthew G. Knepley   sw->gravity = in->gravity;
397de8de29fSMatthew G. Knepley 
398de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
399de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*sw), sw));
400de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
401de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gravity", offsetof(Physics_SW, gravity), 1, "Accelaration due to gravity"));
402de8de29fSMatthew G. Knepley   PetscFunctionReturn(0);
403de8de29fSMatthew G. Knepley }
404de8de29fSMatthew G. Knepley #endif
405de8de29fSMatthew G. Knepley 
406de8de29fSMatthew G. Knepley static PetscErrorCode SetupCEED_SW(DM dm, Physics physics)
407de8de29fSMatthew G. Knepley {
408de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
409de8de29fSMatthew G. Knepley   Ceed                 ceed;
410de8de29fSMatthew G. Knepley   CeedQFunctionContext qfCtx;
411de8de29fSMatthew G. Knepley #endif
412de8de29fSMatthew G. Knepley 
413de8de29fSMatthew G. Knepley   PetscFunctionBegin;
414de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
415de8de29fSMatthew G. Knepley   PetscCall(DMGetCeed(dm, &ceed));
416de8de29fSMatthew G. Knepley   PetscCall(CreateQFunctionContext_SW(physics, ceed, &qfCtx));
417de8de29fSMatthew G. Knepley   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_SW_Rusanov_CEED, PhysicsRiemann_SW_Rusanov_CEED_loc, qfCtx));
418de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
419de8de29fSMatthew G. Knepley #endif
420de8de29fSMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
421de8de29fSMatthew G. Knepley }
422de8de29fSMatthew G. Knepley 
423d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
424d71ae5a4SJacob Faibussowitsch {
425c4762a1bSJed Brown   Physics_SW *sw;
4268d2c18caSMukkund Sunjii   char        sw_riemann[64] = "rusanov";
427c4762a1bSJed Brown 
428c4762a1bSJed Brown   PetscFunctionBeginUser;
429c4762a1bSJed Brown   phys->field_desc = PhysicsFields_SW;
4309566063dSJacob Faibussowitsch   PetscCall(PetscNew(&sw));
431c4762a1bSJed Brown   phys->data     = sw;
432c4762a1bSJed Brown   mod->setupbc   = SetUpBC_SW;
433de8de29fSMatthew G. Knepley   mod->setupCEED = SetupCEED_SW;
434c4762a1bSJed Brown 
4353ba16761SJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov));
4363ba16761SJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL));
437de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
438de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov_ceed", PhysicsRiemann_SW_Rusanov_CEED));
439de8de29fSMatthew G. Knepley #endif
4408d2c18caSMukkund Sunjii 
441d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
442c4762a1bSJed Brown   {
4438d2c18caSMukkund Sunjii     void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
444c4762a1bSJed Brown     sw->gravity = 1.0;
4459566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
4469566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL));
4479566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW));
4488d2c18caSMukkund Sunjii     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_SW;
449c4762a1bSJed Brown   }
450d0609cedSBarry Smith   PetscOptionsHeadEnd();
451c4762a1bSJed Brown   phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */
452c4762a1bSJed Brown 
4539566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
4549566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
4559566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
4569566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
457c4762a1bSJed Brown 
4583ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
459c4762a1bSJed Brown }
460c4762a1bSJed Brown 
461c4762a1bSJed Brown /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
462c4762a1bSJed Brown /* An initial-value and self-similar solutions of the compressible Euler equations */
463c4762a1bSJed Brown /* Ravi Samtaney and D. I. Pullin */
464c4762a1bSJed Brown /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
4659371c9d4SSatish Balay static const struct FieldDescription PhysicsFields_Euler[] = {
4669371c9d4SSatish Balay   {"Density",  1  },
4679371c9d4SSatish Balay   {"Momentum", DIM},
4689371c9d4SSatish Balay   {"Energy",   1  },
4699371c9d4SSatish Balay   {NULL,       0  }
4709371c9d4SSatish Balay };
471c4762a1bSJed Brown 
472c4762a1bSJed Brown /* initial condition */
473c4762a1bSJed Brown int                   initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
474d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
475d71ae5a4SJacob Faibussowitsch {
476c4762a1bSJed Brown   PetscInt       i;
477c4762a1bSJed Brown   Physics        phys = (Physics)ctx;
478c4762a1bSJed Brown   Physics_Euler *eu   = (Physics_Euler *)phys->data;
479c4762a1bSJed Brown   EulerNode     *uu   = (EulerNode *)u;
480de8de29fSMatthew G. Knepley   PetscReal      p0, gamma, c = 0.;
481c4762a1bSJed Brown   PetscFunctionBeginUser;
48254c59aa7SJacob Faibussowitsch   PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
483c4762a1bSJed Brown 
484c4762a1bSJed Brown   for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
485c4762a1bSJed Brown   /* set E and rho */
486de8de29fSMatthew G. Knepley   gamma = eu->gamma;
487c4762a1bSJed Brown 
488c4762a1bSJed Brown   if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) {
489c4762a1bSJed Brown     /******************* Euler Density Shock ********************/
490c4762a1bSJed Brown     /* On initial-value and self-similar solutions of the compressible Euler equations */
491c4762a1bSJed Brown     /* Ravi Samtaney and D. I. Pullin */
492c4762a1bSJed Brown     /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
493c4762a1bSJed Brown     /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity,  */
494c4762a1bSJed Brown     p0 = 1.;
495de8de29fSMatthew G. Knepley     if (x[0] < 0.0 + x[1] * eu->itana) {
496c4762a1bSJed Brown       if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */
497c4762a1bSJed Brown         PetscReal amach, rho, press, gas1, p1;
498de8de29fSMatthew G. Knepley         amach     = eu->amach;
499c4762a1bSJed Brown         rho       = 1.;
500c4762a1bSJed Brown         press     = p0;
501c4762a1bSJed Brown         p1        = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0));
502c4762a1bSJed Brown         gas1      = (gamma - 1.0) / (gamma + 1.0);
503c4762a1bSJed Brown         uu->r     = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0);
504c4762a1bSJed Brown         uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach);
505c4762a1bSJed Brown         uu->E     = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0];
5069371c9d4SSatish Balay       } else {      /* left of discontinuity (0) */
507c4762a1bSJed Brown         uu->r = 1.; /* rho = 1 */
508c4762a1bSJed Brown         uu->E = p0 / (gamma - 1.0);
509c4762a1bSJed Brown       }
5109371c9d4SSatish Balay     } else { /* right of discontinuity (2) */
511de8de29fSMatthew G. Knepley       uu->r = eu->rhoR;
512c4762a1bSJed Brown       uu->E = p0 / (gamma - 1.0);
513c4762a1bSJed Brown     }
5149371c9d4SSatish Balay   } else if (eu->type == EULER_SHOCK_TUBE) {
515c4762a1bSJed Brown     /* For (x<x0) set (rho,u,p)=(8,0,10) and for (x>x0) set (rho,u,p)=(1,0,1). Choose x0 to the midpoint of the domain in the x-direction. */
516c4762a1bSJed Brown     if (x[0] < 0.0) {
517c4762a1bSJed Brown       uu->r = 8.;
518c4762a1bSJed Brown       uu->E = 10. / (gamma - 1.);
5199371c9d4SSatish Balay     } else {
520c4762a1bSJed Brown       uu->r = 1.;
521c4762a1bSJed Brown       uu->E = 1. / (gamma - 1.);
522c4762a1bSJed Brown     }
5239371c9d4SSatish Balay   } else if (eu->type == EULER_LINEAR_WAVE) {
524c4762a1bSJed Brown     initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
5259371c9d4SSatish Balay   } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type);
526c4762a1bSJed Brown 
527c4762a1bSJed Brown   /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */
528de8de29fSMatthew G. Knepley   PetscCall(SpeedOfSound_PG(gamma, uu, &c));
529c4762a1bSJed Brown   c = (uu->ru[0] / uu->r) + c;
530c4762a1bSJed Brown   if (c > phys->maxspeed) phys->maxspeed = c;
531c4762a1bSJed Brown 
5323ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
533c4762a1bSJed Brown }
534c4762a1bSJed Brown 
535c4762a1bSJed Brown /* PetscReal* => EulerNode* conversion */
536d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, void *ctx)
537d71ae5a4SJacob Faibussowitsch {
538c4762a1bSJed Brown   PetscInt         i;
539c4762a1bSJed Brown   const EulerNode *xI   = (const EulerNode *)a_xI;
540c4762a1bSJed Brown   EulerNode       *xG   = (EulerNode *)a_xG;
541c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
542c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
543c4762a1bSJed Brown   PetscFunctionBeginUser;
544c4762a1bSJed Brown   xG->r = xI->r;                                     /* ghost cell density - same */
545c4762a1bSJed Brown   xG->E = xI->E;                                     /* ghost cell energy - same */
546c4762a1bSJed Brown   if (n[1] != 0.) {                                  /* top and bottom */
547c4762a1bSJed Brown     xG->ru[0] = xI->ru[0];                           /* copy tang to wall */
548c4762a1bSJed Brown     xG->ru[1] = -xI->ru[1];                          /* reflect perp to t/b wall */
5499371c9d4SSatish Balay   } else {                                           /* sides */
550c4762a1bSJed Brown     for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
551c4762a1bSJed Brown   }
552c4762a1bSJed Brown   if (eu->type == EULER_LINEAR_WAVE) { /* debug */
553c4762a1bSJed Brown #if 0
55463a3b9bcSJacob Faibussowitsch     PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]);
555c4762a1bSJed Brown #endif
556c4762a1bSJed Brown   }
5573ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
558c4762a1bSJed Brown }
559c4762a1bSJed Brown 
560d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
561d71ae5a4SJacob Faibussowitsch {
562c4762a1bSJed Brown   Physics          phys = (Physics)ctx;
563c4762a1bSJed Brown   Physics_Euler   *eu   = (Physics_Euler *)phys->data;
564c4762a1bSJed Brown   const EulerNode *x    = (const EulerNode *)xx;
565c4762a1bSJed Brown   PetscReal        p;
566c4762a1bSJed Brown 
567c4762a1bSJed Brown   PetscFunctionBeginUser;
568c4762a1bSJed Brown   f[eu->monitor.Density]  = x->r;
569c4762a1bSJed Brown   f[eu->monitor.Momentum] = NormDIM(x->ru);
570c4762a1bSJed Brown   f[eu->monitor.Energy]   = x->E;
571c4762a1bSJed Brown   f[eu->monitor.Speed]    = NormDIM(x->ru) / x->r;
572de8de29fSMatthew G. Knepley   PetscCall(Pressure_PG(eu->gamma, x, &p));
573c4762a1bSJed Brown   f[eu->monitor.Pressure] = p;
5743ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
575c4762a1bSJed Brown }
576c4762a1bSJed Brown 
577d71ae5a4SJacob Faibussowitsch static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys)
578d71ae5a4SJacob Faibussowitsch {
579c4762a1bSJed Brown   Physics_Euler *eu = (Physics_Euler *)phys->data;
58045480ffeSMatthew G. Knepley   DMLabel        label;
58145480ffeSMatthew G. Knepley 
582362febeeSStefano Zampini   PetscFunctionBeginUser;
5839566063dSJacob Faibussowitsch   PetscCall(DMGetLabel(dm, "Face Sets", &label));
584c4762a1bSJed Brown   if (eu->type == EULER_LINEAR_WAVE) {
585c4762a1bSJed Brown     const PetscInt wallids[] = {100, 101};
586dd39110bSPierre Jolivet     PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
5879371c9d4SSatish Balay   } else {
588c4762a1bSJed Brown     const PetscInt wallids[] = {100, 101, 200, 300};
589dd39110bSPierre Jolivet     PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)(void))PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
590c4762a1bSJed Brown   }
5913ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
592c4762a1bSJed Brown }
593c4762a1bSJed Brown 
594de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
595de8de29fSMatthew G. Knepley static PetscErrorCode CreateQFunctionContext_Euler(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
596de8de29fSMatthew G. Knepley {
597de8de29fSMatthew G. Knepley   Physics_Euler *in = (Physics_Euler *)phys->data;
598de8de29fSMatthew G. Knepley   Physics_Euler *eu;
599de8de29fSMatthew G. Knepley 
600de8de29fSMatthew G. Knepley   PetscFunctionBeginUser;
601de8de29fSMatthew G. Knepley   PetscCall(PetscCalloc1(1, &eu));
602de8de29fSMatthew G. Knepley 
603de8de29fSMatthew G. Knepley   eu->gamma = in->gamma;
604de8de29fSMatthew G. Knepley 
605de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
606de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*eu), eu));
607de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
608de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gamma", offsetof(Physics_Euler, gamma), 1, "Heat capacity ratio"));
609de8de29fSMatthew G. Knepley   PetscFunctionReturn(0);
610de8de29fSMatthew G. Knepley }
611de8de29fSMatthew G. Knepley #endif
612de8de29fSMatthew G. Knepley 
613de8de29fSMatthew G. Knepley static PetscErrorCode SetupCEED_Euler(DM dm, Physics physics)
614de8de29fSMatthew G. Knepley {
615de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
616de8de29fSMatthew G. Knepley   Ceed                 ceed;
617de8de29fSMatthew G. Knepley   CeedQFunctionContext qfCtx;
618de8de29fSMatthew G. Knepley #endif
619de8de29fSMatthew G. Knepley 
620de8de29fSMatthew G. Knepley   PetscFunctionBegin;
621de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
622de8de29fSMatthew G. Knepley   PetscCall(DMGetCeed(dm, &ceed));
623de8de29fSMatthew G. Knepley   PetscCall(CreateQFunctionContext_Euler(physics, ceed, &qfCtx));
624de8de29fSMatthew G. Knepley   PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_Euler_Godunov_CEED, PhysicsRiemann_Euler_Godunov_CEED_loc, qfCtx));
625de8de29fSMatthew G. Knepley   PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
626de8de29fSMatthew G. Knepley #endif
627de8de29fSMatthew G. Knepley   PetscFunctionReturn(PETSC_SUCCESS);
628de8de29fSMatthew G. Knepley }
629de8de29fSMatthew G. Knepley 
630d71ae5a4SJacob Faibussowitsch static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems *PetscOptionsObject)
631d71ae5a4SJacob Faibussowitsch {
632c4762a1bSJed Brown   Physics_Euler *eu;
633c4762a1bSJed Brown 
634c4762a1bSJed Brown   PetscFunctionBeginUser;
635c4762a1bSJed Brown   phys->field_desc = PhysicsFields_Euler;
636c4762a1bSJed Brown   phys->riemann    = (PetscRiemannFunc)PhysicsRiemann_Euler_Godunov;
6379566063dSJacob Faibussowitsch   PetscCall(PetscNew(&eu));
638c4762a1bSJed Brown   phys->data     = eu;
639c4762a1bSJed Brown   mod->setupbc   = SetUpBC_Euler;
640de8de29fSMatthew G. Knepley   mod->setupCEED = SetupCEED_Euler;
641de8de29fSMatthew G. Knepley 
642de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov", PhysicsRiemann_Euler_Godunov));
643de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
644de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov_ceed", PhysicsRiemann_Euler_Godunov_CEED));
645de8de29fSMatthew G. Knepley #endif
646de8de29fSMatthew G. Knepley 
647d0609cedSBarry Smith   PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
648c4762a1bSJed Brown   {
649de8de29fSMatthew G. Knepley     void (*PhysicsRiemann_Euler)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
650c4762a1bSJed Brown     PetscReal alpha;
651c4762a1bSJed Brown     char      type[64]       = "linear_wave";
652de8de29fSMatthew G. Knepley     char      eu_riemann[64] = "godunov";
653c4762a1bSJed Brown     PetscBool is;
654de8de29fSMatthew G. Knepley     eu->gamma = 1.4;
655de8de29fSMatthew G. Knepley     eu->amach = 2.02;
656de8de29fSMatthew G. Knepley     eu->rhoR  = 3.0;
657de8de29fSMatthew G. Knepley     eu->itana = 0.57735026918963; /* angle of Euler self similar (SS) shock */
658de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsFList("-eu_riemann", "Riemann solver", "", PhysicsRiemannList_Euler, eu_riemann, eu_riemann, sizeof eu_riemann, NULL));
659de8de29fSMatthew G. Knepley     PetscCall(PetscFunctionListFind(PhysicsRiemannList_Euler, eu_riemann, &PhysicsRiemann_Euler));
660de8de29fSMatthew G. Knepley     phys->riemann = (PetscRiemannFunc)PhysicsRiemann_Euler;
661de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->gamma, &eu->gamma, NULL));
662de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->amach, &eu->amach, NULL));
663de8de29fSMatthew G. Knepley     PetscCall(PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->rhoR, &eu->rhoR, NULL));
664c4762a1bSJed Brown     alpha = 60.;
6659566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL));
66663a3b9bcSJacob Faibussowitsch     PetscCheck(alpha > 0. && alpha <= 90., PETSC_COMM_WORLD, PETSC_ERR_SUP, "Alpha bust be > 0 and <= 90 (%g)", (double)alpha);
667de8de29fSMatthew G. Knepley     eu->itana = 1. / PetscTanReal(alpha * PETSC_PI / 180.0);
6689566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL));
6699566063dSJacob Faibussowitsch     PetscCall(PetscStrcmp(type, "linear_wave", &is));
670c4762a1bSJed Brown     if (is) {
67130602db0SMatthew G. Knepley       /* Remember this should be periodic */
672c4762a1bSJed Brown       eu->type = EULER_LINEAR_WAVE;
6739566063dSJacob Faibussowitsch       PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave"));
6749371c9d4SSatish Balay     } else {
67554c59aa7SJacob Faibussowitsch       PetscCheck(DIM == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "DIM must be 2 unless linear wave test %s", type);
6769566063dSJacob Faibussowitsch       PetscCall(PetscStrcmp(type, "iv_shock", &is));
677c4762a1bSJed Brown       if (is) {
678c4762a1bSJed Brown         eu->type = EULER_IV_SHOCK;
6799566063dSJacob Faibussowitsch         PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock"));
6809371c9d4SSatish Balay       } else {
6819566063dSJacob Faibussowitsch         PetscCall(PetscStrcmp(type, "ss_shock", &is));
682c4762a1bSJed Brown         if (is) {
683c4762a1bSJed Brown           eu->type = EULER_SS_SHOCK;
6849566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock"));
6859371c9d4SSatish Balay         } else {
6869566063dSJacob Faibussowitsch           PetscCall(PetscStrcmp(type, "shock_tube", &is));
687c4762a1bSJed Brown           if (is) eu->type = EULER_SHOCK_TUBE;
68898921bdaSJacob Faibussowitsch           else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type);
6899566063dSJacob Faibussowitsch           PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube"));
690c4762a1bSJed Brown         }
691c4762a1bSJed Brown       }
692c4762a1bSJed Brown     }
693c4762a1bSJed Brown   }
694d0609cedSBarry Smith   PetscOptionsHeadEnd();
695c4762a1bSJed Brown   phys->maxspeed = 0.; /* will get set in solution */
6969566063dSJacob Faibussowitsch   PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
6979566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
6989566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
6999566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
7009566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
7019566063dSJacob Faibussowitsch   PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
702c4762a1bSJed Brown 
7033ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
704c4762a1bSJed Brown }
705c4762a1bSJed Brown 
706d71ae5a4SJacob Faibussowitsch static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, void *ctx)
707d71ae5a4SJacob Faibussowitsch {
708c4762a1bSJed Brown   PetscReal err = 0.;
709c4762a1bSJed Brown   PetscInt  i, j;
710c4762a1bSJed Brown 
711c4762a1bSJed Brown   PetscFunctionBeginUser;
712c4762a1bSJed Brown   for (i = 0; i < numComps; i++) {
713ad540459SPierre Jolivet     for (j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
714c4762a1bSJed Brown   }
715c4762a1bSJed Brown   *error = volume * err;
7163ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
717c4762a1bSJed Brown }
718c4762a1bSJed Brown 
719d71ae5a4SJacob Faibussowitsch PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
720d71ae5a4SJacob Faibussowitsch {
721c4762a1bSJed Brown   PetscSF      sfPoint;
722c4762a1bSJed Brown   PetscSection coordSection;
723c4762a1bSJed Brown   Vec          coordinates;
724c4762a1bSJed Brown   PetscSection sectionCell;
725c4762a1bSJed Brown   PetscScalar *part;
726c4762a1bSJed Brown   PetscInt     cStart, cEnd, c;
727c4762a1bSJed Brown   PetscMPIInt  rank;
728c4762a1bSJed Brown 
729c4762a1bSJed Brown   PetscFunctionBeginUser;
7309566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7319566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7329566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, dmCell));
7339566063dSJacob Faibussowitsch   PetscCall(DMGetPointSF(dm, &sfPoint));
7349566063dSJacob Faibussowitsch   PetscCall(DMSetPointSF(*dmCell, sfPoint));
7359566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
7369566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
7379566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
7389566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionCell));
7399566063dSJacob Faibussowitsch   PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
7409566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
74148a46eb9SPierre Jolivet   for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
7429566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionCell));
7439566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(*dmCell, sectionCell));
7449566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionCell));
7459566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(*dmCell, partition));
7469566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
7479566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*partition, &part));
748c4762a1bSJed Brown   for (c = cStart; c < cEnd; ++c) {
749c4762a1bSJed Brown     PetscScalar *p;
750c4762a1bSJed Brown 
7519566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
752c4762a1bSJed Brown     p[0] = rank;
753c4762a1bSJed Brown   }
7549566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*partition, &part));
7553ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
756c4762a1bSJed Brown }
757c4762a1bSJed Brown 
758d71ae5a4SJacob Faibussowitsch PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
759d71ae5a4SJacob Faibussowitsch {
7603e9753d6SMatthew G. Knepley   DM                 plex, dmMass, dmFace, dmCell, dmCoord;
761c4762a1bSJed Brown   PetscSection       coordSection;
762c4762a1bSJed Brown   Vec                coordinates, facegeom, cellgeom;
763c4762a1bSJed Brown   PetscSection       sectionMass;
764c4762a1bSJed Brown   PetscScalar       *m;
765c4762a1bSJed Brown   const PetscScalar *fgeom, *cgeom, *coords;
766c4762a1bSJed Brown   PetscInt           vStart, vEnd, v;
767c4762a1bSJed Brown 
768c4762a1bSJed Brown   PetscFunctionBeginUser;
7699566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
7709566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateSection(dm, &coordSection));
7719566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
7729566063dSJacob Faibussowitsch   PetscCall(DMClone(dm, &dmMass));
7739566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
7749566063dSJacob Faibussowitsch   PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
7759566063dSJacob Faibussowitsch   PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &sectionMass));
7769566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
7779566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
778c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
779c4762a1bSJed Brown     PetscInt numFaces;
780c4762a1bSJed Brown 
7819566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
7829566063dSJacob Faibussowitsch     PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
783c4762a1bSJed Brown   }
7849566063dSJacob Faibussowitsch   PetscCall(PetscSectionSetUp(sectionMass));
7859566063dSJacob Faibussowitsch   PetscCall(DMSetLocalSection(dmMass, sectionMass));
7869566063dSJacob Faibussowitsch   PetscCall(PetscSectionDestroy(&sectionMass));
7879566063dSJacob Faibussowitsch   PetscCall(DMGetLocalVector(dmMass, massMatrix));
7889566063dSJacob Faibussowitsch   PetscCall(VecGetArray(*massMatrix, &m));
7899566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
7909566063dSJacob Faibussowitsch   PetscCall(VecGetDM(facegeom, &dmFace));
7919566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(facegeom, &fgeom));
7929566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellgeom, &dmCell));
7939566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellgeom, &cgeom));
7949566063dSJacob Faibussowitsch   PetscCall(DMGetCoordinateDM(dm, &dmCoord));
7959566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(coordinates, &coords));
796c4762a1bSJed Brown   for (v = vStart; v < vEnd; ++v) {
797c4762a1bSJed Brown     const PetscInt  *faces;
798c4762a1bSJed Brown     PetscFVFaceGeom *fgA, *fgB, *cg;
799c4762a1bSJed Brown     PetscScalar     *vertex;
800c4762a1bSJed Brown     PetscInt         numFaces, sides[2], f, g;
801c4762a1bSJed Brown 
8029566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
8039566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
8049566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSupport(dmMass, v, &faces));
805c4762a1bSJed Brown     for (f = 0; f < numFaces; ++f) {
806c4762a1bSJed Brown       sides[0] = faces[f];
8079566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
808c4762a1bSJed Brown       for (g = 0; g < numFaces; ++g) {
809c4762a1bSJed Brown         const PetscInt *cells = NULL;
810c4762a1bSJed Brown         PetscReal       area  = 0.0;
811c4762a1bSJed Brown         PetscInt        numCells;
812c4762a1bSJed Brown 
813c4762a1bSJed Brown         sides[1] = faces[g];
8149566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
8159566063dSJacob Faibussowitsch         PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
81654c59aa7SJacob Faibussowitsch         PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
8179566063dSJacob Faibussowitsch         PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
818c4762a1bSJed 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]));
819c4762a1bSJed 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]));
820c4762a1bSJed Brown         m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
8219566063dSJacob Faibussowitsch         PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
822c4762a1bSJed Brown       }
823c4762a1bSJed Brown     }
824c4762a1bSJed Brown   }
8259566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
8269566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
8279566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(coordinates, &coords));
8289566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(*massMatrix, &m));
8299566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dmMass));
8309566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
8313ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
832c4762a1bSJed Brown }
833c4762a1bSJed Brown 
834c4762a1bSJed Brown /* Behavior will be different for multi-physics or when using non-default boundary conditions */
835d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
836d71ae5a4SJacob Faibussowitsch {
837c4762a1bSJed Brown   PetscFunctionBeginUser;
838c4762a1bSJed Brown   mod->solution    = func;
839c4762a1bSJed Brown   mod->solutionctx = ctx;
8403ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
841c4762a1bSJed Brown }
842c4762a1bSJed Brown 
843d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
844d71ae5a4SJacob Faibussowitsch {
845c4762a1bSJed Brown   FunctionalLink link, *ptr;
846c4762a1bSJed Brown   PetscInt       lastoffset = -1;
847c4762a1bSJed Brown 
848c4762a1bSJed Brown   PetscFunctionBeginUser;
849c4762a1bSJed Brown   for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
8509566063dSJacob Faibussowitsch   PetscCall(PetscNew(&link));
8519566063dSJacob Faibussowitsch   PetscCall(PetscStrallocpy(name, &link->name));
852c4762a1bSJed Brown   link->offset = lastoffset + 1;
853c4762a1bSJed Brown   link->func   = func;
854c4762a1bSJed Brown   link->ctx    = ctx;
855c4762a1bSJed Brown   link->next   = NULL;
856c4762a1bSJed Brown   *ptr         = link;
857c4762a1bSJed Brown   *offset      = link->offset;
8583ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
859c4762a1bSJed Brown }
860c4762a1bSJed Brown 
861d71ae5a4SJacob Faibussowitsch static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems *PetscOptionsObject)
862d71ae5a4SJacob Faibussowitsch {
863c4762a1bSJed Brown   PetscInt       i, j;
864c4762a1bSJed Brown   FunctionalLink link;
865c4762a1bSJed Brown   char          *names[256];
866c4762a1bSJed Brown 
867c4762a1bSJed Brown   PetscFunctionBeginUser;
868dd39110bSPierre Jolivet   mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
8699566063dSJacob Faibussowitsch   PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
870c4762a1bSJed Brown   /* Create list of functionals that will be computed somehow */
8719566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
872c4762a1bSJed Brown   /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
8739566063dSJacob Faibussowitsch   PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
874c4762a1bSJed Brown   mod->numCall = 0;
875c4762a1bSJed Brown   for (i = 0; i < mod->numMonitored; i++) {
876c4762a1bSJed Brown     for (link = mod->functionalRegistry; link; link = link->next) {
877c4762a1bSJed Brown       PetscBool match;
8789566063dSJacob Faibussowitsch       PetscCall(PetscStrcasecmp(names[i], link->name, &match));
879c4762a1bSJed Brown       if (match) break;
880c4762a1bSJed Brown     }
88154c59aa7SJacob Faibussowitsch     PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
882c4762a1bSJed Brown     mod->functionalMonitored[i] = link;
883c4762a1bSJed Brown     for (j = 0; j < i; j++) {
884c4762a1bSJed Brown       if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
885c4762a1bSJed Brown     }
886c4762a1bSJed Brown     mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
887c4762a1bSJed Brown   next_name:
8889566063dSJacob Faibussowitsch     PetscCall(PetscFree(names[i]));
889c4762a1bSJed Brown   }
890c4762a1bSJed Brown 
891c4762a1bSJed Brown   /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
892c4762a1bSJed Brown   mod->maxComputed = -1;
893c4762a1bSJed Brown   for (link = mod->functionalRegistry; link; link = link->next) {
894c4762a1bSJed Brown     for (i = 0; i < mod->numCall; i++) {
895c4762a1bSJed Brown       FunctionalLink call = mod->functionalCall[i];
896ad540459SPierre Jolivet       if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
897c4762a1bSJed Brown     }
898c4762a1bSJed Brown   }
8993ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
900c4762a1bSJed Brown }
901c4762a1bSJed Brown 
902d71ae5a4SJacob Faibussowitsch static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
903d71ae5a4SJacob Faibussowitsch {
904c4762a1bSJed Brown   FunctionalLink l, next;
905c4762a1bSJed Brown 
906c4762a1bSJed Brown   PetscFunctionBeginUser;
9073ba16761SJacob Faibussowitsch   if (!link) PetscFunctionReturn(PETSC_SUCCESS);
908c4762a1bSJed Brown   l     = *link;
909c4762a1bSJed Brown   *link = NULL;
910c4762a1bSJed Brown   for (; l; l = next) {
911c4762a1bSJed Brown     next = l->next;
9129566063dSJacob Faibussowitsch     PetscCall(PetscFree(l->name));
9139566063dSJacob Faibussowitsch     PetscCall(PetscFree(l));
914c4762a1bSJed Brown   }
9153ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
916c4762a1bSJed Brown }
917c4762a1bSJed Brown 
918c4762a1bSJed Brown /* put the solution callback into a functional callback */
919d71ae5a4SJacob Faibussowitsch static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
920d71ae5a4SJacob Faibussowitsch {
921c4762a1bSJed Brown   Model mod;
9227510d9b0SBarry Smith   PetscFunctionBeginUser;
923c4762a1bSJed Brown   mod = (Model)modctx;
9249566063dSJacob Faibussowitsch   PetscCall((*mod->solution)(mod, time, x, u, mod->solutionctx));
9253ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
926c4762a1bSJed Brown }
927c4762a1bSJed Brown 
928d71ae5a4SJacob Faibussowitsch PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
929d71ae5a4SJacob Faibussowitsch {
930c4762a1bSJed Brown   PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx);
931c4762a1bSJed Brown   void *ctx[1];
932c4762a1bSJed Brown   Model mod = user->model;
933c4762a1bSJed Brown 
934c4762a1bSJed Brown   PetscFunctionBeginUser;
935c4762a1bSJed Brown   func[0] = SolutionFunctional;
936c4762a1bSJed Brown   ctx[0]  = (void *)mod;
9379566063dSJacob Faibussowitsch   PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X));
9383ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
939c4762a1bSJed Brown }
940c4762a1bSJed Brown 
941d71ae5a4SJacob Faibussowitsch static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
942d71ae5a4SJacob Faibussowitsch {
943c4762a1bSJed Brown   PetscFunctionBeginUser;
9449566063dSJacob Faibussowitsch   PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
9459566063dSJacob Faibussowitsch   PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
9469566063dSJacob Faibussowitsch   PetscCall(PetscViewerFileSetName(*viewer, filename));
9473ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
948c4762a1bSJed Brown }
949c4762a1bSJed Brown 
950d71ae5a4SJacob Faibussowitsch static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
951d71ae5a4SJacob Faibussowitsch {
952c4762a1bSJed Brown   User        user = (User)ctx;
9533e9753d6SMatthew G. Knepley   DM          dm, plex;
954c4762a1bSJed Brown   PetscViewer viewer;
955c4762a1bSJed Brown   char        filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
956c4762a1bSJed Brown   PetscReal   xnorm;
957c4762a1bSJed Brown 
958c4762a1bSJed Brown   PetscFunctionBeginUser;
9599566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "u"));
9609566063dSJacob Faibussowitsch   PetscCall(VecGetDM(X, &dm));
9619566063dSJacob Faibussowitsch   PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));
962c4762a1bSJed Brown 
963ad540459SPierre Jolivet   if (stepnum >= 0) stepnum += user->monitorStepOffset;
964c4762a1bSJed Brown   if (stepnum >= 0) { /* No summary for final time */
965c4762a1bSJed Brown     Model              mod = user->model;
9663e9753d6SMatthew G. Knepley     Vec                cellgeom;
967c4762a1bSJed Brown     PetscInt           c, cStart, cEnd, fcount, i;
968c4762a1bSJed Brown     size_t             ftableused, ftablealloc;
969c4762a1bSJed Brown     const PetscScalar *cgeom, *x;
970c4762a1bSJed Brown     DM                 dmCell;
971c4762a1bSJed Brown     DMLabel            vtkLabel;
972c4762a1bSJed Brown     PetscReal         *fmin, *fmax, *fintegral, *ftmp;
9733e9753d6SMatthew G. Knepley 
9749566063dSJacob Faibussowitsch     PetscCall(DMConvert(dm, DMPLEX, &plex));
9759566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
976c4762a1bSJed Brown     fcount = mod->maxComputed + 1;
9779566063dSJacob Faibussowitsch     PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
978c4762a1bSJed Brown     for (i = 0; i < fcount; i++) {
979c4762a1bSJed Brown       fmin[i]      = PETSC_MAX_REAL;
980c4762a1bSJed Brown       fmax[i]      = PETSC_MIN_REAL;
981c4762a1bSJed Brown       fintegral[i] = 0;
982c4762a1bSJed Brown     }
9839566063dSJacob Faibussowitsch     PetscCall(VecGetDM(cellgeom, &dmCell));
9849566063dSJacob Faibussowitsch     PetscCall(DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd));
9859566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(cellgeom, &cgeom));
9869566063dSJacob Faibussowitsch     PetscCall(VecGetArrayRead(X, &x));
9879566063dSJacob Faibussowitsch     PetscCall(DMGetLabel(dm, "vtk", &vtkLabel));
988c4762a1bSJed Brown     for (c = cStart; c < cEnd; ++c) {
989c4762a1bSJed Brown       PetscFVCellGeom   *cg;
990c4762a1bSJed Brown       const PetscScalar *cx     = NULL;
991c4762a1bSJed Brown       PetscInt           vtkVal = 0;
992c4762a1bSJed Brown 
993c4762a1bSJed Brown       /* not that these two routines as currently implemented work for any dm with a
994c4762a1bSJed Brown        * localSection/globalSection */
9959566063dSJacob Faibussowitsch       PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
9969566063dSJacob Faibussowitsch       PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
9979566063dSJacob Faibussowitsch       if (vtkLabel) PetscCall(DMLabelGetValue(vtkLabel, c, &vtkVal));
998c4762a1bSJed Brown       if (!vtkVal || !cx) continue; /* ghost, or not a global cell */
999c4762a1bSJed Brown       for (i = 0; i < mod->numCall; i++) {
1000c4762a1bSJed Brown         FunctionalLink flink = mod->functionalCall[i];
10019566063dSJacob Faibussowitsch         PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
1002c4762a1bSJed Brown       }
1003c4762a1bSJed Brown       for (i = 0; i < fcount; i++) {
1004c4762a1bSJed Brown         fmin[i] = PetscMin(fmin[i], ftmp[i]);
1005c4762a1bSJed Brown         fmax[i] = PetscMax(fmax[i], ftmp[i]);
1006c4762a1bSJed Brown         fintegral[i] += cg->volume * ftmp[i];
1007c4762a1bSJed Brown       }
1008c4762a1bSJed Brown     }
10099566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
10109566063dSJacob Faibussowitsch     PetscCall(VecRestoreArrayRead(X, &x));
10119566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&plex));
1012712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1013712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1014712fec58SPierre Jolivet     PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
1015c4762a1bSJed Brown 
1016c4762a1bSJed Brown     ftablealloc = fcount * 100;
1017c4762a1bSJed Brown     ftableused  = 0;
10189566063dSJacob Faibussowitsch     PetscCall(PetscMalloc1(ftablealloc, &ftable));
1019c4762a1bSJed Brown     for (i = 0; i < mod->numMonitored; i++) {
1020c4762a1bSJed Brown       size_t         countused;
1021c4762a1bSJed Brown       char           buffer[256], *p;
1022c4762a1bSJed Brown       FunctionalLink flink = mod->functionalMonitored[i];
1023c4762a1bSJed Brown       PetscInt       id    = flink->offset;
1024c4762a1bSJed Brown       if (i % 3) {
10259566063dSJacob Faibussowitsch         PetscCall(PetscArraycpy(buffer, "  ", 2));
1026c4762a1bSJed Brown         p = buffer + 2;
1027c4762a1bSJed Brown       } else if (i) {
1028c4762a1bSJed Brown         char newline[] = "\n";
10299566063dSJacob Faibussowitsch         PetscCall(PetscMemcpy(buffer, newline, sizeof(newline) - 1));
1030c4762a1bSJed Brown         p = buffer + sizeof(newline) - 1;
1031c4762a1bSJed Brown       } else {
1032c4762a1bSJed Brown         p = buffer;
1033c4762a1bSJed Brown       }
10349566063dSJacob 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]));
1035c4762a1bSJed Brown       countused--;
1036c4762a1bSJed Brown       countused += p - buffer;
1037c4762a1bSJed Brown       if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1038c4762a1bSJed Brown         char *ftablenew;
1039c4762a1bSJed Brown         ftablealloc = 2 * ftablealloc + countused;
10409566063dSJacob Faibussowitsch         PetscCall(PetscMalloc(ftablealloc, &ftablenew));
10419566063dSJacob Faibussowitsch         PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
10429566063dSJacob Faibussowitsch         PetscCall(PetscFree(ftable));
1043c4762a1bSJed Brown         ftable = ftablenew;
1044c4762a1bSJed Brown       }
10459566063dSJacob Faibussowitsch       PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
1046c4762a1bSJed Brown       ftableused += countused;
1047c4762a1bSJed Brown       ftable[ftableused] = 0;
1048c4762a1bSJed Brown     }
10499566063dSJacob Faibussowitsch     PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));
1050c4762a1bSJed Brown 
105163a3b9bcSJacob Faibussowitsch     PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT "  time %8.4g  |x| %8.4g  %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
10529566063dSJacob Faibussowitsch     PetscCall(PetscFree(ftable));
1053c4762a1bSJed Brown   }
10543ba16761SJacob Faibussowitsch   if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
1055c4762a1bSJed Brown   if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1056c4762a1bSJed Brown     if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
10579566063dSJacob Faibussowitsch       PetscCall(TSGetStepNumber(ts, &stepnum));
1058c4762a1bSJed Brown     }
105963a3b9bcSJacob Faibussowitsch     PetscCall(PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum));
10609566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, filename, &viewer));
10619566063dSJacob Faibussowitsch     PetscCall(VecView(X, viewer));
10629566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
1063c4762a1bSJed Brown   }
10643ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1065c4762a1bSJed Brown }
1066c4762a1bSJed Brown 
1067d71ae5a4SJacob Faibussowitsch static PetscErrorCode initializeTS(DM dm, User user, TS *ts)
1068d71ae5a4SJacob Faibussowitsch {
1069de8de29fSMatthew G. Knepley   PetscBool useCeed;
1070de8de29fSMatthew G. Knepley 
10717510d9b0SBarry Smith   PetscFunctionBeginUser;
10729566063dSJacob Faibussowitsch   PetscCall(TSCreate(PetscObjectComm((PetscObject)dm), ts));
10739566063dSJacob Faibussowitsch   PetscCall(TSSetType(*ts, TSSSP));
10749566063dSJacob Faibussowitsch   PetscCall(TSSetDM(*ts, dm));
10751baa6e33SBarry Smith   if (user->vtkmon) PetscCall(TSMonitorSet(*ts, MonitorVTK, user, NULL));
1076de8de29fSMatthew G. Knepley   PetscCall(DMPlexGetUseCeed(dm, &useCeed));
10779566063dSJacob Faibussowitsch   PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user));
1078de8de29fSMatthew G. Knepley   if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user));
1079de8de29fSMatthew G. Knepley   else PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
10809566063dSJacob Faibussowitsch   PetscCall(TSSetMaxTime(*ts, 2.0));
10819566063dSJacob Faibussowitsch   PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
10823ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1083c4762a1bSJed Brown }
1084c4762a1bSJed Brown 
108552115386SStefano Zampini typedef struct {
108652115386SStefano Zampini   PetscFV      fvm;
108752115386SStefano Zampini   VecTagger    refineTag;
108852115386SStefano Zampini   VecTagger    coarsenTag;
108952115386SStefano Zampini   DM           adaptedDM;
109052115386SStefano Zampini   User         user;
109152115386SStefano Zampini   PetscReal    cfl;
109252115386SStefano Zampini   PetscLimiter limiter;
109352115386SStefano Zampini   PetscLimiter noneLimiter;
109452115386SStefano Zampini } TransferCtx;
109552115386SStefano Zampini 
109652115386SStefano Zampini static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
1097d71ae5a4SJacob Faibussowitsch {
109852115386SStefano Zampini   TransferCtx       *tctx       = (TransferCtx *)ctx;
109952115386SStefano Zampini   PetscFV            fvm        = tctx->fvm;
110052115386SStefano Zampini   VecTagger          refineTag  = tctx->refineTag;
110152115386SStefano Zampini   VecTagger          coarsenTag = tctx->coarsenTag;
110252115386SStefano Zampini   User               user       = tctx->user;
1103c4762a1bSJed Brown   DM                 dm, gradDM, plex, cellDM, adaptedDM = NULL;
1104c4762a1bSJed Brown   Vec                cellGeom, faceGeom;
1105c4762a1bSJed Brown   PetscBool          isForest, computeGradient;
1106c4762a1bSJed Brown   Vec                grad, locGrad, locX, errVec;
1107c4762a1bSJed Brown   PetscInt           cStart, cEnd, c, dim, nRefine, nCoarsen;
110852115386SStefano Zampini   PetscReal          minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1109c4762a1bSJed Brown   PetscScalar       *errArray;
1110c4762a1bSJed Brown   const PetscScalar *pointVals;
1111c4762a1bSJed Brown   const PetscScalar *pointGrads;
1112c4762a1bSJed Brown   const PetscScalar *pointGeom;
1113c4762a1bSJed Brown   DMLabel            adaptLabel = NULL;
1114c4762a1bSJed Brown   IS                 refineIS, coarsenIS;
1115c4762a1bSJed Brown 
11167510d9b0SBarry Smith   PetscFunctionBeginUser;
111752115386SStefano Zampini   *resize = PETSC_FALSE;
11189566063dSJacob Faibussowitsch   PetscCall(VecGetDM(sol, &dm));
11199566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
112052115386SStefano Zampini   PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
11219566063dSJacob Faibussowitsch   PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
11229566063dSJacob Faibussowitsch   PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
11239566063dSJacob Faibussowitsch   PetscCall(DMIsForest(dm, &isForest));
11249566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
11259566063dSJacob Faibussowitsch   PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
11269566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(plex, &locX));
11279566063dSJacob Faibussowitsch   PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
11289566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
11299566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
11309566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(gradDM, &grad));
11319566063dSJacob Faibussowitsch   PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
11329566063dSJacob Faibussowitsch   PetscCall(DMCreateLocalVector(gradDM, &locGrad));
11339566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
11349566063dSJacob Faibussowitsch   PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
11359566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&grad));
11369566063dSJacob Faibussowitsch   PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
11379566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(locGrad, &pointGrads));
11389566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
11399566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(locX, &pointVals));
11409566063dSJacob Faibussowitsch   PetscCall(VecGetDM(cellGeom, &cellDM));
11419566063dSJacob Faibussowitsch   PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
114277433607SBarry Smith   PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
11439566063dSJacob Faibussowitsch   PetscCall(VecSetUp(errVec));
11449566063dSJacob Faibussowitsch   PetscCall(VecGetArray(errVec, &errArray));
1145c4762a1bSJed Brown   for (c = cStart; c < cEnd; c++) {
1146c4762a1bSJed Brown     PetscReal        errInd = 0.;
1147c4762a1bSJed Brown     PetscScalar     *pointGrad;
1148c4762a1bSJed Brown     PetscScalar     *pointVal;
1149c4762a1bSJed Brown     PetscFVCellGeom *cg;
1150c4762a1bSJed Brown 
11519566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
11529566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
11539566063dSJacob Faibussowitsch     PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));
1154c4762a1bSJed Brown 
11559566063dSJacob Faibussowitsch     PetscCall((user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1156c4762a1bSJed Brown     errArray[c - cStart] = errInd;
1157c4762a1bSJed Brown     minMaxInd[0]         = PetscMin(minMaxInd[0], errInd);
1158c4762a1bSJed Brown     minMaxInd[1]         = PetscMax(minMaxInd[1], errInd);
1159c4762a1bSJed Brown   }
11609566063dSJacob Faibussowitsch   PetscCall(VecRestoreArray(errVec, &errArray));
11619566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(locX, &pointVals));
11629566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
11639566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
11649566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&locGrad));
11659566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&locX));
11669566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1167c4762a1bSJed Brown 
11689566063dSJacob Faibussowitsch   PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
11699566063dSJacob Faibussowitsch   PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
11709566063dSJacob Faibussowitsch   PetscCall(ISGetSize(refineIS, &nRefine));
11719566063dSJacob Faibussowitsch   PetscCall(ISGetSize(coarsenIS, &nCoarsen));
11729566063dSJacob Faibussowitsch   if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
11739566063dSJacob Faibussowitsch   if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
11749566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&coarsenIS));
11759566063dSJacob Faibussowitsch   PetscCall(ISDestroy(&refineIS));
11769566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&errVec));
1177c4762a1bSJed Brown 
11789566063dSJacob Faibussowitsch   PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
117952115386SStefano Zampini   PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1180c4762a1bSJed Brown   minMaxInd[1] = -minMaxInd[1];
1181712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
118252115386SStefano Zampini   PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1183c4762a1bSJed Brown   if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
11849566063dSJacob Faibussowitsch     PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1185c4762a1bSJed Brown   }
11869566063dSJacob Faibussowitsch   PetscCall(DMLabelDestroy(&adaptLabel));
1187c4762a1bSJed Brown   if (adaptedDM) {
118852115386SStefano Zampini     tctx->adaptedDM = adaptedDM;
118952115386SStefano 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));
119052115386SStefano Zampini     *resize = PETSC_TRUE;
1191c4762a1bSJed Brown   } else {
119252115386SStefano Zampini     PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1193c4762a1bSJed Brown   }
11943ba16761SJacob Faibussowitsch   PetscFunctionReturn(PETSC_SUCCESS);
1195c4762a1bSJed Brown }
1196c4762a1bSJed Brown 
119752115386SStefano Zampini static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
119852115386SStefano Zampini {
119952115386SStefano Zampini   TransferCtx *tctx = (TransferCtx *)ctx;
120052115386SStefano Zampini   DM           dm;
120152115386SStefano Zampini   PetscReal    time;
120252115386SStefano Zampini 
120352115386SStefano Zampini   PetscFunctionBeginUser;
120452115386SStefano Zampini   PetscCall(TSGetDM(ts, &dm));
120552115386SStefano Zampini   PetscCall(TSGetTime(ts, &time));
120652115386SStefano Zampini   PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
120752115386SStefano Zampini   for (PetscInt i = 0; i < nv; i++) {
120852115386SStefano Zampini     const char *name;
120952115386SStefano Zampini 
121052115386SStefano Zampini     PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
121152115386SStefano Zampini     PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
121252115386SStefano Zampini     PetscCall(PetscObjectGetName((PetscObject)vecsin[i], &name));
121352115386SStefano Zampini     PetscCall(PetscObjectSetName((PetscObject)vecsout[i], name));
121452115386SStefano Zampini   }
121552115386SStefano Zampini   PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */
121652115386SStefano Zampini 
121752115386SStefano Zampini   Model     mod  = tctx->user->model;
121852115386SStefano Zampini   Physics   phys = mod->physics;
121952115386SStefano Zampini   PetscReal minRadius;
122052115386SStefano Zampini 
122152115386SStefano Zampini   PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius));
122252115386SStefano Zampini   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
122352115386SStefano Zampini   PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed");
122452115386SStefano Zampini 
122552115386SStefano Zampini   PetscReal dt = tctx->cfl * minRadius / mod->maxspeed;
122652115386SStefano Zampini   PetscCall(TSSetTimeStep(ts, dt));
122752115386SStefano Zampini 
122852115386SStefano Zampini   PetscCall(TSSetDM(ts, tctx->adaptedDM));
122952115386SStefano Zampini   PetscCall(DMDestroy(&tctx->adaptedDM));
123052115386SStefano Zampini 
123152115386SStefano Zampini   PetscFunctionReturn(PETSC_SUCCESS);
123252115386SStefano Zampini }
123352115386SStefano Zampini 
1234d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
1235d71ae5a4SJacob Faibussowitsch {
1236c4762a1bSJed Brown   MPI_Comm          comm;
1237c4762a1bSJed Brown   PetscDS           prob;
1238c4762a1bSJed Brown   PetscFV           fvm;
1239c4762a1bSJed Brown   PetscLimiter      limiter = NULL, noneLimiter = NULL;
1240c4762a1bSJed Brown   User              user;
1241c4762a1bSJed Brown   Model             mod;
1242c4762a1bSJed Brown   Physics           phys;
12433e9753d6SMatthew G. Knepley   DM                dm, plex;
1244c4762a1bSJed Brown   PetscReal         ftime, cfl, dt, minRadius;
1245c4762a1bSJed Brown   PetscInt          dim, nsteps;
1246c4762a1bSJed Brown   TS                ts;
1247c4762a1bSJed Brown   TSConvergedReason reason;
1248c4762a1bSJed Brown   Vec               X;
1249c4762a1bSJed Brown   PetscViewer       viewer;
1250b5a892a1SMatthew G. Knepley   PetscBool         vtkCellGeom, useAMR;
125130602db0SMatthew G. Knepley   PetscInt          adaptInterval;
1252c4762a1bSJed Brown   char              physname[256] = "advect";
1253c4762a1bSJed Brown   VecTagger         refineTag = NULL, coarsenTag = NULL;
125452115386SStefano Zampini   TransferCtx       tctx;
1255c4762a1bSJed Brown 
1256327415f7SBarry Smith   PetscFunctionBeginUser;
12579566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1258c4762a1bSJed Brown   comm = PETSC_COMM_WORLD;
1259c4762a1bSJed Brown 
12609566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user));
12619566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model));
12629566063dSJacob Faibussowitsch   PetscCall(PetscNew(&user->model->physics));
1263c4762a1bSJed Brown   mod           = user->model;
1264c4762a1bSJed Brown   phys          = mod->physics;
1265c4762a1bSJed Brown   mod->comm     = comm;
1266c4762a1bSJed Brown   useAMR        = PETSC_FALSE;
1267c4762a1bSJed Brown   adaptInterval = 1;
1268c4762a1bSJed Brown 
1269c4762a1bSJed Brown   /* Register physical models to be available on the command line */
12709566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
12719566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
12729566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1273c4762a1bSJed Brown 
1274d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1275c4762a1bSJed Brown   {
1276c4762a1bSJed Brown     cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
12779566063dSJacob Faibussowitsch     PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1278c4762a1bSJed Brown     user->vtkInterval = 1;
12799566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1280c4762a1bSJed Brown     user->vtkmon = PETSC_TRUE;
12819566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1282c4762a1bSJed Brown     vtkCellGeom = PETSC_FALSE;
1283c6a7a370SJeremy L Thompson     PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
12849566063dSJacob Faibussowitsch     PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
12859566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
12869566063dSJacob Faibussowitsch     PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
12879566063dSJacob Faibussowitsch     PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1288c4762a1bSJed Brown   }
1289d0609cedSBarry Smith   PetscOptionsEnd();
1290c4762a1bSJed Brown 
1291c4762a1bSJed Brown   if (useAMR) {
1292c4762a1bSJed Brown     VecTaggerBox refineBox, coarsenBox;
1293c4762a1bSJed Brown 
1294c4762a1bSJed Brown     refineBox.min = refineBox.max = PETSC_MAX_REAL;
1295c4762a1bSJed Brown     coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;
1296c4762a1bSJed Brown 
12979566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &refineTag));
12989566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
12999566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
13009566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
13019566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(refineTag));
13029566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(refineTag));
13039566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
1304c4762a1bSJed Brown 
13059566063dSJacob Faibussowitsch     PetscCall(VecTaggerCreate(comm, &coarsenTag));
13069566063dSJacob Faibussowitsch     PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
13079566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
13089566063dSJacob Faibussowitsch     PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
13099566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetFromOptions(coarsenTag));
13109566063dSJacob Faibussowitsch     PetscCall(VecTaggerSetUp(coarsenTag));
13119566063dSJacob Faibussowitsch     PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1312c4762a1bSJed Brown   }
1313c4762a1bSJed Brown 
1314d0609cedSBarry Smith   PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1315c4762a1bSJed Brown   {
1316c4762a1bSJed Brown     PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
13179566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
13189566063dSJacob Faibussowitsch     PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
13199566063dSJacob Faibussowitsch     PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
13209566063dSJacob Faibussowitsch     PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1321c4762a1bSJed Brown     /* Count number of fields and dofs */
1322c4762a1bSJed Brown     for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
132354c59aa7SJacob Faibussowitsch     PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
13249566063dSJacob Faibussowitsch     PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1325c4762a1bSJed Brown   }
1326d0609cedSBarry Smith   PetscOptionsEnd();
1327c4762a1bSJed Brown 
1328c4762a1bSJed Brown   /* Create mesh */
1329c4762a1bSJed Brown   {
133030602db0SMatthew G. Knepley     PetscInt i;
133130602db0SMatthew G. Knepley 
13329566063dSJacob Faibussowitsch     PetscCall(DMCreate(comm, &dm));
13339566063dSJacob Faibussowitsch     PetscCall(DMSetType(dm, DMPLEX));
13349566063dSJacob Faibussowitsch     PetscCall(DMSetFromOptions(dm));
13359371c9d4SSatish Balay     for (i = 0; i < DIM; i++) {
13369371c9d4SSatish Balay       mod->bounds[2 * i]     = 0.;
13379371c9d4SSatish Balay       mod->bounds[2 * i + 1] = 1.;
13389371c9d4SSatish Balay     };
1339c4762a1bSJed Brown     dim = DIM;
134030602db0SMatthew G. Knepley     { /* a null name means just do a hex box */
134130602db0SMatthew G. Knepley       PetscInt  cells[3] = {1, 1, 1}, n = 3;
134230602db0SMatthew G. Knepley       PetscBool flg2, skew              = PETSC_FALSE;
1343c4762a1bSJed Brown       PetscInt  nret2 = 2 * DIM;
1344d0609cedSBarry Smith       PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
13459566063dSJacob 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));
13469566063dSJacob Faibussowitsch       PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
13479566063dSJacob Faibussowitsch       PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1348d0609cedSBarry Smith       PetscOptionsEnd();
134930602db0SMatthew G. Knepley       /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1350c4762a1bSJed Brown       if (flg2) {
1351c4762a1bSJed Brown         PetscInt     dimEmbed, i;
1352c4762a1bSJed Brown         PetscInt     nCoords;
1353c4762a1bSJed Brown         PetscScalar *coords;
1354c4762a1bSJed Brown         Vec          coordinates;
1355c4762a1bSJed Brown 
13569566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
13579566063dSJacob Faibussowitsch         PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
13589566063dSJacob Faibussowitsch         PetscCall(VecGetLocalSize(coordinates, &nCoords));
135954c59aa7SJacob Faibussowitsch         PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
13609566063dSJacob Faibussowitsch         PetscCall(VecGetArray(coordinates, &coords));
1361c4762a1bSJed Brown         for (i = 0; i < nCoords; i += dimEmbed) {
1362c4762a1bSJed Brown           PetscInt j;
1363c4762a1bSJed Brown 
1364c4762a1bSJed Brown           PetscScalar *coord = &coords[i];
1365c4762a1bSJed Brown           for (j = 0; j < dimEmbed; j++) {
1366c4762a1bSJed Brown             coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1367c4762a1bSJed Brown             if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1368c4762a1bSJed Brown               if (cells[0] == 2 && i == 8) {
1369c4762a1bSJed Brown                 coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
137030602db0SMatthew G. Knepley               } else if (cells[0] == 3) {
1371c4762a1bSJed Brown                 if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1372c4762a1bSJed Brown                 else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1373c4762a1bSJed Brown                 else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1374c4762a1bSJed Brown               }
1375c4762a1bSJed Brown             }
1376c4762a1bSJed Brown           }
1377c4762a1bSJed Brown         }
13789566063dSJacob Faibussowitsch         PetscCall(VecRestoreArray(coordinates, &coords));
13799566063dSJacob Faibussowitsch         PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1380c4762a1bSJed Brown       }
1381c4762a1bSJed Brown     }
1382c4762a1bSJed Brown   }
13839566063dSJacob Faibussowitsch   PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
13849566063dSJacob Faibussowitsch   PetscCall(DMGetDimension(dm, &dim));
1385c4762a1bSJed Brown 
1386c4762a1bSJed Brown   /* set up BCs, functions, tags */
13879566063dSJacob Faibussowitsch   PetscCall(DMCreateLabel(dm, "Face Sets"));
1388c4762a1bSJed Brown   mod->errorIndicator = ErrorIndicator_Simple;
1389c4762a1bSJed Brown 
1390c4762a1bSJed Brown   {
1391c4762a1bSJed Brown     DM gdm;
1392c4762a1bSJed Brown 
13939566063dSJacob Faibussowitsch     PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
13949566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dm));
1395c4762a1bSJed Brown     dm = gdm;
13969566063dSJacob Faibussowitsch     PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1397c4762a1bSJed Brown   }
1398c4762a1bSJed Brown 
13999566063dSJacob Faibussowitsch   PetscCall(PetscFVCreate(comm, &fvm));
14009566063dSJacob Faibussowitsch   PetscCall(PetscFVSetFromOptions(fvm));
14019566063dSJacob Faibussowitsch   PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
14029566063dSJacob Faibussowitsch   PetscCall(PetscFVSetSpatialDimension(fvm, dim));
14039566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1404c4762a1bSJed Brown   {
1405c4762a1bSJed Brown     PetscInt f, dof;
1406c4762a1bSJed Brown     for (f = 0, dof = 0; f < phys->nfields; f++) {
1407c4762a1bSJed Brown       PetscInt newDof = phys->field_desc[f].dof;
1408c4762a1bSJed Brown 
1409c4762a1bSJed Brown       if (newDof == 1) {
14109566063dSJacob Faibussowitsch         PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
14119371c9d4SSatish Balay       } else {
1412c4762a1bSJed Brown         PetscInt j;
1413c4762a1bSJed Brown 
1414c4762a1bSJed Brown         for (j = 0; j < newDof; j++) {
1415c4762a1bSJed Brown           char compName[256] = "Unknown";
1416c4762a1bSJed Brown 
141763a3b9bcSJacob Faibussowitsch           PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
14189566063dSJacob Faibussowitsch           PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1419c4762a1bSJed Brown         }
1420c4762a1bSJed Brown       }
1421c4762a1bSJed Brown       dof += newDof;
1422c4762a1bSJed Brown     }
1423c4762a1bSJed Brown   }
1424c4762a1bSJed Brown   /* FV is now structured with one field having all physics as components */
14259566063dSJacob Faibussowitsch   PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
14269566063dSJacob Faibussowitsch   PetscCall(DMCreateDS(dm));
14279566063dSJacob Faibussowitsch   PetscCall(DMGetDS(dm, &prob));
14289566063dSJacob Faibussowitsch   PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
14299566063dSJacob Faibussowitsch   PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
14309566063dSJacob Faibussowitsch   PetscCall((*mod->setupbc)(dm, prob, phys));
14319566063dSJacob Faibussowitsch   PetscCall(PetscDSSetFromOptions(prob));
1432c4762a1bSJed Brown   {
1433c4762a1bSJed Brown     char      convType[256];
1434c4762a1bSJed Brown     PetscBool flg;
1435c4762a1bSJed Brown 
1436d0609cedSBarry Smith     PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
14379566063dSJacob Faibussowitsch     PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1438d0609cedSBarry Smith     PetscOptionsEnd();
1439c4762a1bSJed Brown     if (flg) {
1440c4762a1bSJed Brown       DM dmConv;
1441c4762a1bSJed Brown 
14429566063dSJacob Faibussowitsch       PetscCall(DMConvert(dm, convType, &dmConv));
1443c4762a1bSJed Brown       if (dmConv) {
14449566063dSJacob Faibussowitsch         PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
14459566063dSJacob Faibussowitsch         PetscCall(DMDestroy(&dm));
1446c4762a1bSJed Brown         dm = dmConv;
14479566063dSJacob Faibussowitsch         PetscCall(DMSetFromOptions(dm));
1448c4762a1bSJed Brown       }
1449c4762a1bSJed Brown     }
1450c4762a1bSJed Brown   }
1451de8de29fSMatthew G. Knepley #ifdef PETSC_HAVE_LIBCEED
1452de8de29fSMatthew G. Knepley   {
1453de8de29fSMatthew G. Knepley     PetscBool useCeed;
1454de8de29fSMatthew G. Knepley     PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1455de8de29fSMatthew G. Knepley     if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1456de8de29fSMatthew G. Knepley   }
1457de8de29fSMatthew G. Knepley #endif
1458c4762a1bSJed Brown 
14599566063dSJacob Faibussowitsch   PetscCall(initializeTS(dm, user, &ts));
1460c4762a1bSJed Brown 
14619566063dSJacob Faibussowitsch   PetscCall(DMCreateGlobalVector(dm, &X));
14629566063dSJacob Faibussowitsch   PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
14639566063dSJacob Faibussowitsch   PetscCall(SetInitialCondition(dm, X, user));
1464c4762a1bSJed Brown   if (useAMR) {
1465c4762a1bSJed Brown     PetscInt adaptIter;
1466c4762a1bSJed Brown 
1467c4762a1bSJed Brown     /* use no limiting when reconstructing gradients for adaptivity */
14689566063dSJacob Faibussowitsch     PetscCall(PetscFVGetLimiter(fvm, &limiter));
14699566063dSJacob Faibussowitsch     PetscCall(PetscObjectReference((PetscObject)limiter));
14709566063dSJacob Faibussowitsch     PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
14719566063dSJacob Faibussowitsch     PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
1472c4762a1bSJed Brown 
147352115386SStefano Zampini     /* Refinement context */
147452115386SStefano Zampini     tctx.fvm         = fvm;
147552115386SStefano Zampini     tctx.refineTag   = refineTag;
147652115386SStefano Zampini     tctx.coarsenTag  = coarsenTag;
147752115386SStefano Zampini     tctx.adaptedDM   = NULL;
147852115386SStefano Zampini     tctx.user        = user;
147952115386SStefano Zampini     tctx.noneLimiter = noneLimiter;
148052115386SStefano Zampini     tctx.limiter     = limiter;
148152115386SStefano Zampini     tctx.cfl         = cfl;
148252115386SStefano Zampini 
148352115386SStefano Zampini     /* Do some initial refinement steps */
1484c4762a1bSJed Brown     for (adaptIter = 0;; ++adaptIter) {
1485c4762a1bSJed Brown       PetscLogDouble bytes;
148652115386SStefano Zampini       PetscBool      resize;
1487c4762a1bSJed Brown 
14889566063dSJacob Faibussowitsch       PetscCall(PetscMemoryGetCurrentUsage(&bytes));
148963a3b9bcSJacob Faibussowitsch       PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
14909566063dSJacob Faibussowitsch       PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
14919566063dSJacob Faibussowitsch       PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));
1492c4762a1bSJed Brown 
149352115386SStefano Zampini       PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
149452115386SStefano Zampini       if (!resize) break;
14959566063dSJacob Faibussowitsch       PetscCall(DMDestroy(&dm));
14969566063dSJacob Faibussowitsch       PetscCall(VecDestroy(&X));
149752115386SStefano Zampini       dm             = tctx.adaptedDM;
149852115386SStefano Zampini       tctx.adaptedDM = NULL;
149952115386SStefano Zampini       PetscCall(TSSetDM(ts, dm));
15009566063dSJacob Faibussowitsch       PetscCall(DMCreateGlobalVector(dm, &X));
15019566063dSJacob Faibussowitsch       PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
15029566063dSJacob Faibussowitsch       PetscCall(SetInitialCondition(dm, X, user));
1503c4762a1bSJed Brown     }
1504c4762a1bSJed Brown   }
1505c4762a1bSJed Brown 
15069566063dSJacob Faibussowitsch   PetscCall(DMConvert(dm, DMPLEX, &plex));
1507c4762a1bSJed Brown   if (vtkCellGeom) {
1508c4762a1bSJed Brown     DM  dmCell;
1509c4762a1bSJed Brown     Vec cellgeom, partition;
1510c4762a1bSJed Brown 
15119566063dSJacob Faibussowitsch     PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
15129566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
15139566063dSJacob Faibussowitsch     PetscCall(VecView(cellgeom, viewer));
15149566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15159566063dSJacob Faibussowitsch     PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
15169566063dSJacob Faibussowitsch     PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
15179566063dSJacob Faibussowitsch     PetscCall(VecView(partition, viewer));
15189566063dSJacob Faibussowitsch     PetscCall(PetscViewerDestroy(&viewer));
15199566063dSJacob Faibussowitsch     PetscCall(VecDestroy(&partition));
15209566063dSJacob Faibussowitsch     PetscCall(DMDestroy(&dmCell));
1521c4762a1bSJed Brown   }
1522c4762a1bSJed Brown   /* collect max maxspeed from all processes -- todo */
15239566063dSJacob Faibussowitsch   PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
15249566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&plex));
1525712fec58SPierre Jolivet   PetscCall(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
152654c59aa7SJacob Faibussowitsch   PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1527c4762a1bSJed Brown   dt = cfl * minRadius / mod->maxspeed;
15289566063dSJacob Faibussowitsch   PetscCall(TSSetTimeStep(ts, dt));
15299566063dSJacob Faibussowitsch   PetscCall(TSSetFromOptions(ts));
1530c4762a1bSJed Brown 
153152115386SStefano Zampini   /* When using adaptive mesh refinement
153252115386SStefano Zampini      specify callbacks to refine the solution
153352115386SStefano Zampini      and interpolate data from old to new mesh */
153452115386SStefano Zampini   if (useAMR) { PetscCall(TSSetResize(ts, adaptToleranceFVMSetUp, Transfer, &tctx)); }
153552115386SStefano Zampini   PetscCall(TSSetSolution(ts, X));
15369566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&X));
153752115386SStefano Zampini   PetscCall(TSSolve(ts, NULL));
15389566063dSJacob Faibussowitsch   PetscCall(TSGetSolveTime(ts, &ftime));
15399566063dSJacob Faibussowitsch   PetscCall(TSGetStepNumber(ts, &nsteps));
15409566063dSJacob Faibussowitsch   PetscCall(TSGetConvergedReason(ts, &reason));
154163a3b9bcSJacob Faibussowitsch   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
15429566063dSJacob Faibussowitsch   PetscCall(TSDestroy(&ts));
1543c4762a1bSJed Brown 
15449566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&refineTag));
15459566063dSJacob Faibussowitsch   PetscCall(VecTaggerDestroy(&coarsenTag));
15469566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsList));
15479566063dSJacob Faibussowitsch   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1548de8de29fSMatthew G. Knepley   PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
15499566063dSJacob Faibussowitsch   PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
15509566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalMonitored));
15519566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->functionalCall));
15529566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics->data));
15539566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model->physics));
15549566063dSJacob Faibussowitsch   PetscCall(PetscFree(user->model));
15559566063dSJacob Faibussowitsch   PetscCall(PetscFree(user));
15569566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&limiter));
15579566063dSJacob Faibussowitsch   PetscCall(PetscLimiterDestroy(&noneLimiter));
15589566063dSJacob Faibussowitsch   PetscCall(PetscFVDestroy(&fvm));
15599566063dSJacob Faibussowitsch   PetscCall(DMDestroy(&dm));
15609566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
1561b122ec5aSJacob Faibussowitsch   return 0;
1562c4762a1bSJed Brown }
1563c4762a1bSJed Brown 
1564c4762a1bSJed Brown /* Subroutine to set up the initial conditions for the */
1565c4762a1bSJed Brown /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1566c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1567d71ae5a4SJacob Faibussowitsch int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1568d71ae5a4SJacob Faibussowitsch {
1569c4762a1bSJed Brown   int j, k;
1570c4762a1bSJed Brown   /*      Wc=matmul(lv,Ueq) 3 vars */
1571c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1572c4762a1bSJed Brown     wc[k] = 0.;
1573ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1574c4762a1bSJed Brown   }
1575c4762a1bSJed Brown   return 0;
1576c4762a1bSJed Brown }
1577c4762a1bSJed Brown /* ----------------------------------------------------------------------- */
1578d71ae5a4SJacob Faibussowitsch int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1579d71ae5a4SJacob Faibussowitsch {
1580c4762a1bSJed Brown   int k, j;
1581c4762a1bSJed Brown   /*      V=matmul(rv,WC) 3 vars */
1582c4762a1bSJed Brown   for (k = 0; k < 3; ++k) {
1583c4762a1bSJed Brown     v[k] = 0.;
1584ad540459SPierre Jolivet     for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1585c4762a1bSJed Brown   }
1586c4762a1bSJed Brown   return 0;
1587c4762a1bSJed Brown }
1588c4762a1bSJed Brown /* ---------------------------------------------------------------------- */
1589d71ae5a4SJacob Faibussowitsch int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1590d71ae5a4SJacob Faibussowitsch {
1591c4762a1bSJed Brown   int       j, k;
1592c4762a1bSJed Brown   PetscReal rho, csnd, p0;
1593c4762a1bSJed Brown   /* PetscScalar u; */
1594c4762a1bSJed Brown 
15959371c9d4SSatish Balay   for (k = 0; k < 3; ++k)
15969371c9d4SSatish Balay     for (j = 0; j < 3; ++j) {
15979371c9d4SSatish Balay       lv[k][j] = 0.;
15989371c9d4SSatish Balay       rv[k][j] = 0.;
15999371c9d4SSatish Balay     }
1600c4762a1bSJed Brown   rho = ueq[0];
1601c4762a1bSJed Brown   /* u = ueq[1]; */
1602c4762a1bSJed Brown   p0       = ueq[2];
1603c4762a1bSJed Brown   csnd     = PetscSqrtReal(gamma * p0 / rho);
1604c4762a1bSJed Brown   lv[0][1] = rho * .5;
1605c4762a1bSJed Brown   lv[0][2] = -.5 / csnd;
1606c4762a1bSJed Brown   lv[1][0] = csnd;
1607c4762a1bSJed Brown   lv[1][2] = -1. / csnd;
1608c4762a1bSJed Brown   lv[2][1] = rho * .5;
1609c4762a1bSJed Brown   lv[2][2] = .5 / csnd;
1610c4762a1bSJed Brown   rv[0][0] = -1. / csnd;
1611c4762a1bSJed Brown   rv[1][0] = 1. / rho;
1612c4762a1bSJed Brown   rv[2][0] = -csnd;
1613c4762a1bSJed Brown   rv[0][1] = 1. / csnd;
1614c4762a1bSJed Brown   rv[0][2] = 1. / csnd;
1615c4762a1bSJed Brown   rv[1][2] = 1. / rho;
1616c4762a1bSJed Brown   rv[2][2] = csnd;
1617c4762a1bSJed Brown   return 0;
1618c4762a1bSJed Brown }
1619c4762a1bSJed Brown 
1620d71ae5a4SJacob Faibussowitsch int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1621d71ae5a4SJacob Faibussowitsch {
1622c4762a1bSJed Brown   PetscReal p0, u0, wcp[3], wc[3];
1623c4762a1bSJed Brown   PetscReal lv[3][3];
1624c4762a1bSJed Brown   PetscReal vp[3];
1625c4762a1bSJed Brown   PetscReal rv[3][3];
1626c4762a1bSJed Brown   PetscReal eps, ueq[3], rho0, twopi;
1627c4762a1bSJed Brown 
1628c4762a1bSJed Brown   /* Function Body */
1629c4762a1bSJed Brown   twopi  = 2. * PETSC_PI;
1630c4762a1bSJed Brown   eps    = 1e-4;    /* perturbation */
1631c4762a1bSJed Brown   rho0   = 1e3;     /* density of water */
1632c4762a1bSJed Brown   p0     = 101325.; /* init pressure of 1 atm (?) */
1633c4762a1bSJed Brown   u0     = 0.;
1634c4762a1bSJed Brown   ueq[0] = rho0;
1635c4762a1bSJed Brown   ueq[1] = u0;
1636c4762a1bSJed Brown   ueq[2] = p0;
1637c4762a1bSJed Brown   /* Project initial state to characteristic variables */
1638c4762a1bSJed Brown   eigenvectors(rv, lv, ueq, gamma);
1639c4762a1bSJed Brown   projecteqstate(wc, ueq, lv);
1640c4762a1bSJed Brown   wcp[0] = wc[0];
1641c4762a1bSJed Brown   wcp[1] = wc[1];
1642c4762a1bSJed Brown   wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1643c4762a1bSJed Brown   projecttoprim(vp, wcp, rv);
1644c4762a1bSJed Brown   ux->r     = vp[0];         /* density */
1645c4762a1bSJed Brown   ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1646c4762a1bSJed Brown   ux->ru[1] = 0.;
1647c4762a1bSJed Brown #if defined DIM > 2
1648c4762a1bSJed Brown   if (dim > 2) ux->ru[2] = 0.;
1649c4762a1bSJed Brown #endif
1650c4762a1bSJed Brown   /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1651c4762a1bSJed Brown   ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1652c4762a1bSJed Brown   return 0;
1653c4762a1bSJed Brown }
1654c4762a1bSJed Brown 
1655c4762a1bSJed Brown /*TEST
1656c4762a1bSJed Brown 
165730602db0SMatthew G. Knepley   testset:
165830602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0
165930602db0SMatthew G. Knepley 
166030602db0SMatthew G. Knepley     test:
166130602db0SMatthew G. Knepley       suffix: adv_2d_tri_0
166230602db0SMatthew G. Knepley       requires: triangle
166330602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
166430602db0SMatthew 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
166530602db0SMatthew G. Knepley 
166630602db0SMatthew G. Knepley     test:
166730602db0SMatthew G. Knepley       suffix: adv_2d_tri_1
166830602db0SMatthew G. Knepley       requires: triangle
166930602db0SMatthew G. Knepley       TODO: how did this ever get in main when there is no support for this
167030602db0SMatthew 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
167130602db0SMatthew G. Knepley 
167230602db0SMatthew G. Knepley     test:
167330602db0SMatthew G. Knepley       suffix: tut_1
167430602db0SMatthew G. Knepley       requires: exodusii
167530602db0SMatthew G. Knepley       nsize: 1
167646ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -
167730602db0SMatthew G. Knepley 
167830602db0SMatthew G. Knepley     test:
167930602db0SMatthew G. Knepley       suffix: tut_2
168030602db0SMatthew G. Knepley       requires: exodusii
168130602db0SMatthew G. Knepley       nsize: 1
168246ac1a18SMatthew G. Knepley       args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw
168330602db0SMatthew G. Knepley 
168430602db0SMatthew G. Knepley     test:
168530602db0SMatthew G. Knepley       suffix: tut_3
168630602db0SMatthew G. Knepley       requires: exodusii
168730602db0SMatthew G. Knepley       nsize: 4
1688e600fa54SMatthew 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
168930602db0SMatthew G. Knepley 
169030602db0SMatthew G. Knepley     test:
169130602db0SMatthew G. Knepley       suffix: tut_4
169230602db0SMatthew G. Knepley       requires: exodusii
169330602db0SMatthew G. Knepley       nsize: 4
1694e600fa54SMatthew 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
169530602db0SMatthew G. Knepley 
169630602db0SMatthew G. Knepley   testset:
169730602db0SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
169830602db0SMatthew G. Knepley 
1699c4762a1bSJed Brown     # 2D Advection 0-10
1700c4762a1bSJed Brown     test:
1701c4762a1bSJed Brown       suffix: 0
1702c4762a1bSJed Brown       requires: exodusii
170346ac1a18SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1704c4762a1bSJed Brown 
1705c4762a1bSJed Brown     test:
1706c4762a1bSJed Brown       suffix: 1
1707c4762a1bSJed Brown       requires: exodusii
170830602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1709c4762a1bSJed Brown 
1710c4762a1bSJed Brown     test:
1711c4762a1bSJed Brown       suffix: 2
1712c4762a1bSJed Brown       requires: exodusii
1713c4762a1bSJed Brown       nsize: 2
171446ac1a18SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1715c4762a1bSJed Brown 
1716c4762a1bSJed Brown     test:
1717c4762a1bSJed Brown       suffix: 3
1718c4762a1bSJed Brown       requires: exodusii
1719c4762a1bSJed Brown       nsize: 2
1720e600fa54SMatthew G. Knepley       args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1721c4762a1bSJed Brown 
1722c4762a1bSJed Brown     test:
1723c4762a1bSJed Brown       suffix: 4
1724c4762a1bSJed Brown       requires: exodusii
17256c2b77d5SStefano Zampini       nsize: 4
17266c2b77d5SStefano 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
1727c4762a1bSJed Brown 
1728c4762a1bSJed Brown     test:
1729c4762a1bSJed Brown       suffix: 5
1730c4762a1bSJed Brown       requires: exodusii
173146ac1a18SMatthew 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
1732c4762a1bSJed Brown 
1733c4762a1bSJed Brown     test:
1734c4762a1bSJed Brown       suffix: 7
1735c4762a1bSJed Brown       requires: exodusii
173630602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1737c4762a1bSJed Brown 
1738c4762a1bSJed Brown     test:
1739c4762a1bSJed Brown       suffix: 8
1740c4762a1bSJed Brown       requires: exodusii
1741c4762a1bSJed Brown       nsize: 2
1742e600fa54SMatthew 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
1743c4762a1bSJed Brown 
1744c4762a1bSJed Brown     test:
1745c4762a1bSJed Brown       suffix: 9
1746c4762a1bSJed Brown       requires: exodusii
1747c4762a1bSJed Brown       nsize: 8
1748e600fa54SMatthew 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
1749c4762a1bSJed Brown 
1750c4762a1bSJed Brown     test:
1751c4762a1bSJed Brown       suffix: 10
1752c4762a1bSJed Brown       requires: exodusii
175330602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo
1754c4762a1bSJed Brown 
1755c4762a1bSJed Brown   # 2D Shallow water
1756993a5519SMatthew G. Knepley   testset:
1757993a5519SMatthew G. Knepley     args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0
1758993a5519SMatthew G. Knepley 
1759c4762a1bSJed Brown     test:
1760c4762a1bSJed Brown       suffix: sw_0
1761c4762a1bSJed Brown       requires: exodusii
1762993a5519SMatthew G. Knepley       args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1763993a5519SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1764993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1765993a5519SMatthew G. Knepley             -monitor height,energy
1766993a5519SMatthew G. Knepley 
1767993a5519SMatthew G. Knepley     test:
1768de8de29fSMatthew G. Knepley       suffix: sw_ceed
1769de8de29fSMatthew G. Knepley       requires: exodusii libceed
1770de8de29fSMatthew G. Knepley       args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1771de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1772de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1773de8de29fSMatthew G. Knepley             -monitor height,energy
1774de8de29fSMatthew G. Knepley 
1775de8de29fSMatthew G. Knepley     test:
1776993a5519SMatthew G. Knepley       suffix: sw_1
1777993a5519SMatthew G. Knepley       nsize: 2
1778993a5519SMatthew G. Knepley       args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1779993a5519SMatthew 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 \
1780993a5519SMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1781993a5519SMatthew G. Knepley             -monitor height,energy
1782c4762a1bSJed Brown 
17838d2c18caSMukkund Sunjii     test:
17848d2c18caSMukkund Sunjii       suffix: sw_hll
1785993a5519SMatthew G. Knepley       args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1786993a5519SMatthew G. Knepley             -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1787993a5519SMatthew G. Knepley             -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1788993a5519SMatthew G. Knepley             -monitor height,energy
1789993a5519SMatthew G. Knepley 
1790de8de29fSMatthew G. Knepley   # 2D Euler
1791de8de29fSMatthew G. Knepley   testset:
1792de8de29fSMatthew G. Knepley     args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1793de8de29fSMatthew G. Knepley           -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy
1794de8de29fSMatthew G. Knepley 
1795de8de29fSMatthew G. Knepley     test:
1796277e32beSMatthew G. Knepley       suffix: euler_0
1797277e32beSMatthew G. Knepley       requires: exodusii !complex
1798277e32beSMatthew G. Knepley       args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1799277e32beSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1800277e32beSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1801277e32beSMatthew G. Knepley 
1802277e32beSMatthew G. Knepley     test:
1803de8de29fSMatthew G. Knepley       suffix: euler_ceed
1804de8de29fSMatthew G. Knepley       requires: exodusii libceed
1805de8de29fSMatthew G. Knepley       args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1806de8de29fSMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1807de8de29fSMatthew G. Knepley             -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1808de8de29fSMatthew G. Knepley 
1809993a5519SMatthew G. Knepley   testset:
1810993a5519SMatthew G. Knepley     args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
18118d2c18caSMukkund Sunjii 
1812c4762a1bSJed Brown     # 2D Advection: p4est
1813c4762a1bSJed Brown     test:
1814c4762a1bSJed Brown       suffix: p4est_advec_2d
1815c4762a1bSJed Brown       requires: p4est
181630602db0SMatthew 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
1817c4762a1bSJed Brown 
1818c4762a1bSJed Brown     # Advection in a box
1819c4762a1bSJed Brown     test:
1820c4762a1bSJed Brown       suffix: adv_2d_quad_0
1821c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1822c4762a1bSJed Brown 
1823c4762a1bSJed Brown     test:
1824c4762a1bSJed Brown       suffix: adv_2d_quad_1
1825c4762a1bSJed 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
1826c4762a1bSJed Brown       timeoutfactor: 3
1827c4762a1bSJed Brown 
1828c4762a1bSJed Brown     test:
1829c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_0
1830c4762a1bSJed Brown       requires: p4est
1831c4762a1bSJed Brown       args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1832c4762a1bSJed Brown 
1833c4762a1bSJed Brown     test:
1834c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_1
1835c4762a1bSJed Brown       requires: p4est
1836c4762a1bSJed 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
1837c4762a1bSJed Brown       timeoutfactor: 3
1838c4762a1bSJed Brown 
1839*53df731dSPierre Jolivet     test: # broken for quad precision
1840c4762a1bSJed Brown       suffix: adv_2d_quad_p4est_adapt_0
1841*53df731dSPierre Jolivet       requires: p4est !__float128
1842c4762a1bSJed 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
1843c4762a1bSJed Brown       timeoutfactor: 3
1844c4762a1bSJed Brown 
1845c4762a1bSJed Brown     test:
1846c4762a1bSJed Brown       suffix: adv_0
1847c4762a1bSJed Brown       requires: exodusii
184846ac1a18SMatthew 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
1849c4762a1bSJed Brown 
1850c4762a1bSJed Brown     test:
1851c4762a1bSJed Brown       suffix: shock_0
1852c4762a1bSJed Brown       requires: p4est !single !complex
185330602db0SMatthew G. Knepley       args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
185430602db0SMatthew 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 \
185530602db0SMatthew 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 \
185630602db0SMatthew 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. \
185730602db0SMatthew G. Knepley       -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
185830602db0SMatthew G. Knepley       -ts_max_time 0.5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
185930602db0SMatthew G. Knepley       -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
1860c4762a1bSJed Brown       timeoutfactor: 3
1861c4762a1bSJed Brown 
1862c4762a1bSJed Brown     # Test GLVis visualization of PetscFV fields
1863c4762a1bSJed Brown     test:
1864c4762a1bSJed Brown       suffix: glvis_adv_2d_tet
186530602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
186630602db0SMatthew G. Knepley             -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
186730602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1868c4762a1bSJed Brown 
1869c4762a1bSJed Brown     test:
1870c4762a1bSJed Brown       suffix: glvis_adv_2d_quad
187130602db0SMatthew G. Knepley       args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
187230602db0SMatthew G. Knepley             -dm_refine 5 -dm_plex_separate_marker \
187330602db0SMatthew G. Knepley             -ts_monitor_solution glvis: -ts_max_steps 0
1874c4762a1bSJed Brown 
1875c4762a1bSJed Brown TEST*/
1876