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