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: PetscBool useCeed;
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(DMPlexGetUseCeed(dm, &useCeed));
1079: PetscCall(DMTSSetBoundaryLocal(dm, DMPlexTSComputeBoundary, user));
1080: if (useCeed) PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVMCEED, user));
1081: else PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
1082: PetscCall(TSSetMaxTime(*ts, 2.0));
1083: PetscCall(TSSetExactFinalTime(*ts, TS_EXACTFINALTIME_STEPOVER));
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: typedef struct {
1088: PetscFV fvm;
1089: VecTagger refineTag;
1090: VecTagger coarsenTag;
1091: DM adaptedDM;
1092: User user;
1093: PetscReal cfl;
1094: PetscLimiter limiter;
1095: PetscLimiter noneLimiter;
1096: } TransferCtx;
1098: static PetscErrorCode adaptToleranceFVMSetUp(TS ts, PetscInt nstep, PetscReal time, Vec sol, PetscBool *resize, void *ctx)
1099: {
1100: TransferCtx *tctx = (TransferCtx *)ctx;
1101: PetscFV fvm = tctx->fvm;
1102: VecTagger refineTag = tctx->refineTag;
1103: VecTagger coarsenTag = tctx->coarsenTag;
1104: User user = tctx->user;
1105: DM dm, gradDM, plex, cellDM, adaptedDM = NULL;
1106: Vec cellGeom, faceGeom;
1107: PetscBool computeGradient;
1108: Vec grad, locGrad, locX, errVec;
1109: PetscInt cStart, cEnd, c, dim, nRefine, nCoarsen;
1110: PetscReal minMaxInd[2] = {PETSC_MAX_REAL, PETSC_MIN_REAL}, minMaxIndGlobal[2];
1111: PetscScalar *errArray;
1112: const PetscScalar *pointVals;
1113: const PetscScalar *pointGrads;
1114: const PetscScalar *pointGeom;
1115: DMLabel adaptLabel = NULL;
1116: IS refineIS, coarsenIS;
1118: PetscFunctionBeginUser;
1119: *resize = PETSC_FALSE;
1120: PetscCall(VecGetDM(sol, &dm));
1121: PetscCall(DMGetDimension(dm, &dim));
1122: PetscCall(PetscFVSetLimiter(fvm, tctx->noneLimiter));
1123: PetscCall(PetscFVGetComputeGradients(fvm, &computeGradient));
1124: PetscCall(PetscFVSetComputeGradients(fvm, PETSC_TRUE));
1125: PetscCall(DMConvert(dm, DMPLEX, &plex));
1126: PetscCall(DMPlexGetDataFVM(plex, fvm, &cellGeom, &faceGeom, &gradDM));
1127: PetscCall(DMCreateLocalVector(plex, &locX));
1128: PetscCall(DMPlexInsertBoundaryValues(plex, PETSC_TRUE, locX, 0.0, faceGeom, cellGeom, NULL));
1129: PetscCall(DMGlobalToLocalBegin(plex, sol, INSERT_VALUES, locX));
1130: PetscCall(DMGlobalToLocalEnd(plex, sol, INSERT_VALUES, locX));
1131: PetscCall(DMCreateGlobalVector(gradDM, &grad));
1132: PetscCall(DMPlexReconstructGradientsFVM(plex, locX, grad));
1133: PetscCall(DMCreateLocalVector(gradDM, &locGrad));
1134: PetscCall(DMGlobalToLocalBegin(gradDM, grad, INSERT_VALUES, locGrad));
1135: PetscCall(DMGlobalToLocalEnd(gradDM, grad, INSERT_VALUES, locGrad));
1136: PetscCall(VecDestroy(&grad));
1137: PetscCall(DMPlexGetSimplexOrBoxCells(plex, 0, &cStart, &cEnd));
1138: PetscCall(VecGetArrayRead(locGrad, &pointGrads));
1139: PetscCall(VecGetArrayRead(cellGeom, &pointGeom));
1140: PetscCall(VecGetArrayRead(locX, &pointVals));
1141: PetscCall(VecGetDM(cellGeom, &cellDM));
1142: PetscCall(DMLabelCreate(PETSC_COMM_SELF, "adapt", &adaptLabel));
1143: PetscCall(VecCreateFromOptions(PetscObjectComm((PetscObject)plex), NULL, 1, cEnd - cStart, PETSC_DETERMINE, &errVec));
1144: PetscCall(VecSetUp(errVec));
1145: PetscCall(VecGetArray(errVec, &errArray));
1146: for (c = cStart; c < cEnd; c++) {
1147: PetscReal errInd = 0.;
1148: PetscScalar *pointGrad;
1149: PetscScalar *pointVal;
1150: PetscFVCellGeom *cg;
1152: PetscCall(DMPlexPointLocalRead(gradDM, c, pointGrads, &pointGrad));
1153: PetscCall(DMPlexPointLocalRead(cellDM, c, pointGeom, &cg));
1154: PetscCall(DMPlexPointLocalRead(plex, c, pointVals, &pointVal));
1156: PetscCall((*user->model->errorIndicator)(dim, cg->volume, user->model->physics->dof, pointVal, pointGrad, &errInd, user->model->errorCtx));
1157: errArray[c - cStart] = errInd;
1158: minMaxInd[0] = PetscMin(minMaxInd[0], errInd);
1159: minMaxInd[1] = PetscMax(minMaxInd[1], errInd);
1160: }
1161: PetscCall(VecRestoreArray(errVec, &errArray));
1162: PetscCall(VecRestoreArrayRead(locX, &pointVals));
1163: PetscCall(VecRestoreArrayRead(cellGeom, &pointGeom));
1164: PetscCall(VecRestoreArrayRead(locGrad, &pointGrads));
1165: PetscCall(VecDestroy(&locGrad));
1166: PetscCall(VecDestroy(&locX));
1167: PetscCall(DMDestroy(&plex));
1169: PetscCall(VecTaggerComputeIS(refineTag, errVec, &refineIS, NULL));
1170: PetscCall(VecTaggerComputeIS(coarsenTag, errVec, &coarsenIS, NULL));
1171: PetscCall(ISGetSize(refineIS, &nRefine));
1172: PetscCall(ISGetSize(coarsenIS, &nCoarsen));
1173: if (nRefine) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_REFINE, refineIS));
1174: if (nCoarsen) PetscCall(DMLabelSetStratumIS(adaptLabel, DM_ADAPT_COARSEN, coarsenIS));
1175: PetscCall(ISDestroy(&coarsenIS));
1176: PetscCall(ISDestroy(&refineIS));
1177: PetscCall(VecDestroy(&errVec));
1179: PetscCall(PetscFVSetComputeGradients(fvm, computeGradient));
1180: PetscCall(PetscFVSetLimiter(fvm, tctx->limiter));
1181: minMaxInd[1] = -minMaxInd[1];
1182: PetscCallMPI(MPIU_Allreduce(minMaxInd, minMaxIndGlobal, 2, MPIU_REAL, MPI_MIN, PetscObjectComm((PetscObject)dm)));
1183: PetscCall(PetscInfo(ts, "error indicator range (%E, %E)\n", (double)minMaxIndGlobal[0], (double)(-minMaxIndGlobal[1])));
1184: if (nRefine || nCoarsen) { /* at least one cell is over the refinement threshold */
1185: PetscCall(DMAdaptLabel(dm, adaptLabel, &adaptedDM));
1186: }
1187: PetscCall(DMLabelDestroy(&adaptLabel));
1188: if (adaptedDM) {
1189: tctx->adaptedDM = adaptedDM;
1190: PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " adapted mesh, marking %" PetscInt_FMT " cells for refinement, and %" PetscInt_FMT " cells for coarsening\n", nstep, nRefine, nCoarsen));
1191: *resize = PETSC_TRUE;
1192: } else {
1193: PetscCall(PetscInfo(ts, "Step %" PetscInt_FMT " no adaptation\n", nstep));
1194: }
1195: PetscFunctionReturn(PETSC_SUCCESS);
1196: }
1198: static PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vecsin[], Vec vecsout[], void *ctx)
1199: {
1200: TransferCtx *tctx = (TransferCtx *)ctx;
1201: DM dm;
1202: PetscReal time;
1204: PetscFunctionBeginUser;
1205: PetscCall(TSGetDM(ts, &dm));
1206: PetscCall(TSGetTime(ts, &time));
1207: PetscCheck(tctx->adaptedDM, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Missing adaptedDM");
1208: for (PetscInt i = 0; i < nv; i++) {
1209: PetscCall(DMCreateGlobalVector(tctx->adaptedDM, &vecsout[i]));
1210: PetscCall(DMForestTransferVec(dm, vecsin[i], tctx->adaptedDM, vecsout[i], PETSC_TRUE, time));
1211: }
1212: PetscCall(DMForestSetAdaptivityForest(tctx->adaptedDM, NULL)); /* clear internal references to the previous dm */
1214: Model mod = tctx->user->model;
1215: Physics phys = mod->physics;
1216: PetscReal minRadius;
1218: PetscCall(DMPlexGetGeometryFVM(tctx->adaptedDM, NULL, NULL, &minRadius));
1219: PetscCallMPI(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1220: PetscCheck(mod->maxspeed > 0, PetscObjectComm((PetscObject)ts), PETSC_ERR_ARG_WRONGSTATE, "Physics did not set maxspeed");
1222: PetscReal dt = tctx->cfl * minRadius / mod->maxspeed;
1223: PetscCall(TSSetTimeStep(ts, dt));
1225: PetscCall(TSSetDM(ts, tctx->adaptedDM));
1226: PetscCall(DMDestroy(&tctx->adaptedDM));
1227: PetscFunctionReturn(PETSC_SUCCESS);
1228: }
1230: int main(int argc, char **argv)
1231: {
1232: MPI_Comm comm;
1233: PetscDS prob;
1234: PetscFV fvm;
1235: PetscLimiter limiter = NULL, noneLimiter = NULL;
1236: User user;
1237: Model mod;
1238: Physics phys;
1239: DM dm, plex;
1240: PetscReal ftime, cfl, dt, minRadius;
1241: PetscInt dim, nsteps;
1242: TS ts;
1243: TSConvergedReason reason;
1244: Vec X;
1245: PetscViewer viewer;
1246: PetscBool vtkCellGeom, useAMR;
1247: PetscInt adaptInterval;
1248: char physname[256] = "advect";
1249: VecTagger refineTag = NULL, coarsenTag = NULL;
1250: TransferCtx tctx;
1252: PetscFunctionBeginUser;
1253: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1254: comm = PETSC_COMM_WORLD;
1256: PetscCall(PetscNew(&user));
1257: PetscCall(PetscNew(&user->model));
1258: PetscCall(PetscNew(&user->model->physics));
1259: mod = user->model;
1260: phys = mod->physics;
1261: mod->comm = comm;
1262: useAMR = PETSC_FALSE;
1263: adaptInterval = 1;
1265: /* Register physical models to be available on the command line */
1266: PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
1267: PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
1268: PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1270: PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1271: {
1272: cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1273: PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1274: user->vtkInterval = 1;
1275: PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1276: user->vtkmon = PETSC_TRUE;
1277: PetscCall(PetscOptionsBool("-ufv_vtk_monitor", "Use VTKMonitor routine", "", user->vtkmon, &user->vtkmon, NULL));
1278: vtkCellGeom = PETSC_FALSE;
1279: PetscCall(PetscStrncpy(user->outputBasename, "ex11", sizeof(user->outputBasename)));
1280: PetscCall(PetscOptionsString("-ufv_vtk_basename", "VTK output basename", "", user->outputBasename, user->outputBasename, sizeof(user->outputBasename), NULL));
1281: PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
1282: PetscCall(PetscOptionsBool("-ufv_use_amr", "use local adaptive mesh refinement", "", useAMR, &useAMR, NULL));
1283: PetscCall(PetscOptionsInt("-ufv_adapt_interval", "time steps between AMR", "", adaptInterval, &adaptInterval, NULL));
1284: }
1285: PetscOptionsEnd();
1287: if (useAMR) {
1288: VecTaggerBox refineBox, coarsenBox;
1290: refineBox.min = refineBox.max = PETSC_MAX_REAL;
1291: coarsenBox.min = coarsenBox.max = PETSC_MIN_REAL;
1293: PetscCall(VecTaggerCreate(comm, &refineTag));
1294: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)refineTag, "refine_"));
1295: PetscCall(VecTaggerSetType(refineTag, VECTAGGERABSOLUTE));
1296: PetscCall(VecTaggerAbsoluteSetBox(refineTag, &refineBox));
1297: PetscCall(VecTaggerSetFromOptions(refineTag));
1298: PetscCall(VecTaggerSetUp(refineTag));
1299: PetscCall(PetscObjectViewFromOptions((PetscObject)refineTag, NULL, "-tag_view"));
1301: PetscCall(VecTaggerCreate(comm, &coarsenTag));
1302: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)coarsenTag, "coarsen_"));
1303: PetscCall(VecTaggerSetType(coarsenTag, VECTAGGERABSOLUTE));
1304: PetscCall(VecTaggerAbsoluteSetBox(coarsenTag, &coarsenBox));
1305: PetscCall(VecTaggerSetFromOptions(coarsenTag));
1306: PetscCall(VecTaggerSetUp(coarsenTag));
1307: PetscCall(PetscObjectViewFromOptions((PetscObject)coarsenTag, NULL, "-tag_view"));
1308: }
1310: PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1311: {
1312: PetscErrorCode (*physcreate)(Model, Physics, PetscOptionItems *);
1313: PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
1314: PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
1315: PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
1316: PetscCall((*physcreate)(mod, phys, PetscOptionsObject));
1317: /* Count number of fields and dofs */
1318: for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1319: PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1320: PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1321: }
1322: PetscOptionsEnd();
1324: /* Create mesh */
1325: {
1326: PetscInt i;
1328: PetscCall(DMCreate(comm, &dm));
1329: PetscCall(DMSetType(dm, DMPLEX));
1330: PetscCall(DMSetFromOptions(dm));
1331: for (i = 0; i < DIM; i++) {
1332: mod->bounds[2 * i] = 0.;
1333: mod->bounds[2 * i + 1] = 1.;
1334: }
1335: dim = DIM;
1336: { /* a null name means just do a hex box */
1337: PetscInt cells[3] = {1, 1, 1}, n = 3;
1338: PetscBool flg2, skew = PETSC_FALSE;
1339: PetscInt nret2 = 2 * DIM;
1340: PetscOptionsBegin(comm, NULL, "Rectangular mesh options", "");
1341: 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));
1342: PetscCall(PetscOptionsBool("-grid_skew_60", "Skew grid for 60 degree shock mesh", "", skew, &skew, NULL));
1343: PetscCall(PetscOptionsIntArray("-dm_plex_box_faces", "Number of faces along each dimension", "", cells, &n, NULL));
1344: PetscOptionsEnd();
1345: /* TODO Rewrite this with Mark, and remove grid_bounds at that time */
1346: if (flg2) {
1347: PetscInt dimEmbed, i;
1348: PetscInt nCoords;
1349: PetscScalar *coords;
1350: Vec coordinates;
1352: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1353: PetscCall(DMGetCoordinateDim(dm, &dimEmbed));
1354: PetscCall(VecGetLocalSize(coordinates, &nCoords));
1355: PetscCheck(!(nCoords % dimEmbed), PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Coordinate vector the wrong size");
1356: PetscCall(VecGetArray(coordinates, &coords));
1357: for (i = 0; i < nCoords; i += dimEmbed) {
1358: PetscInt j;
1360: PetscScalar *coord = &coords[i];
1361: for (j = 0; j < dimEmbed; j++) {
1362: coord[j] = mod->bounds[2 * j] + coord[j] * (mod->bounds[2 * j + 1] - mod->bounds[2 * j]);
1363: if (dim == 2 && cells[1] == 1 && j == 0 && skew) {
1364: if (cells[0] == 2 && i == 8) {
1365: coord[j] = .57735026918963; /* hack to get 60 deg skewed mesh */
1366: } else if (cells[0] == 3) {
1367: if (i == 2 || i == 10) coord[j] = mod->bounds[1] / 4.;
1368: else if (i == 4) coord[j] = mod->bounds[1] / 2.;
1369: else if (i == 12) coord[j] = 1.57735026918963 * mod->bounds[1] / 2.;
1370: }
1371: }
1372: }
1373: }
1374: PetscCall(VecRestoreArray(coordinates, &coords));
1375: PetscCall(DMSetCoordinatesLocal(dm, coordinates));
1376: }
1377: }
1378: }
1379: PetscCall(DMViewFromOptions(dm, NULL, "-orig_dm_view"));
1380: PetscCall(DMGetDimension(dm, &dim));
1382: /* set up BCs, functions, tags */
1383: PetscCall(DMCreateLabel(dm, "Face Sets"));
1384: mod->errorIndicator = ErrorIndicator_Simple;
1386: {
1387: DM gdm;
1389: PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
1390: PetscCall(DMDestroy(&dm));
1391: dm = gdm;
1392: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1393: }
1395: PetscCall(PetscFVCreate(comm, &fvm));
1396: PetscCall(PetscFVSetFromOptions(fvm));
1397: PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
1398: PetscCall(PetscFVSetSpatialDimension(fvm, dim));
1399: PetscCall(PetscObjectSetName((PetscObject)fvm, ""));
1400: {
1401: PetscInt f, dof;
1402: for (f = 0, dof = 0; f < phys->nfields; f++) {
1403: PetscInt newDof = phys->field_desc[f].dof;
1405: if (newDof == 1) {
1406: PetscCall(PetscFVSetComponentName(fvm, dof, phys->field_desc[f].name));
1407: } else {
1408: PetscInt j;
1410: for (j = 0; j < newDof; j++) {
1411: char compName[256] = "Unknown";
1413: PetscCall(PetscSNPrintf(compName, sizeof(compName), "%s_%" PetscInt_FMT, phys->field_desc[f].name, j));
1414: PetscCall(PetscFVSetComponentName(fvm, dof + j, compName));
1415: }
1416: }
1417: dof += newDof;
1418: }
1419: }
1420: /* FV is now structured with one field having all physics as components */
1421: PetscCall(DMAddField(dm, NULL, (PetscObject)fvm));
1422: PetscCall(DMCreateDS(dm));
1423: PetscCall(DMGetDS(dm, &prob));
1424: PetscCall(PetscDSSetRiemannSolver(prob, 0, user->model->physics->riemann));
1425: PetscCall(PetscDSSetContext(prob, 0, user->model->physics));
1426: PetscCall((*mod->setupbc)(dm, prob, phys));
1427: PetscCall(PetscDSSetFromOptions(prob));
1428: {
1429: char convType[256];
1430: PetscBool flg;
1432: PetscOptionsBegin(comm, "", "Mesh conversion options", "DMPLEX");
1433: PetscCall(PetscOptionsFList("-dm_type", "Convert DMPlex to another format", "ex12", DMList, DMPLEX, convType, 256, &flg));
1434: PetscOptionsEnd();
1435: if (flg) {
1436: DM dmConv;
1438: PetscCall(DMConvert(dm, convType, &dmConv));
1439: if (dmConv) {
1440: PetscCall(DMViewFromOptions(dmConv, NULL, "-dm_conv_view"));
1441: PetscCall(DMDestroy(&dm));
1442: dm = dmConv;
1443: PetscCall(DMSetFromOptions(dm));
1444: }
1445: }
1446: }
1447: #ifdef PETSC_HAVE_LIBCEED
1448: {
1449: PetscBool useCeed;
1450: PetscCall(DMPlexGetUseCeed(dm, &useCeed));
1451: if (useCeed) PetscCall((*user->model->setupCEED)(dm, user->model->physics));
1452: }
1453: #endif
1455: PetscCall(initializeTS(dm, user, &ts));
1457: PetscCall(DMCreateGlobalVector(dm, &X));
1458: PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1459: PetscCall(SetInitialCondition(dm, X, user));
1460: if (useAMR) {
1461: PetscInt adaptIter;
1463: /* use no limiting when reconstructing gradients for adaptivity */
1464: PetscCall(PetscFVGetLimiter(fvm, &limiter));
1465: PetscCall(PetscObjectReference((PetscObject)limiter));
1466: PetscCall(PetscLimiterCreate(PetscObjectComm((PetscObject)fvm), &noneLimiter));
1467: PetscCall(PetscLimiterSetType(noneLimiter, PETSCLIMITERNONE));
1469: /* Refinement context */
1470: tctx.fvm = fvm;
1471: tctx.refineTag = refineTag;
1472: tctx.coarsenTag = coarsenTag;
1473: tctx.adaptedDM = NULL;
1474: tctx.user = user;
1475: tctx.noneLimiter = noneLimiter;
1476: tctx.limiter = limiter;
1477: tctx.cfl = cfl;
1479: /* Do some initial refinement steps */
1480: for (adaptIter = 0;; ++adaptIter) {
1481: PetscLogDouble bytes;
1482: PetscBool resize;
1484: PetscCall(PetscMemoryGetCurrentUsage(&bytes));
1485: PetscCall(PetscInfo(ts, "refinement loop %" PetscInt_FMT ": memory used %g\n", adaptIter, (double)bytes));
1486: PetscCall(DMViewFromOptions(dm, NULL, "-initial_dm_view"));
1487: PetscCall(VecViewFromOptions(X, NULL, "-initial_vec_view"));
1489: PetscCall(adaptToleranceFVMSetUp(ts, -1, 0.0, X, &resize, &tctx));
1490: if (!resize) break;
1491: PetscCall(DMDestroy(&dm));
1492: PetscCall(VecDestroy(&X));
1493: dm = tctx.adaptedDM;
1494: tctx.adaptedDM = NULL;
1495: PetscCall(TSSetDM(ts, dm));
1496: PetscCall(DMCreateGlobalVector(dm, &X));
1497: PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1498: PetscCall(SetInitialCondition(dm, X, user));
1499: }
1500: }
1502: PetscCall(DMConvert(dm, DMPLEX, &plex));
1503: if (vtkCellGeom) {
1504: DM dmCell;
1505: Vec cellgeom, partition;
1507: PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1508: PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
1509: PetscCall(VecView(cellgeom, viewer));
1510: PetscCall(PetscViewerDestroy(&viewer));
1511: PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
1512: PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
1513: PetscCall(VecView(partition, viewer));
1514: PetscCall(PetscViewerDestroy(&viewer));
1515: PetscCall(VecDestroy(&partition));
1516: PetscCall(DMDestroy(&dmCell));
1517: }
1518: /* collect max maxspeed from all processes -- todo */
1519: PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
1520: PetscCall(DMDestroy(&plex));
1521: PetscCallMPI(MPIU_Allreduce(&phys->maxspeed, &mod->maxspeed, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1522: PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1523: dt = cfl * minRadius / mod->maxspeed;
1524: PetscCall(TSSetTimeStep(ts, dt));
1525: PetscCall(TSSetFromOptions(ts));
1527: /* When using adaptive mesh refinement
1528: specify callbacks to refine the solution
1529: and interpolate data from old to new mesh
1530: When mesh adaption is requested, the step will be rolled back
1531: */
1532: if (useAMR) PetscCall(TSSetResize(ts, PETSC_TRUE, adaptToleranceFVMSetUp, Transfer, &tctx));
1533: PetscCall(TSSetSolution(ts, X));
1534: PetscCall(VecDestroy(&X));
1535: PetscCall(TSSolve(ts, NULL));
1536: PetscCall(TSGetSolveTime(ts, &ftime));
1537: PetscCall(TSGetStepNumber(ts, &nsteps));
1538: PetscCall(TSGetConvergedReason(ts, &reason));
1539: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
1540: PetscCall(TSDestroy(&ts));
1542: PetscCall(VecTaggerDestroy(&refineTag));
1543: PetscCall(VecTaggerDestroy(&coarsenTag));
1544: PetscCall(PetscFunctionListDestroy(&PhysicsList));
1545: PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_SW));
1546: PetscCall(PetscFunctionListDestroy(&PhysicsRiemannList_Euler));
1547: PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
1548: PetscCall(PetscFree(user->model->functionalMonitored));
1549: PetscCall(PetscFree(user->model->functionalCall));
1550: PetscCall(PetscFree(user->model->physics->data));
1551: PetscCall(PetscFree(user->model->physics));
1552: PetscCall(PetscFree(user->model));
1553: PetscCall(PetscFree(user));
1554: PetscCall(PetscLimiterDestroy(&limiter));
1555: PetscCall(PetscLimiterDestroy(&noneLimiter));
1556: PetscCall(PetscFVDestroy(&fvm));
1557: PetscCall(DMDestroy(&dm));
1558: PetscCall(PetscFinalize());
1559: return 0;
1560: }
1562: /* Subroutine to set up the initial conditions for the */
1563: /* Shock Interface interaction or linear wave (Ravi Samtaney,Mark Adams). */
1564: /* ----------------------------------------------------------------------- */
1565: int projecteqstate(PetscReal wc[], const PetscReal ueq[], PetscReal lv[][3])
1566: {
1567: int j, k;
1568: /* Wc=matmul(lv,Ueq) 3 vars */
1569: for (k = 0; k < 3; ++k) {
1570: wc[k] = 0.;
1571: for (j = 0; j < 3; ++j) wc[k] += lv[k][j] * ueq[j];
1572: }
1573: return 0;
1574: }
1575: /* ----------------------------------------------------------------------- */
1576: int projecttoprim(PetscReal v[], const PetscReal wc[], PetscReal rv[][3])
1577: {
1578: int k, j;
1579: /* V=matmul(rv,WC) 3 vars */
1580: for (k = 0; k < 3; ++k) {
1581: v[k] = 0.;
1582: for (j = 0; j < 3; ++j) v[k] += rv[k][j] * wc[j];
1583: }
1584: return 0;
1585: }
1586: /* ---------------------------------------------------------------------- */
1587: int eigenvectors(PetscReal rv[][3], PetscReal lv[][3], const PetscReal ueq[], PetscReal gamma)
1588: {
1589: int j, k;
1590: PetscReal rho, csnd, p0;
1591: /* PetscScalar u; */
1593: for (k = 0; k < 3; ++k)
1594: for (j = 0; j < 3; ++j) {
1595: lv[k][j] = 0.;
1596: rv[k][j] = 0.;
1597: }
1598: rho = ueq[0];
1599: /* u = ueq[1]; */
1600: p0 = ueq[2];
1601: csnd = PetscSqrtReal(gamma * p0 / rho);
1602: lv[0][1] = rho * .5;
1603: lv[0][2] = -.5 / csnd;
1604: lv[1][0] = csnd;
1605: lv[1][2] = -1. / csnd;
1606: lv[2][1] = rho * .5;
1607: lv[2][2] = .5 / csnd;
1608: rv[0][0] = -1. / csnd;
1609: rv[1][0] = 1. / rho;
1610: rv[2][0] = -csnd;
1611: rv[0][1] = 1. / csnd;
1612: rv[0][2] = 1. / csnd;
1613: rv[1][2] = 1. / rho;
1614: rv[2][2] = csnd;
1615: return 0;
1616: }
1618: int initLinearWave(EulerNode *ux, const PetscReal gamma, const PetscReal coord[], const PetscReal Lx)
1619: {
1620: PetscReal p0, u0, wcp[3], wc[3];
1621: PetscReal lv[3][3];
1622: PetscReal vp[3];
1623: PetscReal rv[3][3];
1624: PetscReal eps, ueq[3], rho0, twopi;
1626: /* Function Body */
1627: twopi = 2. * PETSC_PI;
1628: eps = 1e-4; /* perturbation */
1629: rho0 = 1e3; /* density of water */
1630: p0 = 101325.; /* init pressure of 1 atm (?) */
1631: u0 = 0.;
1632: ueq[0] = rho0;
1633: ueq[1] = u0;
1634: ueq[2] = p0;
1635: /* Project initial state to characteristic variables */
1636: eigenvectors(rv, lv, ueq, gamma);
1637: projecteqstate(wc, ueq, lv);
1638: wcp[0] = wc[0];
1639: wcp[1] = wc[1];
1640: wcp[2] = wc[2] + eps * PetscCosReal(coord[0] * 2. * twopi / Lx);
1641: projecttoprim(vp, wcp, rv);
1642: ux->r = vp[0]; /* density */
1643: ux->ru[0] = vp[0] * vp[1]; /* x momentum */
1644: ux->ru[1] = 0.;
1645: #if defined DIM > 2
1646: if (dim > 2) ux->ru[2] = 0.;
1647: #endif
1648: /* E = rho * e + rho * v^2/2 = p/(gam-1) + rho*v^2/2 */
1649: ux->E = vp[2] / (gamma - 1.) + 0.5 * vp[0] * vp[1] * vp[1];
1650: return 0;
1651: }
1653: /*TEST
1655: testset:
1656: args: -dm_plex_adj_cone -dm_plex_adj_closure 0
1658: test:
1659: suffix: adv_2d_tri_0
1660: requires: triangle
1661: TODO: how did this ever get in main when there is no support for this
1662: 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
1664: test:
1665: suffix: adv_2d_tri_1
1666: requires: triangle
1667: TODO: how did this ever get in main when there is no support for this
1668: 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
1670: test:
1671: suffix: tut_1
1672: requires: exodusii
1673: nsize: 1
1674: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -
1676: test:
1677: suffix: tut_2
1678: requires: exodusii
1679: nsize: 1
1680: args: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw
1682: test:
1683: suffix: tut_3
1684: requires: exodusii
1685: nsize: 4
1686: 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
1688: test:
1689: suffix: tut_4
1690: requires: exodusii !single
1691: nsize: 4
1692: 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
1694: testset:
1695: args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
1697: # 2D Advection 0-10
1698: test:
1699: suffix: 0
1700: requires: exodusii
1701: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1703: test:
1704: suffix: 1
1705: requires: exodusii
1706: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1708: test:
1709: suffix: 2
1710: requires: exodusii
1711: nsize: 2
1712: args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo
1714: test:
1715: suffix: 3
1716: requires: exodusii
1717: nsize: 2
1718: args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo
1720: test:
1721: suffix: 4
1722: requires: exodusii
1723: nsize: 4
1724: args: -dm_distribute_overlap 1 -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo -petscpartitioner_type simple
1726: test:
1727: suffix: 5
1728: requires: exodusii
1729: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside.exo -ts_type rosw -ts_adapt_reject_safety 1
1731: test:
1732: suffix: 7
1733: requires: exodusii
1734: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad-15.exo -dm_refine 1
1736: test:
1737: suffix: 8
1738: requires: exodusii
1739: nsize: 2
1740: 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
1742: test:
1743: suffix: 9
1744: requires: exodusii
1745: nsize: 8
1746: 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
1748: test:
1749: suffix: 10
1750: requires: exodusii
1751: args: -ufv_vtk_interval 0 -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/sevenside-quad.exo
1753: # 2D Shallow water
1754: testset:
1755: args: -physics sw -ufv_vtk_interval 0 -dm_plex_adj_cone -dm_plex_adj_closure 0
1757: test:
1758: suffix: sw_0
1759: requires: exodusii
1760: args: -bc_wall 100,101 -ufv_cfl 5 -petscfv_type leastsquares -petsclimiter_type sin \
1761: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1762: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1763: -monitor height,energy
1765: test:
1766: suffix: sw_ceed
1767: requires: exodusii libceed
1768: args: -sw_riemann rusanov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1769: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1770: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1771: -monitor height,energy
1773: test:
1774: suffix: sw_ceed_small
1775: requires: exodusii libceed
1776: args: -sw_riemann rusanov_ceed -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin -dm_plex_use_ceed \
1777: -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 \
1778: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1779: -monitor height,energy
1781: test:
1782: suffix: sw_1
1783: nsize: 2
1784: args: -bc_wall 1,3 -ufv_cfl 5 -petsclimiter_type sin \
1785: -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 \
1786: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1787: -monitor height,energy
1789: test:
1790: suffix: sw_hll
1791: args: -sw_riemann hll -bc_wall 1,2,3,4 -ufv_cfl 3 -petscfv_type leastsquares -petsclimiter_type sin \
1792: -grid_bounds 0,5,0,5 -dm_plex_simplex 0 -dm_plex_box_faces 25,25 \
1793: -ts_max_steps 5 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1794: -monitor height,energy
1796: # 2D Euler
1797: testset:
1798: args: -physics euler -eu_type linear_wave -eu_gamma 1.4 -dm_plex_adj_cone -dm_plex_adj_closure 0 \
1799: -ufv_vtk_interval 0 -ufv_vtk_basename ${wPETSC_DIR}/ex11 -monitor density,energy
1801: test:
1802: suffix: euler_0
1803: requires: exodusii !complex !single
1804: args: -eu_riemann godunov -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1805: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo \
1806: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1808: test:
1809: suffix: euler_ceed
1810: requires: exodusii libceed
1811: args: -eu_riemann godunov_ceed -bc_wall 100,101 -ufv_cfl 5 -petsclimiter_type sin \
1812: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/annulus-20.exo -dm_plex_use_ceed \
1813: -ts_max_time 1 -ts_ssp_type rks2 -ts_ssp_nstages 10
1815: testset:
1816: args: -dm_plex_adj_cone -dm_plex_adj_closure 0 -dm_plex_simplex 0 -dm_plex_box_faces 1,1,1
1818: # 2D Advection: p4est
1819: test:
1820: suffix: p4est_advec_2d
1821: requires: p4est
1822: 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
1824: # Advection in a box
1825: test:
1826: suffix: adv_2d_quad_0
1827: args: -ufv_vtk_interval 0 -dm_refine 3 -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1829: test:
1830: suffix: adv_2d_quad_1
1831: 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
1832: timeoutfactor: 3
1834: test:
1835: suffix: adv_2d_quad_p4est_0
1836: requires: p4est
1837: args: -ufv_vtk_interval 0 -dm_refine 5 -dm_type p4est -dm_plex_separate_marker -bc_inflow 1,2,4 -bc_outflow 3
1839: test:
1840: suffix: adv_2d_quad_p4est_1
1841: requires: p4est
1842: 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
1843: timeoutfactor: 3
1845: test: # broken for quad precision
1846: suffix: adv_2d_quad_p4est_adapt_0
1847: requires: p4est !__float128
1848: 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
1849: timeoutfactor: 3
1851: test:
1852: suffix: adv_0
1853: requires: exodusii
1854: 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
1856: # Run with -dm_forest_maximum_refinement 6 -ts_max_time 0.5 instead to get the full movie
1857: test:
1858: suffix: shock_0
1859: requires: p4est !single !complex
1860: args: -dm_plex_box_faces 2,1 -grid_bounds -1,1.,0.,1 -grid_skew_60 \
1861: -dm_type p4est -dm_forest_partition_overlap 1 -dm_forest_maximum_refinement 2 -dm_forest_minimum_refinement 2 -dm_forest_initial_refinement 2 \
1862: -ufv_use_amr -refine_vec_tagger_box 0.5,inf -coarsen_vec_tagger_box 0,1.e-2 -refine_tag_view -coarsen_tag_view \
1863: -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. \
1864: -petscfv_type leastsquares -petsclimiter_type minmod -petscfv_compute_gradients 0 \
1865: -ts_max_steps 3 -ts_ssp_type rks2 -ts_ssp_nstages 10 \
1866: -ufv_vtk_basename ${wPETSC_DIR}/ex11 -ufv_vtk_interval 0 -monitor density,energy
1868: # Test GLVis visualization of PetscFV fields
1869: test:
1870: suffix: glvis_adv_2d_tri
1871: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1872: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1873: -ts_monitor_solution glvis: -ts_max_steps 0
1875: test:
1876: suffix: glvis_adv_2d_quad
1877: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
1878: -dm_refine 5 -dm_plex_separate_marker \
1879: -ts_monitor_solution glvis: -ts_max_steps 0
1881: # Test CGNS file writing for PetscFV fields
1882: test:
1883: suffix: cgns_adv_2d_tri
1884: requires: cgns
1885: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1886: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1887: -ts_monitor_solution cgns:sol.cgns -ts_max_steps 0
1889: test:
1890: suffix: cgns_adv_2d_quad
1891: requires: cgns
1892: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 -bc_inflow 1,2,4 -bc_outflow 3 \
1893: -dm_refine 5 -dm_plex_separate_marker \
1894: -ts_monitor_solution cgns:sol.cgns -ts_max_steps 0
1896: # Test VTK file writing for PetscFV fields with -ts_monitor_solution_vtk
1897: test:
1898: suffix: vtk_adv_2d_tri
1899: args: -ufv_vtk_interval 0 -ufv_vtk_monitor 0 \
1900: -dm_plex_filename ${wPETSC_DIR}/share/petsc/datafiles/meshes/square_periodic.msh -dm_plex_gmsh_periodic 0 \
1901: -ts_monitor_solution_vtk 'bar-%03d.vtu' -ts_monitor_solution_vtk_interval 2 -ts_max_steps 5
1903: TEST*/