Actual source code: ex11.c
1: static char help[] = "Second Order TVD Finite Volume Example.\n";
2: /*F
4: We use a second order TVD finite volume method to evolve a system of PDEs. Our simple upwinded residual evaluation loops
5: over all mesh faces and uses a Riemann solver to produce the flux given the face geometry and cell values,
6: \begin{equation}
7: f_i = \mathrm{riemann}(\mathrm{phys}, p_\mathrm{centroid}, \hat n, x^L, x^R)
8: \end{equation}
9: and then update the cell values given the cell volume.
10: \begin{eqnarray}
11: f^L_i &-=& \frac{f_i}{vol^L} \\
12: f^R_i &+=& \frac{f_i}{vol^R}
13: \end{eqnarray}
15: As an example, we can consider the shallow water wave equation,
16: \begin{eqnarray}
17: h_t + \nabla\cdot \left( uh \right) &=& 0 \\
18: (uh)_t + \nabla\cdot \left( u\otimes uh + \frac{g h^2}{2} I \right) &=& 0
19: \end{eqnarray}
20: where $h$ is wave height, $u$ is wave velocity, and $g$ is the acceleration due to gravity.
22: A representative Riemann solver for the shallow water equations is given in the PhysicsRiemann_SW() function,
23: \begin{eqnarray}
24: f^{L,R}_h &=& uh^{L,R} \cdot \hat n \\
25: f^{L,R}_{uh} &=& \frac{f^{L,R}_h}{h^{L,R}} uh^{L,R} + g (h^{L,R})^2 \hat n \\
26: c^{L,R} &=& \sqrt{g h^{L,R}} \\
27: s &=& \max\left( \left|\frac{uh^L \cdot \hat n}{h^L}\right| + c^L, \left|\frac{uh^R \cdot \hat n}{h^R}\right| + c^R \right) \\
28: f_i &=& \frac{A_\mathrm{face}}{2} \left( f^L_i + f^R_i + s \left( x^L_i - x^R_i \right) \right)
29: \end{eqnarray}
30: where $c$ is the local gravity wave speed and $f_i$ is a Rusanov flux.
32: The more sophisticated residual evaluation in RHSFunctionLocal_LS() uses a least-squares fit to a quadratic polynomial
33: over a neighborhood of the given element.
35: The mesh is read in from an ExodusII file, usually generated by Cubit.
37: The example also shows how to handle AMR in a time-dependent TS solver.
38: F*/
39: #include <petscdmplex.h>
40: #include <petscdmforest.h>
41: #include <petscds.h>
42: #include <petscts.h>
44: #include "ex11.h"
46: static PetscFunctionList PhysicsList, PhysicsRiemannList_SW, PhysicsRiemannList_Euler;
48: /* 'User' implements a discretization of a continuous model. */
49: typedef struct _n_User *User;
50: typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
51: typedef PetscErrorCode (*SetUpBCFunction)(DM, PetscDS, Physics);
52: typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
53: typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
54: static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
55: static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
56: static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);
58: typedef struct _n_FunctionalLink *FunctionalLink;
59: struct _n_FunctionalLink {
60: char *name;
61: FunctionalFunction func;
62: void *ctx;
63: PetscInt offset;
64: FunctionalLink next;
65: };
67: struct _n_Model {
68: MPI_Comm comm; /* Does not do collective communication, but some error conditions can be collective */
69: Physics physics;
70: FunctionalLink functionalRegistry;
71: PetscInt maxComputed;
72: PetscInt numMonitored;
73: FunctionalLink *functionalMonitored;
74: PetscInt numCall;
75: FunctionalLink *functionalCall;
76: SolutionFunction solution;
77: SetUpBCFunction setupbc;
78: void *solutionctx;
79: PetscReal maxspeed; /* estimate of global maximum speed (for CFL calculation) */
80: PetscReal bounds[2 * DIM];
81: PetscErrorCode (*errorIndicator)(PetscInt, PetscReal, PetscInt, const PetscScalar[], const PetscScalar[], PetscReal *, void *);
82: void *errorCtx;
83: PetscErrorCode (*setupCEED)(DM, Physics);
84: };
86: struct _n_User {
87: PetscInt vtkInterval; /* For monitor */
88: char outputBasename[PETSC_MAX_PATH_LEN]; /* Basename for output files */
89: PetscInt monitorStepOffset;
90: Model model;
91: PetscBool vtkmon;
92: };
94: #ifdef PETSC_HAVE_LIBCEED
95: // Free a plain data context that was allocated using PETSc; returning libCEED error codes
96: static int FreeContextPetsc(void *data)
97: {
98: if (PetscFree(data)) return CeedError(NULL, CEED_ERROR_ACCESS, "PetscFree failed");
99: return CEED_ERROR_SUCCESS;
100: }
101: #endif
103: /******************* Advect ********************/
104: typedef enum {
105: ADVECT_SOL_TILTED,
106: ADVECT_SOL_BUMP,
107: ADVECT_SOL_BUMP_CAVITY
108: } AdvectSolType;
109: static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "BUMP_CAVITY", "AdvectSolType", "ADVECT_SOL_", 0};
110: typedef enum {
111: ADVECT_SOL_BUMP_CONE,
112: ADVECT_SOL_BUMP_COS
113: } AdvectSolBumpType;
114: static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};
116: typedef struct {
117: PetscReal wind[DIM];
118: } Physics_Advect_Tilted;
119: typedef struct {
120: PetscReal center[DIM];
121: PetscReal radius;
122: AdvectSolBumpType type;
123: } Physics_Advect_Bump;
125: typedef struct {
126: PetscReal inflowState;
127: AdvectSolType soltype;
128: union
129: {
130: Physics_Advect_Tilted tilted;
131: Physics_Advect_Bump bump;
132: } sol;
133: struct {
134: PetscInt Solution;
135: PetscInt Error;
136: } functional;
137: } Physics_Advect;
139: static const struct FieldDescription PhysicsFields_Advect[] = {
140: {"U", 1},
141: {NULL, 0}
142: };
144: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, PetscCtx ctx)
145: {
146: Physics phys = (Physics)ctx;
147: Physics_Advect *advect = (Physics_Advect *)phys->data;
149: PetscFunctionBeginUser;
150: xG[0] = advect->inflowState;
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, PetscCtx ctx)
155: {
156: PetscFunctionBeginUser;
157: xG[0] = xI[0];
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: 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)
162: {
163: Physics_Advect *advect = (Physics_Advect *)phys->data;
164: PetscReal wind[DIM], wn;
166: switch (advect->soltype) {
167: case ADVECT_SOL_TILTED: {
168: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
169: wind[0] = tilted->wind[0];
170: wind[1] = tilted->wind[1];
171: } break;
172: case ADVECT_SOL_BUMP:
173: wind[0] = -qp[1];
174: wind[1] = qp[0];
175: break;
176: case ADVECT_SOL_BUMP_CAVITY: {
177: PetscInt i;
178: PetscReal comp2[3] = {0., 0., 0.}, rad2;
180: rad2 = 0.;
181: for (i = 0; i < dim; i++) {
182: comp2[i] = qp[i] * qp[i];
183: rad2 += comp2[i];
184: }
186: wind[0] = -qp[1];
187: wind[1] = qp[0];
188: if (rad2 > 1.) {
189: PetscInt maxI = 0;
190: PetscReal maxComp2 = comp2[0];
192: for (i = 1; i < dim; i++) {
193: if (comp2[i] > maxComp2) {
194: maxI = i;
195: maxComp2 = comp2[i];
196: }
197: }
198: wind[maxI] = 0.;
199: }
200: } break;
201: default: {
202: for (PetscInt i = 0; i < DIM; ++i) wind[i] = 0.0;
203: }
204: /* default: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for solution type %s",AdvectSolBumpTypes[advect->soltype]); */
205: }
206: wn = Dot2Real(wind, n);
207: flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
208: }
210: static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, PetscCtx ctx)
211: {
212: Physics phys = (Physics)ctx;
213: Physics_Advect *advect = (Physics_Advect *)phys->data;
215: PetscFunctionBeginUser;
216: switch (advect->soltype) {
217: case ADVECT_SOL_TILTED: {
218: PetscReal x0[DIM];
219: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
220: Waxpy2Real(-time, tilted->wind, x, x0);
221: if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
222: else u[0] = advect->inflowState;
223: } break;
224: case ADVECT_SOL_BUMP_CAVITY:
225: case ADVECT_SOL_BUMP: {
226: Physics_Advect_Bump *bump = &advect->sol.bump;
227: PetscReal x0[DIM], v[DIM], r, cost, sint;
228: cost = PetscCosReal(time);
229: sint = PetscSinReal(time);
230: x0[0] = cost * x[0] + sint * x[1];
231: x0[1] = -sint * x[0] + cost * x[1];
232: Waxpy2Real(-1, bump->center, x0, v);
233: r = Norm2Real(v);
234: switch (bump->type) {
235: case ADVECT_SOL_BUMP_CONE:
236: u[0] = PetscMax(1 - r / bump->radius, 0);
237: break;
238: case ADVECT_SOL_BUMP_COS:
239: u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
240: break;
241: }
242: } break;
243: default:
244: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
245: }
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscReal *x, const PetscScalar *y, PetscReal *f, PetscCtx ctx)
250: {
251: Physics phys = (Physics)ctx;
252: Physics_Advect *advect = (Physics_Advect *)phys->data;
253: PetscScalar yexact[1] = {0.0};
255: PetscFunctionBeginUser;
256: PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
257: f[advect->functional.Solution] = PetscRealPart(y[0]);
258: f[advect->functional.Error] = PetscAbsScalar(y[0] - yexact[0]);
259: PetscFunctionReturn(PETSC_SUCCESS);
260: }
262: static PetscErrorCode SetUpBC_Advect(DM dm, PetscDS prob, Physics phys)
263: {
264: const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
265: DMLabel label;
267: PetscFunctionBeginUser;
268: /* Register "canned" boundary conditions and defaults for where to apply. */
269: PetscCall(DMGetLabel(dm, "Face Sets", &label));
270: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (PetscVoidFn *)PhysicsBoundary_Advect_Inflow, NULL, phys, NULL));
271: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (PetscVoidFn *)PhysicsBoundary_Advect_Outflow, NULL, phys, NULL));
272: PetscFunctionReturn(PETSC_SUCCESS);
273: }
275: static PetscErrorCode PhysicsCreate_Advect(Model mod, Physics phys, PetscOptionItems PetscOptionsObject)
276: {
277: Physics_Advect *advect;
279: PetscFunctionBeginUser;
280: phys->field_desc = PhysicsFields_Advect;
281: phys->riemann = (PetscRiemannFn *)(PetscVoidFn *)PhysicsRiemann_Advect;
282: PetscCall(PetscNew(&advect));
283: phys->data = advect;
284: mod->setupbc = SetUpBC_Advect;
286: PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
287: {
288: PetscInt two = 2, dof = 1;
289: advect->soltype = ADVECT_SOL_TILTED;
290: PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL));
291: switch (advect->soltype) {
292: case ADVECT_SOL_TILTED: {
293: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
294: two = 2;
295: tilted->wind[0] = 0.0;
296: tilted->wind[1] = 1.0;
297: PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
298: advect->inflowState = -2.0;
299: PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL));
300: phys->maxspeed = Norm2Real(tilted->wind);
301: } break;
302: case ADVECT_SOL_BUMP_CAVITY:
303: case ADVECT_SOL_BUMP: {
304: Physics_Advect_Bump *bump = &advect->sol.bump;
305: two = 2;
306: bump->center[0] = 2.;
307: bump->center[1] = 0.;
308: PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
309: bump->radius = 0.9;
310: PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL));
311: bump->type = ADVECT_SOL_BUMP_CONE;
312: PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL));
313: phys->maxspeed = 3.; /* radius of mesh, kludge */
314: } break;
315: }
316: }
317: PetscOptionsHeadEnd();
318: /* Initial/transient solution with default boundary conditions */
319: PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys));
320: /* Register "canned" functionals */
321: PetscCall(ModelFunctionalRegister(mod, "Solution", &advect->functional.Solution, PhysicsFunctional_Advect, phys));
322: PetscCall(ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /******************* Shallow Water ********************/
327: static const struct FieldDescription PhysicsFields_SW[] = {
328: {"Height", 1 },
329: {"Momentum", DIM},
330: {NULL, 0 }
331: };
333: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, PetscCtx ctx)
334: {
335: PetscFunctionBeginUser;
336: xG[0] = xI[0];
337: xG[1] = -xI[1];
338: xG[2] = -xI[2];
339: PetscFunctionReturn(PETSC_SUCCESS);
340: }
342: static PetscErrorCode PhysicsSolution_SW(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, PetscCtx ctx)
343: {
344: PetscReal dx[2], r, sigma;
346: PetscFunctionBeginUser;
347: PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
348: dx[0] = x[0] - 1.5;
349: dx[1] = x[1] - 1.0;
350: r = Norm2Real(dx);
351: sigma = 0.5;
352: u[0] = 1 + 2 * PetscExpReal(-PetscSqr(r) / (2 * PetscSqr(sigma)));
353: u[1] = 0.0;
354: u[2] = 0.0;
355: PetscFunctionReturn(PETSC_SUCCESS);
356: }
358: static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, PetscCtx ctx)
359: {
360: Physics phys = (Physics)ctx;
361: Physics_SW *sw = (Physics_SW *)phys->data;
362: const SWNode *x = (const SWNode *)xx;
363: PetscReal u[2];
364: PetscReal h;
366: PetscFunctionBeginUser;
367: h = x->h;
368: Scale2Real(1. / x->h, x->uh, u);
369: f[sw->functional.Height] = h;
370: f[sw->functional.Speed] = Norm2Real(u) + PetscSqrtReal(sw->gravity * h);
371: f[sw->functional.Energy] = 0.5 * (Dot2Real(x->uh, u) + sw->gravity * PetscSqr(h));
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: static PetscErrorCode SetUpBC_SW(DM dm, PetscDS prob, Physics phys)
376: {
377: const PetscInt wallids[] = {100, 101, 200, 300};
378: DMLabel label;
380: PetscFunctionBeginUser;
381: PetscCall(DMGetLabel(dm, "Face Sets", &label));
382: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (PetscVoidFn *)PhysicsBoundary_SW_Wall, NULL, phys, NULL));
383: PetscFunctionReturn(PETSC_SUCCESS);
384: }
386: #ifdef PETSC_HAVE_LIBCEED
387: static PetscErrorCode CreateQFunctionContext_SW(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
388: {
389: Physics_SW *in = (Physics_SW *)phys->data;
390: Physics_SW *sw;
392: PetscFunctionBeginUser;
393: PetscCall(PetscCalloc1(1, &sw));
395: sw->gravity = in->gravity;
397: PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
398: PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*sw), sw));
399: PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
400: PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gravity", offsetof(Physics_SW, gravity), 1, "Acceleration due to gravity"));
401: PetscFunctionReturn(PETSC_SUCCESS);
402: }
403: #endif
405: static PetscErrorCode SetupCEED_SW(DM dm, Physics physics)
406: {
407: #ifdef PETSC_HAVE_LIBCEED
408: Ceed ceed;
409: CeedQFunctionContext qfCtx;
410: #endif
412: PetscFunctionBegin;
413: #ifdef PETSC_HAVE_LIBCEED
414: PetscCall(DMGetCeed(dm, &ceed));
415: PetscCall(CreateQFunctionContext_SW(physics, ceed, &qfCtx));
416: PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_SW_Rusanov_CEED, PhysicsRiemann_SW_Rusanov_CEED_loc, qfCtx));
417: PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
418: #endif
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode PhysicsCreate_SW(Model mod, Physics phys, PetscOptionItems PetscOptionsObject)
423: {
424: Physics_SW *sw;
425: char sw_riemann[64] = "rusanov";
427: PetscFunctionBeginUser;
428: phys->field_desc = PhysicsFields_SW;
429: PetscCall(PetscNew(&sw));
430: phys->data = sw;
431: mod->setupbc = SetUpBC_SW;
432: mod->setupCEED = SetupCEED_SW;
434: PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov", PhysicsRiemann_SW_Rusanov));
435: PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "hll", PhysicsRiemann_SW_HLL));
436: #ifdef PETSC_HAVE_LIBCEED
437: PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_SW, "rusanov_ceed", PhysicsRiemann_SW_Rusanov_CEED));
438: #endif
440: PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
441: {
442: void (*PhysicsRiemann_SW)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
443: sw->gravity = 1.0;
444: PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
445: PetscCall(PetscOptionsFList("-sw_riemann", "Riemann solver", "", PhysicsRiemannList_SW, sw_riemann, sw_riemann, sizeof sw_riemann, NULL));
446: PetscCall(PetscFunctionListFind(PhysicsRiemannList_SW, sw_riemann, &PhysicsRiemann_SW));
447: phys->riemann = (PetscRiemannFn *)(PetscVoidFn *)PhysicsRiemann_SW;
448: }
449: PetscOptionsHeadEnd();
450: phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */
452: PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
453: PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
454: PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
455: PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
456: PetscFunctionReturn(PETSC_SUCCESS);
457: }
459: /******************* Euler Density Shock (EULER_IV_SHOCK,EULER_SS_SHOCK) ********************/
460: /* An initial-value and self-similar solutions of the compressible Euler equations */
461: /* Ravi Samtaney and D. I. Pullin */
462: /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
463: static const struct FieldDescription PhysicsFields_Euler[] = {
464: {"Density", 1 },
465: {"Momentum", DIM},
466: {"Energy", 1 },
467: {NULL, 0 }
468: };
470: /* initial condition */
471: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx);
472: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, PetscCtx ctx)
473: {
474: PetscInt i;
475: Physics phys = (Physics)ctx;
476: Physics_Euler *eu = (Physics_Euler *)phys->data;
477: EulerNode *uu = (EulerNode *)u;
478: PetscReal p0, gamma, c = 0.;
480: PetscFunctionBeginUser;
481: PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
483: for (i = 0; i < DIM; i++) uu->ru[i] = 0.0; /* zero out initial velocity */
484: /* set E and rho */
485: gamma = eu->gamma;
487: if (eu->type == EULER_IV_SHOCK || eu->type == EULER_SS_SHOCK) {
488: /******************* Euler Density Shock ********************/
489: /* On initial-value and self-similar solutions of the compressible Euler equations */
490: /* Ravi Samtaney and D. I. Pullin */
491: /* Phys. Fluids 8, 2650 (1996); http://dx.doi.org/10.1063/1.869050 */
492: /* initial conditions 1: left of shock, 0: left of discontinuity 2: right of discontinuity, */
493: p0 = 1.;
494: if (x[0] < 0.0 + x[1] * eu->itana) {
495: if (x[0] < mod->bounds[0] * 0.5) { /* left of shock (1) */
496: PetscReal amach, rho, press, gas1, p1;
497: amach = eu->amach;
498: rho = 1.;
499: press = p0;
500: p1 = press * (1.0 + 2.0 * gamma / (gamma + 1.0) * (amach * amach - 1.0));
501: gas1 = (gamma - 1.0) / (gamma + 1.0);
502: uu->r = rho * (p1 / press + gas1) / (gas1 * p1 / press + 1.0);
503: uu->ru[0] = ((uu->r - rho) * PetscSqrtReal(gamma * press / rho) * amach);
504: uu->E = p1 / (gamma - 1.0) + .5 / uu->r * uu->ru[0] * uu->ru[0];
505: } else { /* left of discontinuity (0) */
506: uu->r = 1.; /* rho = 1 */
507: uu->E = p0 / (gamma - 1.0);
508: }
509: } else { /* right of discontinuity (2) */
510: uu->r = eu->rhoR;
511: uu->E = p0 / (gamma - 1.0);
512: }
513: } else if (eu->type == EULER_SHOCK_TUBE) {
514: /* 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. */
515: if (x[0] < 0.0) {
516: uu->r = 8.;
517: uu->E = 10. / (gamma - 1.);
518: } else {
519: uu->r = 1.;
520: uu->E = 1. / (gamma - 1.);
521: }
522: } else if (eu->type == EULER_LINEAR_WAVE) {
523: initLinearWave(uu, gamma, x, mod->bounds[1] - mod->bounds[0]);
524: } else SETERRQ(mod->comm, PETSC_ERR_SUP, "Unknown type %d", eu->type);
526: /* set phys->maxspeed: (mod->maxspeed = phys->maxspeed) in main; */
527: PetscCall(SpeedOfSound_PG(gamma, uu, &c));
528: c = (uu->ru[0] / uu->r) + c;
529: if (c > phys->maxspeed) phys->maxspeed = c;
530: PetscFunctionReturn(PETSC_SUCCESS);
531: }
533: /* PetscReal* => EulerNode* conversion */
534: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *a_xI, PetscScalar *a_xG, PetscCtx ctx)
535: {
536: PetscInt i;
537: const EulerNode *xI = (const EulerNode *)a_xI;
538: EulerNode *xG = (EulerNode *)a_xG;
539: Physics phys = (Physics)ctx;
540: Physics_Euler *eu = (Physics_Euler *)phys->data;
542: PetscFunctionBeginUser;
543: xG->r = xI->r; /* ghost cell density - same */
544: xG->E = xI->E; /* ghost cell energy - same */
545: if (n[1] != 0.) { /* top and bottom */
546: xG->ru[0] = xI->ru[0]; /* copy tang to wall */
547: xG->ru[1] = -xI->ru[1]; /* reflect perp to t/b wall */
548: } else { /* sides */
549: for (i = 0; i < DIM; i++) xG->ru[i] = xI->ru[i]; /* copy */
550: }
551: if (eu->type == EULER_LINEAR_WAVE) { /* debug */
552: #if 0
553: PetscPrintf(PETSC_COMM_WORLD,"%s coord=%g,%g\n",PETSC_FUNCTION_NAME,(double)c[0],(double)c[1]);
554: #endif
555: }
556: PetscFunctionReturn(PETSC_SUCCESS);
557: }
559: static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, PetscCtx ctx)
560: {
561: Physics phys = (Physics)ctx;
562: Physics_Euler *eu = (Physics_Euler *)phys->data;
563: const EulerNode *x = (const EulerNode *)xx;
564: PetscReal p;
566: PetscFunctionBeginUser;
567: f[eu->monitor.Density] = x->r;
568: f[eu->monitor.Momentum] = NormDIM(x->ru);
569: f[eu->monitor.Energy] = x->E;
570: f[eu->monitor.Speed] = NormDIM(x->ru) / x->r;
571: PetscCall(Pressure_PG(eu->gamma, x, &p));
572: f[eu->monitor.Pressure] = p;
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: static PetscErrorCode SetUpBC_Euler(DM dm, PetscDS prob, Physics phys)
577: {
578: Physics_Euler *eu = (Physics_Euler *)phys->data;
579: DMLabel label;
581: PetscFunctionBeginUser;
582: PetscCall(DMGetLabel(dm, "Face Sets", &label));
583: if (eu->type == EULER_LINEAR_WAVE) {
584: const PetscInt wallids[] = {100, 101};
585: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (PetscVoidFn *)PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
586: } else {
587: const PetscInt wallids[] = {100, 101, 200, 300};
588: PetscCall(PetscDSAddBoundary(prob, DM_BC_NATURAL_RIEMANN, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (PetscVoidFn *)PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
589: }
590: PetscFunctionReturn(PETSC_SUCCESS);
591: }
593: #ifdef PETSC_HAVE_LIBCEED
594: static PetscErrorCode CreateQFunctionContext_Euler(Physics phys, Ceed ceed, CeedQFunctionContext *qfCtx)
595: {
596: Physics_Euler *in = (Physics_Euler *)phys->data;
597: Physics_Euler *eu;
599: PetscFunctionBeginUser;
600: PetscCall(PetscCalloc1(1, &eu));
602: eu->gamma = in->gamma;
604: PetscCallCEED(CeedQFunctionContextCreate(ceed, qfCtx));
605: PetscCallCEED(CeedQFunctionContextSetData(*qfCtx, CEED_MEM_HOST, CEED_USE_POINTER, sizeof(*eu), eu));
606: PetscCallCEED(CeedQFunctionContextSetDataDestroy(*qfCtx, CEED_MEM_HOST, FreeContextPetsc));
607: PetscCallCEED(CeedQFunctionContextRegisterDouble(*qfCtx, "gamma", offsetof(Physics_Euler, gamma), 1, "Heat capacity ratio"));
608: PetscFunctionReturn(PETSC_SUCCESS);
609: }
610: #endif
612: static PetscErrorCode SetupCEED_Euler(DM dm, Physics physics)
613: {
614: #ifdef PETSC_HAVE_LIBCEED
615: Ceed ceed;
616: CeedQFunctionContext qfCtx;
617: #endif
619: PetscFunctionBegin;
620: #ifdef PETSC_HAVE_LIBCEED
621: PetscCall(DMGetCeed(dm, &ceed));
622: PetscCall(CreateQFunctionContext_Euler(physics, ceed, &qfCtx));
623: PetscCall(DMCeedCreateFVM(dm, PETSC_TRUE, PhysicsRiemann_Euler_Godunov_CEED, PhysicsRiemann_Euler_Godunov_CEED_loc, qfCtx));
624: PetscCallCEED(CeedQFunctionContextDestroy(&qfCtx));
625: #endif
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: static PetscErrorCode PhysicsCreate_Euler(Model mod, Physics phys, PetscOptionItems PetscOptionsObject)
630: {
631: Physics_Euler *eu;
633: PetscFunctionBeginUser;
634: phys->field_desc = PhysicsFields_Euler;
635: phys->riemann = (PetscRiemannFn *)(PetscVoidFn *)PhysicsRiemann_Euler_Godunov;
636: PetscCall(PetscNew(&eu));
637: phys->data = eu;
638: mod->setupbc = SetUpBC_Euler;
639: mod->setupCEED = SetupCEED_Euler;
641: PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov", PhysicsRiemann_Euler_Godunov));
642: #ifdef PETSC_HAVE_LIBCEED
643: PetscCall(PetscFunctionListAdd(&PhysicsRiemannList_Euler, "godunov_ceed", PhysicsRiemann_Euler_Godunov_CEED));
644: #endif
646: PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
647: {
648: void (*PhysicsRiemann_Euler)(PetscInt, PetscInt, const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscInt, const PetscScalar, PetscScalar *, Physics);
649: PetscReal alpha;
650: char type[64] = "linear_wave";
651: char eu_riemann[64] = "godunov";
652: PetscBool is;
653: eu->gamma = 1.4;
654: eu->amach = 2.02;
655: eu->rhoR = 3.0;
656: eu->itana = 0.57735026918963; /* angle of Euler self similar (SS) shock */
657: PetscCall(PetscOptionsFList("-eu_riemann", "Riemann solver", "", PhysicsRiemannList_Euler, eu_riemann, eu_riemann, sizeof eu_riemann, NULL));
658: PetscCall(PetscFunctionListFind(PhysicsRiemannList_Euler, eu_riemann, &PhysicsRiemann_Euler));
659: phys->riemann = (PetscRiemannFn *)(PetscVoidFn *)PhysicsRiemann_Euler;
660: PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->gamma, &eu->gamma, NULL));
661: PetscCall(PetscOptionsReal("-eu_amach", "Shock speed (Mach)", "", eu->amach, &eu->amach, NULL));
662: PetscCall(PetscOptionsReal("-eu_rho2", "Density right of discontinuity", "", eu->rhoR, &eu->rhoR, NULL));
663: alpha = 60.;
664: PetscCall(PetscOptionsRangeReal("-eu_alpha", "Angle of discontinuity", "", alpha, &alpha, NULL, 0.0, 90.0));
665: eu->itana = 1. / PetscTanReal(alpha * PETSC_PI / 180.0);
666: PetscCall(PetscOptionsString("-eu_type", "Type of Euler test", "", type, type, sizeof(type), NULL));
667: PetscCall(PetscStrcmp(type, "linear_wave", &is));
668: if (is) {
669: /* Remember this should be periodic */
670: eu->type = EULER_LINEAR_WAVE;
671: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "linear_wave"));
672: } else {
673: PetscCheck(DIM == 2, PETSC_COMM_WORLD, PETSC_ERR_SUP, "DIM must be 2 unless linear wave test %s", type);
674: PetscCall(PetscStrcmp(type, "iv_shock", &is));
675: if (is) {
676: eu->type = EULER_IV_SHOCK;
677: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "iv_shock"));
678: } else {
679: PetscCall(PetscStrcmp(type, "ss_shock", &is));
680: if (is) {
681: eu->type = EULER_SS_SHOCK;
682: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "ss_shock"));
683: } else {
684: PetscCall(PetscStrcmp(type, "shock_tube", &is));
685: PetscCheck(is, PETSC_COMM_WORLD, PETSC_ERR_SUP, "Unknown Euler type %s", type);
686: eu->type = EULER_SHOCK_TUBE;
687: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s set Euler type: %s\n", PETSC_FUNCTION_NAME, "shock_tube"));
688: }
689: }
690: }
691: }
692: PetscOptionsHeadEnd();
693: phys->maxspeed = 0.; /* will get set in solution */
694: PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
695: PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
696: PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
697: PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
698: PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
699: PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
700: PetscFunctionReturn(PETSC_SUCCESS);
701: }
703: static PetscErrorCode ErrorIndicator_Simple(PetscInt dim, PetscReal volume, PetscInt numComps, const PetscScalar u[], const PetscScalar grad[], PetscReal *error, PetscCtx ctx)
704: {
705: PetscReal err = 0.;
707: PetscFunctionBeginUser;
708: for (PetscInt i = 0; i < numComps; i++) {
709: for (PetscInt j = 0; j < dim; j++) err += PetscSqr(PetscRealPart(grad[i * dim + j]));
710: }
711: *error = volume * err;
712: PetscFunctionReturn(PETSC_SUCCESS);
713: }
715: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
716: {
717: PetscSF sfPoint;
718: PetscSection coordSection;
719: Vec coordinates;
720: PetscSection sectionCell;
721: PetscScalar *part;
722: PetscInt cStart, cEnd, c;
723: PetscMPIInt rank;
725: PetscFunctionBeginUser;
726: PetscCall(DMGetCoordinateSection(dm, &coordSection));
727: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
728: PetscCall(DMClone(dm, dmCell));
729: PetscCall(DMGetPointSF(dm, &sfPoint));
730: PetscCall(DMSetPointSF(*dmCell, sfPoint));
731: PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
732: PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
733: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
734: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell));
735: PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
736: PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
737: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
738: PetscCall(PetscSectionSetUp(sectionCell));
739: PetscCall(DMSetLocalSection(*dmCell, sectionCell));
740: PetscCall(PetscSectionDestroy(§ionCell));
741: PetscCall(DMCreateLocalVector(*dmCell, partition));
742: PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
743: PetscCall(VecGetArray(*partition, &part));
744: for (c = cStart; c < cEnd; ++c) {
745: PetscScalar *p;
747: PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
748: p[0] = rank;
749: }
750: PetscCall(VecRestoreArray(*partition, &part));
751: PetscFunctionReturn(PETSC_SUCCESS);
752: }
754: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
755: {
756: DM plex, dmMass, dmFace, dmCell, dmCoord;
757: PetscSection coordSection;
758: Vec coordinates, facegeom, cellgeom;
759: PetscSection sectionMass;
760: PetscScalar *m;
761: const PetscScalar *fgeom, *cgeom, *coords;
762: PetscInt vStart, vEnd, v;
764: PetscFunctionBeginUser;
765: PetscCall(DMConvert(dm, DMPLEX, &plex));
766: PetscCall(DMGetCoordinateSection(dm, &coordSection));
767: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
768: PetscCall(DMClone(dm, &dmMass));
769: PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
770: PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
771: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionMass));
772: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
773: PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
774: for (v = vStart; v < vEnd; ++v) {
775: PetscInt numFaces;
777: PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
778: PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
779: }
780: PetscCall(PetscSectionSetUp(sectionMass));
781: PetscCall(DMSetLocalSection(dmMass, sectionMass));
782: PetscCall(PetscSectionDestroy(§ionMass));
783: PetscCall(DMGetLocalVector(dmMass, massMatrix));
784: PetscCall(VecGetArray(*massMatrix, &m));
785: PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
786: PetscCall(VecGetDM(facegeom, &dmFace));
787: PetscCall(VecGetArrayRead(facegeom, &fgeom));
788: PetscCall(VecGetDM(cellgeom, &dmCell));
789: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
790: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
791: PetscCall(VecGetArrayRead(coordinates, &coords));
792: for (v = vStart; v < vEnd; ++v) {
793: const PetscInt *faces;
794: PetscFVFaceGeom *fgA, *fgB, *cg;
795: PetscScalar *vertex;
796: PetscInt numFaces, sides[2], f, g;
798: PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
799: PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
800: PetscCall(DMPlexGetSupport(dmMass, v, &faces));
801: for (f = 0; f < numFaces; ++f) {
802: sides[0] = faces[f];
803: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
804: for (g = 0; g < numFaces; ++g) {
805: const PetscInt *cells = NULL;
806: PetscReal area = 0.0;
807: PetscInt numCells;
809: sides[1] = faces[g];
810: PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
811: PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
812: PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
813: PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
814: area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
815: area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
816: m[f * numFaces + g] = Dot2Real(fgA->normal, fgB->normal) * area * 0.5;
817: PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
818: }
819: }
820: }
821: PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
822: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
823: PetscCall(VecRestoreArrayRead(coordinates, &coords));
824: PetscCall(VecRestoreArray(*massMatrix, &m));
825: PetscCall(DMDestroy(&dmMass));
826: PetscCall(DMDestroy(&plex));
827: PetscFunctionReturn(PETSC_SUCCESS);
828: }
830: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
831: static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, PetscCtx ctx)
832: {
833: PetscFunctionBeginUser;
834: mod->solution = func;
835: mod->solutionctx = ctx;
836: PetscFunctionReturn(PETSC_SUCCESS);
837: }
839: static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, PetscCtx ctx)
840: {
841: FunctionalLink link, *ptr;
842: PetscInt lastoffset = -1;
844: PetscFunctionBeginUser;
845: for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
846: PetscCall(PetscNew(&link));
847: PetscCall(PetscStrallocpy(name, &link->name));
848: link->offset = lastoffset + 1;
849: link->func = func;
850: link->ctx = ctx;
851: link->next = NULL;
852: *ptr = link;
853: *offset = link->offset;
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptionItems PetscOptionsObject)
858: {
859: PetscInt i;
860: FunctionalLink link;
861: char *names[256];
863: PetscFunctionBeginUser;
864: mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
865: PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
866: /* Create list of functionals that will be computed somehow */
867: PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
868: /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
869: PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
870: mod->numCall = 0;
871: for (i = 0; i < mod->numMonitored; i++) {
872: for (link = mod->functionalRegistry; link; link = link->next) {
873: PetscBool match;
874: PetscCall(PetscStrcasecmp(names[i], link->name, &match));
875: if (match) break;
876: }
877: PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
878: mod->functionalMonitored[i] = link;
879: for (PetscInt j = 0; j < i; j++) {
880: if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
881: }
882: mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
883: next_name:
884: PetscCall(PetscFree(names[i]));
885: }
887: /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
888: mod->maxComputed = -1;
889: for (link = mod->functionalRegistry; link; link = link->next) {
890: for (i = 0; i < mod->numCall; i++) {
891: FunctionalLink call = mod->functionalCall[i];
892: if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
893: }
894: }
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
899: {
900: FunctionalLink l, next;
902: PetscFunctionBeginUser;
903: if (!link) PetscFunctionReturn(PETSC_SUCCESS);
904: l = *link;
905: *link = NULL;
906: for (; l; l = next) {
907: next = l->next;
908: PetscCall(PetscFree(l->name));
909: PetscCall(PetscFree(l));
910: }
911: PetscFunctionReturn(PETSC_SUCCESS);
912: }
914: /* put the solution callback into a functional callback */
915: static PetscErrorCode SolutionFunctional(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *modctx)
916: {
917: Model mod;
919: PetscFunctionBeginUser;
920: mod = (Model)modctx;
921: PetscCall((*mod->solution)(mod, time, x, u, mod->solutionctx));
922: PetscFunctionReturn(PETSC_SUCCESS);
923: }
925: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
926: {
927: PetscErrorCode (*func[1])(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nf, PetscScalar *u, PetscCtx ctx);
928: PetscCtx ctx[1];
929: Model mod = user->model;
931: PetscFunctionBeginUser;
932: func[0] = SolutionFunctional;
933: ctx[0] = (void *)mod;
934: PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, X));
935: PetscFunctionReturn(PETSC_SUCCESS);
936: }
938: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
939: {
940: PetscFunctionBeginUser;
941: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
942: PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
943: PetscCall(PetscViewerFileSetName(*viewer, filename));
944: PetscFunctionReturn(PETSC_SUCCESS);
945: }
947: static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, PetscCtx ctx)
948: {
949: User user = (User)ctx;
950: DM dm, plex;
951: PetscViewer viewer;
952: char filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
953: PetscReal xnorm;
954: PetscBool rollback;
956: PetscFunctionBeginUser;
957: PetscCall(TSGetStepRollBack(ts, &rollback));
958: if (rollback) PetscFunctionReturn(PETSC_SUCCESS);
959: PetscCall(PetscObjectSetName((PetscObject)X, "u"));
960: PetscCall(VecGetDM(X, &dm));
961: PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));
963: if (stepnum >= 0) stepnum += user->monitorStepOffset;
964: if (stepnum >= 0) { /* No summary for final time */
965: Model mod = user->model;
966: Vec cellgeom;
967: PetscInt c, cStart, cEnd, fcount, i;
968: size_t ftableused, ftablealloc;
969: const PetscScalar *cgeom, *x;
970: DM dmCell;
971: DMLabel vtkLabel;
972: PetscReal *fmin, *fmax, *fintegral, *ftmp;
974: PetscCall(DMConvert(dm, DMPLEX, &plex));
975: PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
976: fcount = mod->maxComputed + 1;
977: PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
978: for (i = 0; i < fcount; i++) {
979: fmin[i] = PETSC_MAX_REAL;
980: fmax[i] = PETSC_MIN_REAL;
981: fintegral[i] = 0;
982: }
983: PetscCall(VecGetDM(cellgeom, &dmCell));
984: PetscCall(DMPlexGetSimplexOrBoxCells(dmCell, 0, &cStart, &cEnd));
985: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
986: PetscCall(VecGetArrayRead(X, &x));
987: PetscCall(DMGetLabel(dm, "vtk", &vtkLabel));
988: for (c = cStart; c < cEnd; ++c) {
989: PetscFVCellGeom *cg;
990: const PetscScalar *cx = NULL;
991: PetscInt vtkVal = 0;
993: /* not that these two routines as currently implemented work for any dm with a
994: * localSection/globalSection */
995: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
996: PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
997: if (vtkLabel) PetscCall(DMLabelGetValue(vtkLabel, c, &vtkVal));
998: if (!vtkVal || !cx) continue; /* ghost, or not a global cell */
999: for (i = 0; i < mod->numCall; i++) {
1000: FunctionalLink flink = mod->functionalCall[i];
1001: PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
1002: }
1003: for (i = 0; i < fcount; i++) {
1004: fmin[i] = PetscMin(fmin[i], ftmp[i]);
1005: fmax[i] = PetscMax(fmax[i], ftmp[i]);
1006: fintegral[i] += cg->volume * ftmp[i];
1007: }
1008: }
1009: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
1010: PetscCall(VecRestoreArrayRead(X, &x));
1011: PetscCall(DMDestroy(&plex));
1012: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmin, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1013: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fmax, (PetscMPIInt)fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1014: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, fintegral, (PetscMPIInt)fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
1016: ftablealloc = fcount * 100;
1017: ftableused = 0;
1018: PetscCall(PetscMalloc1(ftablealloc, &ftable));
1019: for (i = 0; i < mod->numMonitored; i++) {
1020: size_t countused;
1021: char buffer[256], *p;
1022: FunctionalLink flink = mod->functionalMonitored[i];
1023: PetscInt id = flink->offset;
1024: if (i % 3) {
1025: PetscCall(PetscArraycpy(buffer, " ", 2));
1026: p = buffer + 2;
1027: } else if (i) {
1028: char newline[] = "\n";
1029: PetscCall(PetscMemcpy(buffer, newline, sizeof(newline) - 1));
1030: p = buffer + sizeof(newline) - 1;
1031: } else {
1032: p = buffer;
1033: }
1034: 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]));
1035: countused--;
1036: countused += p - buffer;
1037: if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1038: char *ftablenew;
1039: ftablealloc = 2 * ftablealloc + countused;
1040: PetscCall(PetscMalloc(ftablealloc, &ftablenew));
1041: PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
1042: PetscCall(PetscFree(ftable));
1043: ftable = ftablenew;
1044: }
1045: PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
1046: ftableused += countused;
1047: ftable[ftableused] = 0;
1048: }
1049: PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));
1051: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT " time %8.4g |x| %8.4g %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
1052: PetscCall(PetscFree(ftable));
1053: }
1054: if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
1055: if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1056: if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
1057: PetscCall(TSGetStepNumber(ts, &stepnum));
1058: }
1059: PetscCall(PetscSNPrintf(filename, sizeof filename, "%s-%03" PetscInt_FMT ".vtu", user->outputBasename, stepnum));
1060: PetscCall(OutputVTK(dm, filename, &viewer));
1061: PetscCall(VecView(X, viewer));
1062: PetscCall(PetscViewerDestroy(&viewer));
1063: }
1064: PetscFunctionReturn(PETSC_SUCCESS);
1065: }
1067: static PetscErrorCode initializeTS(DM dm, User user, TS *ts)
1068: {
1069: #ifdef PETSC_HAVE_LIBCEED
1070: PetscBool useCeed;
1071: #endif
1073: PetscFunctionBeginUser;
1074: PetscCall(TSCreate(PetscObjectComm((PetscObject)dm), ts));
1075: PetscCall(TSSetType(*ts, TSSSP));
1076: PetscCall(TSSetDM(*ts, dm));
1077: if (user->vtkmon) PetscCall(TSMonitorSet(*ts, MonitorVTK, user, NULL));
1078: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user));
1079: #ifdef PETSC_HAVE_LIBCEED
1080: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1081: if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user));
1082: else
1083: #endif
1084: PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
1085: PetscCall(TSSetMaxTime(*ts, 2.0));
1086: PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
1087: PetscFunctionReturn(PETSC_SUCCESS);
1088: }
1090: typedef struct {
1091: PetscFV fvm;
1092: VecTagger refineTag;
1093: VecTagger coarsenTag;
1094: DM adaptedDM;
1095: User user;
1096: PetscReal cfl;
1097: PetscLimiter limiter;
1098: PetscLimiter noneLimiter;
1099: } TransferCtx;
1101: static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, PetscCtx ctx)
1102: {
1103: TransferCtx *tctx = (TransferCtx *)ctx;
1104: PetscFV fvm = tctx->fvm;
1105: VecTagger refineTag = tctx->refineTag;
1106: VecTagger coarsenTag = tctx->coarsenTag;
1107: User user = tctx->user;
1108: DM dm, gradDM, plex, cellDM, adaptedDM = NULL;
1109: Vec cellGeom, faceGeom;
1110: PetscBool computeGradient;
1111: Vec grad, locGrad, locX, errVec;
1112: PetscInt cStart, cEnd, c, dim, nRefine, nCoarsen;
1113: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1114: PetscScalar *errArray;
1115: const PetscScalar *pointVals;
1116: const PetscScalar *pointGrads;
1117: const PetscScalar *pointGeom;
1118: DMLabel adaptLabel = NULL;
1119: IS refineIS, coarsenIS;
1121: PetscFunctionBeginUser;
1122: *resize = PETSC_FALSE;
1123: PetscCall(VecGetDM(sol, &dm));
1124: PetscCall(DMGetDimension(dm, &dim));
1125: PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
1126: PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
1127: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
1128: PetscCall(DMConvert(dm, DMPLEX, &plex));
1129: PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
1130: PetscCall(DMCreateLocalVector(plex, &locX));
1131: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
1132: PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
1133: PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
1134: PetscCall(DMCreateGlobalVector(gradDM, &grad));
1135: PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
1136: PetscCall(DMCreateLocalVector(gradDM, &locGrad));
1137: PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
1138: PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
1139: PetscCall(VecDestroy(&grad));
1140: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
1141: PetscCall(VecGetArrayRead(locGrad, &pointGrads));
1142: PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
1143: PetscCall(VecGetArrayRead(locX, &pointVals));
1144: PetscCall(VecGetDM(cellGeom, &cellDM));
1145: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
1146: PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
1147: PetscCall(VecSetUp(errVec));
1148: PetscCall(VecGetArray(errVec, &errArray));
1149: for (c = cStart; c < cEnd; c++) {
1150: PetscReal errInd = 0.;
1151: PetscScalar *pointGrad;
1152: PetscScalar *pointVal;
1153: PetscFVCellGeom *cg;
1155: PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
1156: PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
1157: PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));
1159: PetscCall((*user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1160: errArray[c - cStart] = errInd;
1161: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
1162: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
1163: }
1164: PetscCall(VecRestoreArray(errVec, &errArray));
1165: PetscCall(VecRestoreArrayRead(locX, &pointVals));
1166: PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
1167: PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
1168: PetscCall(VecDestroy(&locGrad));
1169: PetscCall(VecDestroy(&locX));
1170: PetscCall(DMDestroy(&plex));
1172: PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
1173: PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
1174: PetscCall(ISGetSize(refineIS, &nRefine));
1175: PetscCall(ISGetSize(coarsenIS, &nCoarsen));
1176: if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
1177: if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
1178: PetscCall(ISDestroy(&coarsenIS));
1179: PetscCall(ISDestroy(&refineIS));
1180: PetscCall(VecDestroy(&errVec));
1182: PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
1183: PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1184: minMaxInd[1] = -minMaxInd[1];
1185: PetscCallMPI(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1186: PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1187: if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
1188: PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1189: }
1190: PetscCall(DMLabelDestroy(&adaptLabel));
1191: if (adaptedDM) {
1192: tctx->adaptedDM = adaptedDM;
1193: PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nstep, nRefine, nCoarsen));
1194: *resize = PETSC_TRUE;
1195: } else {
1196: PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1197: }
1198: PetscFunctionReturn(PETSC_SUCCESS);
1199: }
1201: static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], PetscCtx ctx)
1202: {
1203: TransferCtx *tctx = (TransferCtx *)ctx;
1204: DM dm;
1205: PetscReal time;
1207: PetscFunctionBeginUser;
1208: PetscCall(TSGetDM(ts, &dm));
1209: PetscCall(TSGetTime(ts, &time));
1210: PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
1211: for (PetscInt i = 0; i < nv; i++) {
1212: PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
1213: PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
1214: }
1215: PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */
1217: Model mod = tctx->user->model;
1218: Physics phys = mod->physics;
1219: PetscReal minRadius;
1221: PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius));
1222: PetscCallMPI(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1223: PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed");
1225: PetscReal dt = tctx->cfl * minRadius / mod->maxspeed;
1226: PetscCall(TSSetTimeStep(ts, dt));
1228: PetscCall(TSSetDM(ts, tctx->adaptedDM));
1229: PetscCall(DMDestroy(&tctx->adaptedDM));
1230: PetscFunctionReturn(PETSC_SUCCESS);
1231: }
1233: int main(int argc, char **argv)
1234: {
1235: MPI_Comm comm;
1236: PetscDS prob;
1237: PetscFV fvm;
1238: PetscLimiter limiter = NULL, noneLimiter = NULL;
1239: User user;
1240: Model mod;
1241: Physics phys;
1242: DM dm, plex;
1243: PetscReal ftime, cfl, dt, minRadius;
1244: PetscInt dim, nsteps;
1245: TS ts;
1246: TSConvergedReason reason;
1247: Vec X;
1248: PetscViewer viewer;
1249: PetscBool vtkCellGeom, useAMR;
1250: PetscInt adaptInterval;
1251: char physname[256] = "advect";
1252: VecTagger refineTag = NULL, coarsenTag = NULL;
1253: TransferCtx tctx;
1255: PetscFunctionBeginUser;
1256: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1257: comm = PETSC_COMM_WORLD;
1259: PetscCall(PetscNew(&user));
1260: PetscCall(PetscNew(&user->model));
1261: PetscCall(PetscNew(&user->model->physics));
1262: mod = user->model;
1263: phys = mod->physics;
1264: mod->comm = comm;
1265: useAMR = PETSC_FALSE;
1266: adaptInterval = 1;
1268: /* Register physical models to be available on the command line */
1269: PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
1270: PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
1271: PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1273: PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1274: {
1275: cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1276: PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1277: user->vtkInterval = 1;
1278: PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1279: user->vtkmon = PETSC_TRUE;
1280: PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1281: vtkCellGeom = PETSC_FALSE;
1282: PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
1283: PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
1284: PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
1285: PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
1286: PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1287: }
1288: PetscOptionsEnd();
1290: if (useAMR) {
1291: VecTaggerBox refineBox, coarsenBox;
1293: refineBox.min = refineBox.max = PETSC_MAX_REAL;
1294: coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;
1296: PetscCall(VecTaggerCreate(comm, &refineTag));
1297: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
1298: PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
1299: PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
1300: PetscCall(VecTaggerSetFromOptions(refineTag));
1301: PetscCall(VecTaggerSetUp(refineTag));
1302: PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
1304: PetscCall(VecTaggerCreate(comm, &coarsenTag));
1305: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
1306: PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
1307: PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
1308: PetscCall(VecTaggerSetFromOptions(coarsenTag));
1309: PetscCall(VecTaggerSetUp(coarsenTag));
1310: PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1311: }
1313: PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1314: {
1315: PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems);
1316: PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
1317: PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
1318: PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
1319: PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1320: /* Count number of fields and dofs */
1321: for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1322: PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1323: PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1324: }
1325: PetscOptionsEnd();
1327: /* Create mesh */
1328: {
1329: PetscInt i;
1331: PetscCall(DMCreate(comm, &dm));
1332: PetscCall(DMSetType(dm, DMPLEX));
1333: PetscCall(DMSetFromOptions(dm));
1334: for (i = 0; i < DIM; i++) {
1335: mod->bounds[2 * i] = 0.;
1336: mod->bounds[2 * i + 1] = 1.;
1337: }
1338: dim = DIM;
1339: { /* a null name means just do a hex box */
1340: PetscInt cells[3] = {1, 1, 1}, n = 3;
1341: PetscBool flg2, skew = PETSC_FALSE;
1342: PetscInt nret2 = 2 * DIM;
1343: PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
1344: 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));
1345: PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
1346: PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1347: PetscOptionsEnd();
1348: /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1349: if (flg2) {
1350: PetscInt dimEmbed;
1351: PetscInt nCoords;
1352: PetscScalar *coords;
1353: Vec coordinates;
1355: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1356: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1357: PetscCall(VecGetLocalSize(coordinates, &nCoords));
1358: PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
1359: PetscCall(VecGetArray(coordinates, &coords));
1360: for (PetscInt i = 0; i < nCoords; i += dimEmbed) {
1361: PetscScalar *coord = &coords[i];
1362: for (PetscInt j = 0; j < dimEmbed; j++) {
1363: coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1364: if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1365: if (cells[0] == 2 && i == 8) {
1366: coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
1367: } else if (cells[0] == 3) {
1368: if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1369: else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1370: else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1371: }
1372: }
1373: }
1374: }
1375: PetscCall(VecRestoreArray(coordinates, &coords));
1376: PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1377: }
1378: }
1379: }
1380: PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
1381: PetscCall(DMGetDimension(dm, &dim));
1383: /* set up BCs, functions, tags */
1384: PetscCall(DMCreateLabel(dm, "Face Sets"));
1385: mod->errorIndicator = ErrorIndicator_Simple;
1387: {
1388: DM gdm;
1390: PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
1391: PetscCall(DMDestroy(&dm));
1392: dm = gdm;
1393: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1394: }
1396: PetscCall(PetscFVCreate(comm, &fvm));
1397: PetscCall(PetscFVSetFromOptions(fvm));
1398: PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
1399: PetscCall(PetscFVSetSpatialDimension(fvm, dim));
1400: PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1401: {
1402: PetscInt f, dof;
1403: for (f = 0, dof = 0; f < phys->nfields; f++) {
1404: PetscInt newDof = phys->field_desc[f].dof;
1406: if (newDof == 1) PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
1407: else {
1408: for (PetscInt j = 0; j < newDof; j++) {
1409: char compName[256] = "Unknown";
1411: PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
1412: PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1413: }
1414: }
1415: dof += newDof;
1416: }
1417: }
1418: /* FV is now structured with one field having all physics as components */
1419: PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
1420: PetscCall(DMCreateDS(dm));
1421: PetscCall(DMGetDS(dm, &prob));
1422: PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
1423: PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
1424: PetscCall((*mod->setupbc)(dm, prob, phys));
1425: PetscCall(PetscDSSetFromOptions(prob));
1426: {
1427: char convType[256];
1428: PetscBool flg;
1430: PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
1431: PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1432: PetscOptionsEnd();
1433: if (flg) {
1434: DM dmConv;
1436: PetscCall(DMConvert(dm, convType, &dmConv));
1437: if (dmConv) {
1438: PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
1439: PetscCall(DMDestroy(&dm));
1440: dm = dmConv;
1441: PetscCall(DMSetFromOptions(dm));
1442: }
1443: }
1444: }
1445: #ifdef PETSC_HAVE_LIBCEED
1446: {
1447: PetscBool useCeed;
1448: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1449: if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1450: }
1451: #endif
1453: PetscCall(initializeTS(dm, user, &ts));
1455: PetscCall(DMCreateGlobalVector(dm, &X));
1456: PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1457: PetscCall(SetInitialCondition(dm, X, user));
1458: if (useAMR) {
1459: /* use no limiting when reconstructing gradients for adaptivity */
1460: PetscCall(PetscFVGetLimiter(fvm, &limiter));
1461: PetscCall(PetscObjectReference((PetscObject)limiter));
1462: PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
1463: PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
1465: /* Refinement context */
1466: tctx.fvm = fvm;
1467: tctx.refineTag = refineTag;
1468: tctx.coarsenTag = coarsenTag;
1469: tctx.adaptedDM = NULL;
1470: tctx.user = user;
1471: tctx.noneLimiter = noneLimiter;
1472: tctx.limiter = limiter;
1473: tctx.cfl = cfl;
1475: /* Do some initial refinement steps */
1476: for (PetscInt adaptIter = 0;; ++adaptIter) {
1477: PetscLogDouble bytes;
1478: PetscBool resize;
1480: PetscCall(PetscMemoryGetCurrentUsage(&bytes));
1481: PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, bytes));
1482: PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
1483: PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));
1485: PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
1486: if (!resize) break;
1487: PetscCall(DMDestroy(&dm));
1488: PetscCall(VecDestroy(&X));
1489: dm = tctx.adaptedDM;
1490: tctx.adaptedDM = NULL;
1491: PetscCall(TSSetDM(ts, dm));
1492: PetscCall(DMCreateGlobalVector(dm, &X));
1493: PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1494: PetscCall(SetInitialCondition(dm, X, user));
1495: }
1496: }
1498: PetscCall(DMConvert(dm, DMPLEX, &plex));
1499: if (vtkCellGeom) {
1500: DM dmCell;
1501: Vec cellgeom, partition;
1503: PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1504: PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
1505: PetscCall(VecView(cellgeom, viewer));
1506: PetscCall(PetscViewerDestroy(&viewer));
1507: PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
1508: PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
1509: PetscCall(VecView(partition, viewer));
1510: PetscCall(PetscViewerDestroy(&viewer));
1511: PetscCall(VecDestroy(&partition));
1512: PetscCall(DMDestroy(&dmCell));
1513: }
1514: /* collect max maxspeed from all processes -- todo */
1515: PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
1516: PetscCall(DMDestroy(&plex));
1517: PetscCallMPI(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1518: PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1519: dt = cfl * minRadius / mod->maxspeed;
1520: PetscCall(TSSetTimeStep(ts, dt));
1521: PetscCall(TSSetFromOptions(ts));
1523: /* When using adaptive mesh refinement
1524: specify callbacks to refine the solution
1525: and interpolate data from old to new mesh
1526: When mesh adaption is requested, the step will be rolled back
1527: */
1528: if (useAMR) PetscCall(TSSetResize(ts, PETSC_TRUE, adaptToleranceFVMSetUp, Transfer, &tctx));
1529: PetscCall(TSSetSolution(ts, X));
1530: PetscCall(VecDestroy(&X));
1531: PetscCall(TSSolve(ts, NULL));
1532: PetscCall(TSGetSolveTime(ts, &ftime));
1533: PetscCall(TSGetStepNumber(ts, &nsteps));
1534: PetscCall(TSGetConvergedReason(ts, &reason));
1535: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
1536: PetscCall(TSDestroy(&ts));
1538: PetscCall(VecTaggerDestroy(&refineTag));
1539: PetscCall(VecTaggerDestroy(&coarsenTag));
1540: PetscCall(PetscFunctionListDestroy(&PhysicsList));
1541: PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1542: PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
1543: PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
1544: PetscCall(PetscFree(user->model->functionalMonitored));
1545: PetscCall(PetscFree(user->model->functionalCall));
1546: PetscCall(PetscFree(user->model->physics->data));
1547: PetscCall(PetscFree(user->model->physics));
1548: PetscCall(PetscFree(user->model));
1549: PetscCall(PetscFree(user));
1550: PetscCall(PetscLimiterDestroy(&limiter));
1551: PetscCall(PetscLimiterDestroy(&noneLimiter));
1552: PetscCall(PetscFVDestroy(&fvm));
1553: PetscCall(DMDestroy(&dm));
1554: PetscCall(PetscFinalize());
1555: return 0;
1556: }
1558: /* Subroutine to set up the initial conditions for the */
1559: /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1560: /* ----------------------------------------------------------------------- */
1561: int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1562: {
1563: int j, k;
1564: /* Wc=matmul(lv,Ueq) 3 vars */
1565: for (k = 0; k < 3; ++k) {
1566: wc[k] = 0.;
1567: for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1568: }
1569: return 0;
1570: }
1571: /* ----------------------------------------------------------------------- */
1572: int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1573: {
1574: int k, j;
1575: /* V=matmul(rv,WC) 3 vars */
1576: for (k = 0; k < 3; ++k) {
1577: v[k] = 0.;
1578: for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1579: }
1580: return 0;
1581: }
1582: /* ---------------------------------------------------------------------- */
1583: int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1584: {
1585: int j, k;
1586: PetscReal rho, csnd, p0;
1587: /* PetscScalar u; */
1589: for (k = 0; k < 3; ++k)
1590: for (j = 0; j < 3; ++j) {
1591: lv[k][j] = 0.;
1592: rv[k][j] = 0.;
1593: }
1594: rho = ueq[0];
1595: /* u = ueq[1]; */
1596: p0 = ueq[2];
1597: csnd = PetscSqrtReal(gamma * p0 / rho);
1598: lv[0][1] = rho * .5;
1599: lv[0][2] = -.5 / csnd;
1600: lv[1][0] = csnd;
1601: lv[1][2] = -1. / csnd;
1602: lv[2][1] = rho * .5;
1603: lv[2][2] = .5 / csnd;
1604: rv[0][0] = -1. / csnd;
1605: rv[1][0] = 1. / rho;
1606: rv[2][0] = -csnd;
1607: rv[0][1] = 1. / csnd;
1608: rv[0][2] = 1. / csnd;
1609: rv[1][2] = 1. / rho;
1610: rv[2][2] = csnd;
1611: return 0;
1612: }
1614: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1615: {
1616: PetscReal p0, u0, wcp[3], wc[3];
1617: PetscReal lv[3][3];
1618: PetscReal vp[3];
1619: PetscReal rv[3][3];
1620: PetscReal eps, ueq[3], rho0, twopi;
1622: /* Function Body */
1623: twopi = 2. * PETSC_PI;
1624: eps = 1e-4; /* perturbation */
1625: rho0 = 1e3; /* density of water */
1626: p0 = 101325.; /* init pressure of 1 atm (?) */
1627: u0 = 0.;
1628: ueq[0] = rho0;
1629: ueq[1] = u0;
1630: ueq[2] = p0;
1631: /* Project initial state to characteristic variables */
1632: eigenvectors(rv, lv, ueq, gamma);
1633: projecteqstate(wc, ueq, lv);
1634: wcp[0] = wc[0];
1635: wcp[1] = wc[1];
1636: wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1637: projecttoprim(vp, wcp, rv);
1638: ux->r = vp[0]; /* density */
1639: ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1640: ux->ru[1] = 0.;
1641: /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1642: ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1643: return 0;
1644: }
1646: /*TEST
1648: testset:
1649: args: -dm_plex_adj_cone -dm_plex_adj_closure 0
1651: test:
1652: suffix: adv_2d_tri_0
1653: requires: triangle
1654: TODO: how did this ever get in main when there is no support for this
1655: 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
1657: test:
1658: suffix: adv_2d_tri_1
1659: requires: triangle
1660: TODO: how did this ever get in main when there is no support for this
1661: 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
1663: test:
1664: suffix: tut_1
1665: requires: exodusii
1666: nsize: 1
1667: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -
1669: test:
1670: suffix: tut_2
1671: requires: exodusii
1672: nsize: 1
1673: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw
1675: test:
1676: suffix: tut_3
1677: requires: exodusii
1678: nsize: 4
1679: 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
1681: test:
1682: suffix: tut_4
1683: requires: exodusii !single
1684: nsize: 4
1685: 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
1687: testset:
1688: args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
1690: # 2D Advection 0-10
1691: test:
1692: suffix: 0
1693: requires: exodusii
1694: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1696: test:
1697: suffix: 1
1698: requires: exodusii
1699: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1701: test:
1702: suffix: 2
1703: requires: exodusii
1704: nsize: 2
1705: args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1707: test:
1708: suffix: 3
1709: requires: exodusii
1710: nsize: 2
1711: args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1713: test:
1714: suffix: 4
1715: requires: exodusii
1716: nsize: 4
1717: args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple
1719: test:
1720: suffix: 5
1721: requires: exodusii
1722: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1
1724: test:
1725: suffix: 7
1726: requires: exodusii
1727: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1729: test:
1730: suffix: 8
1731: requires: exodusii
1732: nsize: 2
1733: 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
1735: test:
1736: suffix: 9
1737: requires: exodusii
1738: nsize: 8
1739: 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
1741: test:
1742: suffix: 10
1743: requires: exodusii
1744: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo
1746: # 2D Shallow water
1747: testset:
1748: args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0
1750: test:
1751: suffix: sw_0
1752: requires: exodusii
1753: args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1754: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1755: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_num_stages 10 \
1756: -monitor height,energy
1758: test:
1759: suffix: sw_ceed
1760: requires: exodusii libceed
1761: args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1762: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1763: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_num_stages 10 \
1764: -monitor height,energy
1766: test:
1767: suffix: sw_ceed_small
1768: requires: exodusii libceed
1769: args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \
1770: -dm_plex_shape annulus -dm_plex_simplex 0 -dm_plex_box_lower 0,1 -dm_plex_box_upper 6.28,3 -dm_plex_box_faces 8,2 \
1771: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_num_stages 10 \
1772: -monitor height,energy
1774: test:
1775: suffix: sw_1
1776: nsize: 2
1777: args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1778: -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 \
1779: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_num_stages 10 \
1780: -monitor height,energy
1782: test:
1783: suffix: sw_hll
1784: args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1785: -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1786: -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_num_stages 10 \
1787: -monitor height,energy
1789: # 2D Euler
1790: testset:
1791: args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1792: -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy
1794: test:
1795: suffix: euler_0
1796: requires: exodusii !complex !single
1797: args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1798: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1799: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_num_stages 10
1801: test:
1802: suffix: euler_ceed
1803: requires: exodusii libceed
1804: args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1805: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1806: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_num_stages 10
1808: testset:
1809: args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
1811: # 2D Advection: p4est
1812: test:
1813: suffix: p4est_advec_2d
1814: requires: p4est
1815: 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
1817: # Advection in a box
1818: test:
1819: suffix: adv_2d_quad_0
1820: args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1822: test:
1823: suffix: adv_2d_quad_1
1824: 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
1825: timeoutfactor: 3
1827: test:
1828: suffix: adv_2d_quad_p4est_0
1829: requires: p4est
1830: args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1832: test:
1833: suffix: adv_2d_quad_p4est_1
1834: requires: p4est
1835: 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
1836: timeoutfactor: 3
1838: test: # broken for quad precision
1839: suffix: adv_2d_quad_p4est_adapt_0
1840: requires: p4est !__float128
1841: 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
1842: timeoutfactor: 3
1844: test:
1845: suffix: adv_0
1846: requires: exodusii
1847: 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
1849: # Run with -dm_forest_maximum_refinement 6 -ts_max_time 0.5 instead to get the full movie
1850: test:
1851: suffix: shock_0
1852: requires: p4est !single !complex
1853: args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
1854: -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 2 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
1855: -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
1856: -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. \
1857: -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
1858: -ts_max_steps 3 -ts_ssp_type rks2 -ts_ssp_num_stages 10 \
1859: -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
1861: # Test GLVis visualization of PetscFV fields
1862: test:
1863: suffix: glvis_adv_2d_tri
1864: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1865: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1866: -ts_monitor_solution glvis: -ts_max_steps 0
1868: test:
1869: suffix: glvis_adv_2d_quad
1870: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
1871: -dm_refine 5 -dm_plex_separate_marker \
1872: -ts_monitor_solution glvis: -ts_max_steps 0
1874: # Test CGNS file writing for PetscFV fields
1875: test:
1876: suffix: cgns_adv_2d_tri
1877: requires: cgns
1878: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1879: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1880: -ts_monitor_solution cgns:sol.cgns -ts_max_steps 0
1882: test:
1883: suffix: cgns_adv_2d_quad
1884: requires: cgns
1885: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
1886: -dm_refine 5 -dm_plex_separate_marker \
1887: -ts_monitor_solution cgns:sol.cgns -ts_max_steps 0
1889: # Test CGNS file writing, cgns_batch_size, ts_run_steps, ts_monitor_skip_initial
1890: test:
1891: suffix: cgns_adv_2d_tri_monitor
1892: requires: cgns
1893: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1894: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1895: -ts_monitor_solution cgns:sol-%d.cgns -ts_run_steps 4 -ts_monitor_solution_interval 2 \
1896: -viewer_cgns_batch_size 1 -ts_monitor_solution_skip_initial
1898: # Test VTK file writing for PetscFV fields with -ts_monitor_solution_vtk
1899: test:
1900: suffix: vtk_adv_2d_tri
1901: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1902: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1903: -ts_monitor_solution_vtk 'bar-%03d.vtu' -ts_monitor_solution_vtk_interval 2 -ts_max_steps 5
1905: TEST*/