Actual source code: ex11_sa.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.
36: F*/
37: #include <petscts.h>
38: #include <petscfv.h>
39: #include <petscdmplex.h>
40: #include <petscsf.h>
41: #include <petscblaslapack.h>
43: #define DIM 2 /* Geometric dimension */
45: static PetscFunctionList PhysicsList;
47: /* Represents continuum physical equations. */
48: typedef struct _n_Physics *Physics;
50: /* Physical model includes boundary conditions, initial conditions, and functionals of interest. It is
51: * discretization-independent, but its members depend on the scenario being solved. */
52: typedef struct _n_Model *Model;
54: /* 'User' implements a discretization of a continuous model. */
55: typedef struct _n_User *User;
57: typedef PetscErrorCode (*RiemannFunction)(const PetscReal *, const PetscReal *, const PetscScalar *, const PetscScalar *, PetscScalar *, void *);
58: typedef PetscErrorCode (*SolutionFunction)(Model, PetscReal, const PetscReal *, PetscScalar *, void *);
59: typedef PetscErrorCode (*FunctionalFunction)(Model, PetscReal, const PetscReal *, const PetscScalar *, PetscReal *, void *);
60: typedef PetscErrorCode (*SetupFields)(Physics, PetscSection);
61: static PetscErrorCode ModelSolutionSetDefault(Model, SolutionFunction, void *);
62: static PetscErrorCode ModelFunctionalRegister(Model, const char *, PetscInt *, FunctionalFunction, void *);
63: static PetscErrorCode OutputVTK(DM, const char *, PetscViewer *);
65: struct FieldDescription {
66: const char *name;
67: PetscInt dof;
68: };
70: typedef struct _n_FunctionalLink *FunctionalLink;
71: struct _n_FunctionalLink {
72: char *name;
73: FunctionalFunction func;
74: void *ctx;
75: PetscInt offset;
76: FunctionalLink next;
77: };
79: struct _n_Physics {
80: RiemannFunction riemann;
81: PetscInt dof; /* number of degrees of freedom per cell */
82: PetscReal maxspeed; /* kludge to pick initial time step, need to add monitoring and step control */
83: void *data;
84: PetscInt nfields;
85: const struct FieldDescription *field_desc;
86: };
88: struct _n_Model {
89: MPI_Comm comm; /* Does not do collective communication, but some error conditions can be collective */
90: Physics physics;
91: FunctionalLink functionalRegistry;
92: PetscInt maxComputed;
93: PetscInt numMonitored;
94: FunctionalLink *functionalMonitored;
95: PetscInt numCall;
96: FunctionalLink *functionalCall;
97: SolutionFunction solution;
98: void *solutionctx;
99: PetscReal maxspeed; /* estimate of global maximum speed (for CFL calculation) */
100: };
102: struct _n_User {
103: PetscInt numSplitFaces;
104: PetscInt vtkInterval; /* For monitor */
105: Model model;
106: };
108: static inline PetscScalar DotDIM(const PetscScalar *x, const PetscScalar *y)
109: {
110: PetscInt i;
111: PetscScalar prod = 0.0;
113: for (i = 0; i < DIM; i++) prod += x[i] * y[i];
114: return prod;
115: }
116: static inline PetscReal NormDIM(const PetscScalar *x)
117: {
118: return PetscSqrtReal(PetscAbsScalar(DotDIM(x, x)));
119: }
120: static inline void axDIM(const PetscScalar a, PetscScalar *x)
121: {
122: PetscInt i;
123: for (i = 0; i < DIM; i++) x[i] *= a;
124: }
125: static inline void waxDIM(const PetscScalar a, const PetscScalar *x, PetscScalar *w)
126: {
127: PetscInt i;
128: for (i = 0; i < DIM; i++) w[i] = x[i] * a;
129: }
130: static inline void NormalSplitDIM(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
131: { /* Split x into normal and tangential components */
132: PetscInt i;
133: PetscScalar c;
134: c = DotDIM(x, n) / DotDIM(n, n);
135: for (i = 0; i < DIM; i++) {
136: xn[i] = c * n[i];
137: xt[i] = x[i] - xn[i];
138: }
139: }
141: static inline PetscScalar Dot2(const PetscScalar *x, const PetscScalar *y)
142: {
143: return x[0] * y[0] + x[1] * y[1];
144: }
145: static inline PetscReal Norm2(const PetscScalar *x)
146: {
147: return PetscSqrtReal(PetscAbsScalar(Dot2(x, x)));
148: }
149: static inline void Normalize2(PetscScalar *x)
150: {
151: PetscReal a = 1. / Norm2(x);
152: x[0] *= a;
153: x[1] *= a;
154: }
155: static inline void Waxpy2(PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
156: {
157: w[0] = a * x[0] + y[0];
158: w[1] = a * x[1] + y[1];
159: }
160: static inline void Scale2(PetscScalar a, const PetscScalar *x, PetscScalar *y)
161: {
162: y[0] = a * x[0];
163: y[1] = a * x[1];
164: }
166: static inline void WaxpyD(PetscInt dim, PetscScalar a, const PetscScalar *x, const PetscScalar *y, PetscScalar *w)
167: {
168: PetscInt d;
169: for (d = 0; d < dim; ++d) w[d] = a * x[d] + y[d];
170: }
171: static inline PetscScalar DotD(PetscInt dim, const PetscScalar *x, const PetscScalar *y)
172: {
173: PetscScalar sum = 0.0;
174: PetscInt d;
175: for (d = 0; d < dim; ++d) sum += x[d] * y[d];
176: return sum;
177: }
178: static inline PetscReal NormD(PetscInt dim, const PetscScalar *x)
179: {
180: return PetscSqrtReal(PetscAbsScalar(DotD(dim, x, x)));
181: }
183: static inline void NormalSplit(const PetscReal *n, const PetscScalar *x, PetscScalar *xn, PetscScalar *xt)
184: { /* Split x into normal and tangential components */
185: Scale2(Dot2(x, n) / Dot2(n, n), n, xn);
186: Waxpy2(-1, xn, x, xt);
187: }
189: /******************* Advect ********************/
190: typedef enum {
191: ADVECT_SOL_TILTED,
192: ADVECT_SOL_BUMP
193: } AdvectSolType;
194: static const char *const AdvectSolTypes[] = {"TILTED", "BUMP", "AdvectSolType", "ADVECT_SOL_", 0};
195: typedef enum {
196: ADVECT_SOL_BUMP_CONE,
197: ADVECT_SOL_BUMP_COS
198: } AdvectSolBumpType;
199: static const char *const AdvectSolBumpTypes[] = {"CONE", "COS", "AdvectSolBumpType", "ADVECT_SOL_BUMP_", 0};
201: typedef struct {
202: PetscReal wind[DIM];
203: } Physics_Advect_Tilted;
204: typedef struct {
205: PetscReal center[DIM];
206: PetscReal radius;
207: AdvectSolBumpType type;
208: } Physics_Advect_Bump;
210: typedef struct {
211: PetscReal inflowState;
212: AdvectSolType soltype;
213: union
214: {
215: Physics_Advect_Tilted tilted;
216: Physics_Advect_Bump bump;
217: } sol;
218: struct {
219: PetscInt Error;
220: } functional;
221: } Physics_Advect;
223: static const struct FieldDescription PhysicsFields_Advect[] = {
224: {"U", 1},
225: {NULL, 0}
226: };
228: static PetscErrorCode PhysicsBoundary_Advect_Inflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
229: {
230: Physics phys = (Physics)ctx;
231: Physics_Advect *advect = (Physics_Advect *)phys->data;
233: PetscFunctionBeginUser;
234: xG[0] = advect->inflowState;
235: PetscFunctionReturn(PETSC_SUCCESS);
236: }
238: static PetscErrorCode PhysicsBoundary_Advect_Outflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
239: {
240: PetscFunctionBeginUser;
241: xG[0] = xI[0];
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: static PetscErrorCode PhysicsRiemann_Advect(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
246: {
247: Physics_Advect *advect = (Physics_Advect *)phys->data;
248: PetscReal wind[DIM], wn;
250: PetscFunctionBeginUser;
251: switch (advect->soltype) {
252: case ADVECT_SOL_TILTED: {
253: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
254: wind[0] = tilted->wind[0];
255: wind[1] = tilted->wind[1];
256: } break;
257: case ADVECT_SOL_BUMP:
258: wind[0] = -qp[1];
259: wind[1] = qp[0];
260: break;
261: default:
262: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No support for solution type %s", AdvectSolBumpTypes[advect->soltype]);
263: }
264: wn = Dot2(wind, n);
265: flux[0] = (wn > 0 ? xL[0] : xR[0]) * wn;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: }
269: static PetscErrorCode PhysicsSolution_Advect(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
270: {
271: Physics phys = (Physics)ctx;
272: Physics_Advect *advect = (Physics_Advect *)phys->data;
274: PetscFunctionBeginUser;
275: switch (advect->soltype) {
276: case ADVECT_SOL_TILTED: {
277: PetscReal x0[DIM];
278: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
279: Waxpy2(-time, tilted->wind, x, x0);
280: if (x0[1] > 0) u[0] = 1. * x[0] + 3. * x[1];
281: else u[0] = advect->inflowState;
282: } break;
283: case ADVECT_SOL_BUMP: {
284: Physics_Advect_Bump *bump = &advect->sol.bump;
285: PetscReal x0[DIM], v[DIM], r, cost, sint;
286: cost = PetscCosReal(time);
287: sint = PetscSinReal(time);
288: x0[0] = cost * x[0] + sint * x[1];
289: x0[1] = -sint * x[0] + cost * x[1];
290: Waxpy2(-1, bump->center, x0, v);
291: r = Norm2(v);
292: switch (bump->type) {
293: case ADVECT_SOL_BUMP_CONE:
294: u[0] = PetscMax(1 - r / bump->radius, 0);
295: break;
296: case ADVECT_SOL_BUMP_COS:
297: u[0] = 0.5 + 0.5 * PetscCosReal(PetscMin(r / bump->radius, 1) * PETSC_PI);
298: break;
299: }
300: } break;
301: default:
302: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Unknown solution type");
303: }
304: PetscFunctionReturn(PETSC_SUCCESS);
305: }
307: static PetscErrorCode PhysicsFunctional_Advect(Model mod, PetscReal time, const PetscScalar *x, const PetscScalar *y, PetscReal *f, void *ctx)
308: {
309: Physics phys = (Physics)ctx;
310: Physics_Advect *advect = (Physics_Advect *)phys->data;
311: PetscScalar yexact[1];
313: PetscFunctionBeginUser;
314: PetscCall(PhysicsSolution_Advect(mod, time, x, yexact, phys));
315: f[advect->functional.Error] = PetscAbsScalar(y[0] - yexact[0]);
316: PetscFunctionReturn(PETSC_SUCCESS);
317: }
319: static PetscErrorCode PhysicsCreate_Advect(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
320: {
321: Physics_Advect *advect;
323: PetscFunctionBeginUser;
324: phys->field_desc = PhysicsFields_Advect;
325: phys->riemann = (RiemannFunction)PhysicsRiemann_Advect;
326: PetscCall(PetscNew(&advect));
327: phys->data = advect;
328: PetscOptionsHeadBegin(PetscOptionsObject, "Advect options");
329: {
330: PetscInt two = 2, dof = 1;
331: advect->soltype = ADVECT_SOL_TILTED;
332: PetscCall(PetscOptionsEnum("-advect_sol_type", "solution type", "", AdvectSolTypes, (PetscEnum)advect->soltype, (PetscEnum *)&advect->soltype, NULL));
333: switch (advect->soltype) {
334: case ADVECT_SOL_TILTED: {
335: Physics_Advect_Tilted *tilted = &advect->sol.tilted;
336: two = 2;
337: tilted->wind[0] = 0.0;
338: tilted->wind[1] = 1.0;
339: PetscCall(PetscOptionsRealArray("-advect_tilted_wind", "background wind vx,vy", "", tilted->wind, &two, NULL));
340: advect->inflowState = -2.0;
341: PetscCall(PetscOptionsRealArray("-advect_tilted_inflow", "Inflow state", "", &advect->inflowState, &dof, NULL));
342: phys->maxspeed = Norm2(tilted->wind);
343: } break;
344: case ADVECT_SOL_BUMP: {
345: Physics_Advect_Bump *bump = &advect->sol.bump;
346: two = 2;
347: bump->center[0] = 2.;
348: bump->center[1] = 0.;
349: PetscCall(PetscOptionsRealArray("-advect_bump_center", "location of center of bump x,y", "", bump->center, &two, NULL));
350: bump->radius = 0.9;
351: PetscCall(PetscOptionsReal("-advect_bump_radius", "radius of bump", "", bump->radius, &bump->radius, NULL));
352: bump->type = ADVECT_SOL_BUMP_CONE;
353: PetscCall(PetscOptionsEnum("-advect_bump_type", "type of bump", "", AdvectSolBumpTypes, (PetscEnum)bump->type, (PetscEnum *)&bump->type, NULL));
354: phys->maxspeed = 3.; /* radius of mesh, kludge */
355: } break;
356: }
357: }
358: PetscOptionsHeadEnd();
359: {
360: const PetscInt inflowids[] = {100, 200, 300}, outflowids[] = {101};
361: DMLabel label;
363: PetscCall(DMGetLabel(dm, "Face Sets", &label));
364: /* Register "canned" boundary conditions and defaults for where to apply. */
365: PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Inflow, NULL, phys, NULL));
366: PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)())PhysicsBoundary_Advect_Outflow, NULL, phys, NULL));
367: /* Initial/transient solution with default boundary conditions */
368: PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Advect, phys));
369: /* Register "canned" functionals */
370: PetscCall(ModelFunctionalRegister(mod, "Error", &advect->functional.Error, PhysicsFunctional_Advect, phys));
371: }
372: PetscFunctionReturn(PETSC_SUCCESS);
373: }
375: /******************* Shallow Water ********************/
376: typedef struct {
377: PetscReal gravity;
378: PetscReal boundaryHeight;
379: struct {
380: PetscInt Height;
381: PetscInt Speed;
382: PetscInt Energy;
383: } functional;
384: } Physics_SW;
385: typedef struct {
386: PetscScalar vals[0];
387: PetscScalar h;
388: PetscScalar uh[DIM];
389: } SWNode;
391: static const struct FieldDescription PhysicsFields_SW[] = {
392: {"Height", 1 },
393: {"Momentum", DIM},
394: {NULL, 0 }
395: };
397: /*
398: * h_t + div(uh) = 0
399: * (uh)_t + div (u\otimes uh + g h^2 / 2 I) = 0
400: *
401: * */
402: static PetscErrorCode SWFlux(Physics phys, const PetscReal *n, const SWNode *x, SWNode *f)
403: {
404: Physics_SW *sw = (Physics_SW *)phys->data;
405: PetscScalar uhn, u[DIM];
406: PetscInt i;
408: PetscFunctionBeginUser;
409: Scale2(1. / x->h, x->uh, u);
410: uhn = Dot2(x->uh, n);
411: f->h = uhn;
412: for (i = 0; i < DIM; i++) f->uh[i] = u[i] * uhn + sw->gravity * PetscSqr(x->h) * n[i];
413: PetscFunctionReturn(PETSC_SUCCESS);
414: }
416: static PetscErrorCode PhysicsBoundary_SW_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
417: {
418: PetscFunctionBeginUser;
419: xG[0] = xI[0];
420: xG[1] = -xI[1];
421: xG[2] = -xI[2];
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: static PetscErrorCode PhysicsRiemann_SW(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
426: {
427: Physics_SW *sw = (Physics_SW *)phys->data;
428: PetscReal cL, cR, speed, nn[DIM];
429: const SWNode *uL = (const SWNode *)xL, *uR = (const SWNode *)xR;
430: SWNode fL, fR;
431: PetscInt i;
433: PetscFunctionBeginUser;
434: PetscCheck(uL->h >= 0 && uR->h >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
435: nn[0] = n[0];
436: nn[1] = n[1];
437: Normalize2(nn);
438: SWFlux(phys, nn, uL, &fL);
439: SWFlux(phys, nn, uR, &fR);
440: cL = PetscSqrtReal(sw->gravity * PetscRealPart(uL->h));
441: cR = PetscSqrtReal(sw->gravity * PetscRealPart(uR->h)); /* gravity wave speed */
442: speed = PetscMax(PetscAbsScalar(Dot2(uL->uh, nn) / uL->h) + cL, PetscAbsScalar(Dot2(uR->uh, nn) / uR->h) + cR);
443: for (i = 0; i < 1 + DIM; i++) flux[i] = (0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i])) * Norm2(n);
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode PhysicsSolution_SW(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
448: {
449: PetscReal dx[2], r, sigma;
451: PetscFunctionBeginUser;
452: PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
453: dx[0] = x[0] - 1.5;
454: dx[1] = x[1] - 1.0;
455: r = Norm2(dx);
456: sigma = 0.5;
457: u[0] = 1 + 2 * PetscExpScalar(-PetscSqr(r) / (2 * PetscSqr(sigma)));
458: u[1] = 0.0;
459: u[2] = 0.0;
460: PetscFunctionReturn(PETSC_SUCCESS);
461: }
463: static PetscErrorCode PhysicsFunctional_SW(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
464: {
465: Physics phys = (Physics)ctx;
466: Physics_SW *sw = (Physics_SW *)phys->data;
467: const SWNode *x = (const SWNode *)xx;
468: PetscScalar u[2];
469: PetscReal h;
471: PetscFunctionBeginUser;
472: h = PetscRealPart(x->h);
473: Scale2(1. / x->h, x->uh, u);
474: f[sw->functional.Height] = h;
475: f[sw->functional.Speed] = Norm2(u) + PetscSqrtReal(sw->gravity * h);
476: f[sw->functional.Energy] = 0.5 * (Dot2(x->uh, u) + sw->gravity * PetscSqr(h));
477: PetscFunctionReturn(PETSC_SUCCESS);
478: }
480: static PetscErrorCode PhysicsCreate_SW(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
481: {
482: Physics_SW *sw;
484: PetscFunctionBeginUser;
485: phys->field_desc = PhysicsFields_SW;
486: phys->riemann = (RiemannFunction)PhysicsRiemann_SW;
487: PetscCall(PetscNew(&sw));
488: phys->data = sw;
489: PetscOptionsHeadBegin(PetscOptionsObject, "SW options");
490: {
491: sw->gravity = 1.0;
492: PetscCall(PetscOptionsReal("-sw_gravity", "Gravitational constant", "", sw->gravity, &sw->gravity, NULL));
493: }
494: PetscOptionsHeadEnd();
495: phys->maxspeed = PetscSqrtReal(2.0 * sw->gravity); /* Mach 1 for depth of 2 */
497: {
498: const PetscInt wallids[] = {100, 101, 200, 300};
499: DMLabel label;
501: PetscCall(DMGetLabel(dm, "Face Sets", &label));
502: PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_SW_Wall, NULL, phys, NULL));
503: PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_SW, phys));
504: PetscCall(ModelFunctionalRegister(mod, "Height", &sw->functional.Height, PhysicsFunctional_SW, phys));
505: PetscCall(ModelFunctionalRegister(mod, "Speed", &sw->functional.Speed, PhysicsFunctional_SW, phys));
506: PetscCall(ModelFunctionalRegister(mod, "Energy", &sw->functional.Energy, PhysicsFunctional_SW, phys));
507: }
508: PetscFunctionReturn(PETSC_SUCCESS);
509: }
511: /******************* Euler ********************/
512: typedef struct {
513: PetscScalar vals[0];
514: PetscScalar r;
515: PetscScalar ru[DIM];
516: PetscScalar e;
517: } EulerNode;
518: typedef PetscErrorCode (*EquationOfState)(const PetscReal *, const EulerNode *, PetscScalar *);
519: typedef struct {
520: PetscInt npars;
521: PetscReal pars[DIM];
522: EquationOfState pressure;
523: EquationOfState sound;
524: struct {
525: PetscInt Density;
526: PetscInt Momentum;
527: PetscInt Energy;
528: PetscInt Pressure;
529: PetscInt Speed;
530: } monitor;
531: } Physics_Euler;
533: static const struct FieldDescription PhysicsFields_Euler[] = {
534: {"Density", 1 },
535: {"Momentum", DIM},
536: {"Energy", 1 },
537: {NULL, 0 }
538: };
540: static PetscErrorCode Pressure_PG(const PetscReal *pars, const EulerNode *x, PetscScalar *p)
541: {
542: PetscScalar ru2;
544: PetscFunctionBeginUser;
545: ru2 = DotDIM(x->ru, x->ru);
546: ru2 /= x->r;
547: /* kinematic dof = params[0] */
548: (*p) = 2.0 * (x->e - 0.5 * ru2) / pars[0];
549: PetscFunctionReturn(PETSC_SUCCESS);
550: }
552: static PetscErrorCode SpeedOfSound_PG(const PetscReal *pars, const EulerNode *x, PetscScalar *c)
553: {
554: PetscScalar p;
556: PetscFunctionBeginUser;
557: /* TODO remove direct usage of Pressure_PG */
558: Pressure_PG(pars, x, &p);
559: /* TODO check the sign of p */
560: /* pars[1] = heat capacity ratio */
561: (*c) = PetscSqrtScalar(pars[1] * p / x->r);
562: PetscFunctionReturn(PETSC_SUCCESS);
563: }
565: /*
566: * x = (rho,rho*(u_1),...,rho*e)^T
567: * x_t+div(f_1(x))+...+div(f_DIM(x)) = 0
568: *
569: * f_i(x) = u_i*x+(0,0,...,p,...,p*u_i)^T
570: *
571: * */
572: static PetscErrorCode EulerFlux(Physics phys, const PetscReal *n, const EulerNode *x, EulerNode *f)
573: {
574: Physics_Euler *eu = (Physics_Euler *)phys->data;
575: PetscScalar u, nu, p;
576: PetscInt i;
578: PetscFunctionBeginUser;
579: u = DotDIM(x->ru, x->ru);
580: u /= (x->r * x->r);
581: nu = DotDIM(x->ru, n);
582: /* TODO check the sign of p */
583: eu->pressure(eu->pars, x, &p);
584: f->r = nu * x->r;
585: for (i = 0; i < DIM; i++) f->ru[i] = nu * x->ru[i] + n[i] * p;
586: f->e = nu * (x->e + p);
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: /* PetscReal* => EulerNode* conversion */
591: static PetscErrorCode PhysicsBoundary_Euler_Wall(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *xI, PetscScalar *xG, void *ctx)
592: {
593: PetscInt i;
594: PetscScalar xn[DIM], xt[DIM];
596: PetscFunctionBeginUser;
597: xG[0] = xI[0];
598: NormalSplitDIM(n, xI + 1, xn, xt);
599: for (i = 0; i < DIM; i++) xG[i + 1] = -xn[i] + xt[i];
600: xG[DIM + 1] = xI[DIM + 1];
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: /* PetscReal* => EulerNode* conversion */
605: static PetscErrorCode PhysicsRiemann_Euler_Rusanov(const PetscReal *qp, const PetscReal *n, const PetscScalar *xL, const PetscScalar *xR, PetscScalar *flux, Physics phys)
606: {
607: Physics_Euler *eu = (Physics_Euler *)phys->data;
608: PetscScalar cL, cR, speed;
609: const EulerNode *uL = (const EulerNode *)xL, *uR = (const EulerNode *)xR;
610: EulerNode fL, fR;
611: PetscInt i;
613: PetscFunctionBeginUser;
614: PetscCheck(uL->r >= 0 && uR->r >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
615: EulerFlux(phys, n, uL, &fL);
616: EulerFlux(phys, n, uR, &fR);
617: eu->sound(eu->pars, uL, &cL);
618: eu->sound(eu->pars, uR, &cR);
619: speed = PetscMax(cL, cR) + PetscMax(PetscAbsScalar(DotDIM(uL->ru, n) / NormDIM(n)), PetscAbsScalar(DotDIM(uR->ru, n) / NormDIM(n)));
620: for (i = 0; i < 2 + DIM; i++) flux[i] = 0.5 * (fL.vals[i] + fR.vals[i]) + 0.5 * speed * (xL[i] - xR[i]);
621: PetscFunctionReturn(PETSC_SUCCESS);
622: }
624: static PetscErrorCode PhysicsSolution_Euler(Model mod, PetscReal time, const PetscReal *x, PetscScalar *u, void *ctx)
625: {
626: PetscInt i;
628: PetscFunctionBeginUser;
629: PetscCheck(time == 0.0, mod->comm, PETSC_ERR_SUP, "No solution known for time %g", (double)time);
630: u[0] = 1.0;
631: u[DIM + 1] = 1.0 + PetscAbsReal(x[0]);
632: for (i = 1; i < DIM + 1; i++) u[i] = 0.0;
633: PetscFunctionReturn(PETSC_SUCCESS);
634: }
636: static PetscErrorCode PhysicsFunctional_Euler(Model mod, PetscReal time, const PetscReal *coord, const PetscScalar *xx, PetscReal *f, void *ctx)
637: {
638: Physics phys = (Physics)ctx;
639: Physics_Euler *eu = (Physics_Euler *)phys->data;
640: const EulerNode *x = (const EulerNode *)xx;
641: PetscScalar p;
643: PetscFunctionBeginUser;
644: f[eu->monitor.Density] = x->r;
645: f[eu->monitor.Momentum] = NormDIM(x->ru);
646: f[eu->monitor.Energy] = x->e;
647: f[eu->monitor.Speed] = NormDIM(x->ru) / x->r;
648: eu->pressure(eu->pars, x, &p);
649: f[eu->monitor.Pressure] = p;
650: PetscFunctionReturn(PETSC_SUCCESS);
651: }
653: static PetscErrorCode PhysicsCreate_Euler(PetscDS prob, Model mod, Physics phys, PetscOptions *PetscOptionsObject)
654: {
655: Physics_Euler *eu;
657: PetscFunctionBeginUser;
658: phys->field_desc = PhysicsFields_Euler;
659: phys->riemann = (RiemannFunction)PhysicsRiemann_Euler_Rusanov;
660: PetscCall(PetscNew(&eu));
661: phys->data = eu;
662: PetscOptionsHeadBegin(PetscOptionsObject, "Euler options");
663: {
664: eu->pars[0] = 3.0;
665: eu->pars[1] = 1.67;
666: PetscCall(PetscOptionsReal("-eu_f", "Degrees of freedom", "", eu->pars[0], &eu->pars[0], NULL));
667: PetscCall(PetscOptionsReal("-eu_gamma", "Heat capacity ratio", "", eu->pars[1], &eu->pars[1], NULL));
668: }
669: PetscOptionsHeadEnd();
670: eu->pressure = Pressure_PG;
671: eu->sound = SpeedOfSound_PG;
672: phys->maxspeed = 1.0;
673: {
674: const PetscInt wallids[] = {100, 101, 200, 300};
675: DMLabel label;
677: PetscCall(DMGetLabel(dm, "Face Sets", &label));
678: PetscCall(PetscDSAddBoundary(prob, PETSC_TRUE, "wall", label, PETSC_STATIC_ARRAY_LENGTH(wallids), wallids, 0, 0, NULL, (void (*)())PhysicsBoundary_Euler_Wall, NULL, phys, NULL));
679: PetscCall(ModelSolutionSetDefault(mod, PhysicsSolution_Euler, phys));
680: PetscCall(ModelFunctionalRegister(mod, "Speed", &eu->monitor.Speed, PhysicsFunctional_Euler, phys));
681: PetscCall(ModelFunctionalRegister(mod, "Energy", &eu->monitor.Energy, PhysicsFunctional_Euler, phys));
682: PetscCall(ModelFunctionalRegister(mod, "Density", &eu->monitor.Density, PhysicsFunctional_Euler, phys));
683: PetscCall(ModelFunctionalRegister(mod, "Momentum", &eu->monitor.Momentum, PhysicsFunctional_Euler, phys));
684: PetscCall(ModelFunctionalRegister(mod, "Pressure", &eu->monitor.Pressure, PhysicsFunctional_Euler, phys));
685: }
686: PetscFunctionReturn(PETSC_SUCCESS);
687: }
689: PetscErrorCode ConstructCellBoundary(DM dm, User user)
690: {
691: const char *name = "Cell Sets";
692: const char *bdname = "split faces";
693: IS regionIS, innerIS;
694: const PetscInt *regions, *cells;
695: PetscInt numRegions, innerRegion, numCells, c;
696: PetscInt cStart, cEnd, cEndInterior, fStart, fEnd;
697: PetscBool hasLabel;
699: PetscFunctionBeginUser;
700: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
701: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
702: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
704: PetscCall(DMHasLabel(dm, name, &hasLabel));
705: if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
706: PetscCall(DMGetLabelSize(dm, name, &numRegions));
707: if (numRegions != 2) PetscFunctionReturn(PETSC_SUCCESS);
708: /* Get the inner id */
709: PetscCall(DMGetLabelIdIS(dm, name, ®ionIS));
710: PetscCall(ISGetIndices(regionIS, ®ions));
711: innerRegion = regions[0];
712: PetscCall(ISRestoreIndices(regionIS, ®ions));
713: PetscCall(ISDestroy(®ionIS));
714: /* Find the faces between cells in different regions, could call DMPlexCreateNeighborCSR() */
715: PetscCall(DMGetStratumIS(dm, name, innerRegion, &innerIS));
716: PetscCall(ISGetLocalSize(innerIS, &numCells));
717: PetscCall(ISGetIndices(innerIS, &cells));
718: PetscCall(DMCreateLabel(dm, bdname));
719: for (c = 0; c < numCells; ++c) {
720: const PetscInt cell = cells[c];
721: const PetscInt *faces;
722: PetscInt numFaces, f;
724: PetscCheck((cell >= cStart) && (cell < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", cell);
725: PetscCall(DMPlexGetConeSize(dm, cell, &numFaces));
726: PetscCall(DMPlexGetCone(dm, cell, &faces));
727: for (f = 0; f < numFaces; ++f) {
728: const PetscInt face = faces[f];
729: const PetscInt *neighbors;
730: PetscInt nC, regionA, regionB;
732: PetscCheck((face >= fStart) && (face < fEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a face", face);
733: PetscCall(DMPlexGetSupportSize(dm, face, &nC));
734: if (nC != 2) continue;
735: PetscCall(DMPlexGetSupport(dm, face, &neighbors));
736: if ((neighbors[0] >= cEndInterior) || (neighbors[1] >= cEndInterior)) continue;
737: PetscCheck((neighbors[0] >= cStart) && (neighbors[0] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[0]);
738: PetscCheck((neighbors[1] >= cStart) && (neighbors[1] < cEnd), PETSC_COMM_SELF, PETSC_ERR_LIB, "Got invalid point %d which is not a cell", neighbors[1]);
739: PetscCall(DMGetLabelValue(dm, name, neighbors[0], ®ionA));
740: PetscCall(DMGetLabelValue(dm, name, neighbors[1], ®ionB));
741: PetscCheck(regionA >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[0]);
742: PetscCheck(regionB >= 0, PetscObjectComm((PetscObject)dm), PETSC_ERR_ARG_WRONG, "Invalid label %s: Cell %d has no value", name, neighbors[1]);
743: if (regionA != regionB) PetscCall(DMSetLabelValue(dm, bdname, faces[f], 1));
744: }
745: }
746: PetscCall(ISRestoreIndices(innerIS, &cells));
747: PetscCall(ISDestroy(&innerIS));
748: {
749: DMLabel label;
751: PetscCall(DMGetLabel(dm, bdname, &label));
752: PetscCall(DMLabelView(label, PETSC_VIEWER_STDOUT_WORLD));
753: }
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: /* Right now, I have just added duplicate faces, which see both cells. We can
758: - Add duplicate vertices and decouple the face cones
759: - Disconnect faces from cells across the rotation gap
760: */
761: PetscErrorCode SplitFaces(DM *dmSplit, const char labelName[], User user)
762: {
763: DM dm = *dmSplit, sdm;
764: PetscSF sfPoint, gsfPoint;
765: PetscSection coordSection, newCoordSection;
766: Vec coordinates;
767: IS idIS;
768: const PetscInt *ids;
769: PetscInt *newpoints;
770: PetscInt dim, depth, maxConeSize, maxSupportSize, numLabels, numGhostCells;
771: PetscInt numFS, fs, pStart, pEnd, p, cEnd, cEndInterior, vStart, vEnd, v, fStart, fEnd, newf, d, l;
772: PetscBool hasLabel;
774: PetscFunctionBeginUser;
775: PetscCall(DMHasLabel(dm, labelName, &hasLabel));
776: if (!hasLabel) PetscFunctionReturn(PETSC_SUCCESS);
777: PetscCall(DMCreate(PetscObjectComm((PetscObject)dm), &sdm));
778: PetscCall(DMSetType(sdm, DMPLEX));
779: PetscCall(DMGetDimension(dm, &dim));
780: PetscCall(DMSetDimension(sdm, dim));
782: PetscCall(DMGetLabelIdIS(dm, labelName, &idIS));
783: PetscCall(ISGetLocalSize(idIS, &numFS));
784: PetscCall(ISGetIndices(idIS, &ids));
786: user->numSplitFaces = 0;
787: for (fs = 0; fs < numFS; ++fs) {
788: PetscInt numBdFaces;
790: PetscCall(DMGetStratumSize(dm, labelName, ids[fs], &numBdFaces));
791: user->numSplitFaces += numBdFaces;
792: }
793: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
794: pEnd += user->numSplitFaces;
795: PetscCall(DMPlexSetChart(sdm, pStart, pEnd));
796: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
797: PetscCall(DMPlexSetGhostCellStratum(sdm, cEndInterior, PETSC_DETERMINE));
798: PetscCall(DMPlexGetHeightStratum(dm, 0, NULL, &cEnd));
799: numGhostCells = cEnd - cEndInterior;
800: /* Set cone and support sizes */
801: PetscCall(DMPlexGetDepth(dm, &depth));
802: for (d = 0; d <= depth; ++d) {
803: PetscCall(DMPlexGetDepthStratum(dm, d, &pStart, &pEnd));
804: for (p = pStart; p < pEnd; ++p) {
805: PetscInt newp = p;
806: PetscInt size;
808: PetscCall(DMPlexGetConeSize(dm, p, &size));
809: PetscCall(DMPlexSetConeSize(sdm, newp, size));
810: PetscCall(DMPlexGetSupportSize(dm, p, &size));
811: PetscCall(DMPlexSetSupportSize(sdm, newp, size));
812: }
813: }
814: PetscCall(DMPlexGetHeightStratum(dm, 1, &fStart, &fEnd));
815: for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
816: IS faceIS;
817: const PetscInt *faces;
818: PetscInt numFaces, f;
820: PetscCall(DMGetStratumIS(dm, labelName, ids[fs], &faceIS));
821: PetscCall(ISGetLocalSize(faceIS, &numFaces));
822: PetscCall(ISGetIndices(faceIS, &faces));
823: for (f = 0; f < numFaces; ++f, ++newf) {
824: PetscInt size;
826: /* Right now I think that both faces should see both cells */
827: PetscCall(DMPlexGetConeSize(dm, faces[f], &size));
828: PetscCall(DMPlexSetConeSize(sdm, newf, size));
829: PetscCall(DMPlexGetSupportSize(dm, faces[f], &size));
830: PetscCall(DMPlexSetSupportSize(sdm, newf, size));
831: }
832: PetscCall(ISRestoreIndices(faceIS, &faces));
833: PetscCall(ISDestroy(&faceIS));
834: }
835: PetscCall(DMSetUp(sdm));
836: /* Set cones and supports */
837: PetscCall(DMPlexGetMaxSizes(dm, &maxConeSize, &maxSupportSize));
838: PetscCall(PetscMalloc1(PetscMax(maxConeSize, maxSupportSize), &newpoints));
839: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
840: for (p = pStart; p < pEnd; ++p) {
841: const PetscInt *points, *orientations;
842: PetscInt size, i, newp = p;
844: PetscCall(DMPlexGetConeSize(dm, p, &size));
845: PetscCall(DMPlexGetCone(dm, p, &points));
846: PetscCall(DMPlexGetConeOrientation(dm, p, &orientations));
847: for (i = 0; i < size; ++i) newpoints[i] = points[i];
848: PetscCall(DMPlexSetCone(sdm, newp, newpoints));
849: PetscCall(DMPlexSetConeOrientation(sdm, newp, orientations));
850: PetscCall(DMPlexGetSupportSize(dm, p, &size));
851: PetscCall(DMPlexGetSupport(dm, p, &points));
852: for (i = 0; i < size; ++i) newpoints[i] = points[i];
853: PetscCall(DMPlexSetSupport(sdm, newp, newpoints));
854: }
855: PetscCall(PetscFree(newpoints));
856: for (fs = 0, newf = fEnd; fs < numFS; ++fs) {
857: IS faceIS;
858: const PetscInt *faces;
859: PetscInt numFaces, f;
861: PetscCall(DMGetStratumIS(dm, labelName, ids[fs], &faceIS));
862: PetscCall(ISGetLocalSize(faceIS, &numFaces));
863: PetscCall(ISGetIndices(faceIS, &faces));
864: for (f = 0; f < numFaces; ++f, ++newf) {
865: const PetscInt *points;
867: PetscCall(DMPlexGetCone(dm, faces[f], &points));
868: PetscCall(DMPlexSetCone(sdm, newf, points));
869: PetscCall(DMPlexGetSupport(dm, faces[f], &points));
870: PetscCall(DMPlexSetSupport(sdm, newf, points));
871: }
872: PetscCall(ISRestoreIndices(faceIS, &faces));
873: PetscCall(ISDestroy(&faceIS));
874: }
875: PetscCall(ISRestoreIndices(idIS, &ids));
876: PetscCall(ISDestroy(&idIS));
877: PetscCall(DMPlexStratify(sdm));
878: /* Convert coordinates */
879: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
880: PetscCall(DMGetCoordinateSection(dm, &coordSection));
881: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &newCoordSection));
882: PetscCall(PetscSectionSetNumFields(newCoordSection, 1));
883: PetscCall(PetscSectionSetFieldComponents(newCoordSection, 0, dim));
884: PetscCall(PetscSectionSetChart(newCoordSection, vStart, vEnd));
885: for (v = vStart; v < vEnd; ++v) {
886: PetscCall(PetscSectionSetDof(newCoordSection, v, dim));
887: PetscCall(PetscSectionSetFieldDof(newCoordSection, v, 0, dim));
888: }
889: PetscCall(PetscSectionSetUp(newCoordSection));
890: PetscCall(DMSetCoordinateSection(sdm, PETSC_DETERMINE, newCoordSection));
891: PetscCall(PetscSectionDestroy(&newCoordSection)); /* relinquish our reference */
892: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
893: PetscCall(DMSetCoordinatesLocal(sdm, coordinates));
894: /* Convert labels */
895: PetscCall(DMGetNumLabels(dm, &numLabels));
896: for (l = 0; l < numLabels; ++l) {
897: const char *lname;
898: PetscBool isDepth;
900: PetscCall(DMGetLabelName(dm, l, &lname));
901: PetscCall(PetscStrcmp(lname, "depth", &isDepth));
902: if (isDepth) continue;
903: PetscCall(DMCreateLabel(sdm, lname));
904: PetscCall(DMGetLabelIdIS(dm, lname, &idIS));
905: PetscCall(ISGetLocalSize(idIS, &numFS));
906: PetscCall(ISGetIndices(idIS, &ids));
907: for (fs = 0; fs < numFS; ++fs) {
908: IS pointIS;
909: const PetscInt *points;
910: PetscInt numPoints;
912: PetscCall(DMGetStratumIS(dm, lname, ids[fs], &pointIS));
913: PetscCall(ISGetLocalSize(pointIS, &numPoints));
914: PetscCall(ISGetIndices(pointIS, &points));
915: for (p = 0; p < numPoints; ++p) {
916: PetscInt newpoint = points[p];
918: PetscCall(DMSetLabelValue(sdm, lname, newpoint, ids[fs]));
919: }
920: PetscCall(ISRestoreIndices(pointIS, &points));
921: PetscCall(ISDestroy(&pointIS));
922: }
923: PetscCall(ISRestoreIndices(idIS, &ids));
924: PetscCall(ISDestroy(&idIS));
925: }
926: /* Convert pointSF */
927: const PetscSFNode *remotePoints;
928: PetscSFNode *gremotePoints;
929: const PetscInt *localPoints;
930: PetscInt *glocalPoints, *newLocation, *newRemoteLocation;
931: PetscInt numRoots, numLeaves;
932: PetscMPIInt size;
934: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)dm), &size));
935: PetscCall(DMGetPointSF(dm, &sfPoint));
936: PetscCall(DMGetPointSF(sdm, &gsfPoint));
937: PetscCall(DMPlexGetChart(dm, &pStart, &pEnd));
938: PetscCall(PetscSFGetGraph(sfPoint, &numRoots, &numLeaves, &localPoints, &remotePoints));
939: if (numRoots >= 0) {
940: PetscCall(PetscMalloc2(numRoots, &newLocation, pEnd - pStart, &newRemoteLocation));
941: for (l = 0; l < numRoots; l++) newLocation[l] = l; /* + (l >= cEnd ? numGhostCells : 0); */
942: PetscCall(PetscSFBcastBegin(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
943: PetscCall(PetscSFBcastEnd(sfPoint, MPIU_INT, newLocation, newRemoteLocation, MPI_REPLACE));
944: PetscCall(PetscMalloc1(numLeaves, &glocalPoints));
945: PetscCall(PetscMalloc1(numLeaves, &gremotePoints));
946: for (l = 0; l < numLeaves; ++l) {
947: glocalPoints[l] = localPoints[l]; /* localPoints[l] >= cEnd ? localPoints[l] + numGhostCells : localPoints[l]; */
948: gremotePoints[l].rank = remotePoints[l].rank;
949: gremotePoints[l].index = newRemoteLocation[localPoints[l]];
950: }
951: PetscCall(PetscFree2(newLocation, newRemoteLocation));
952: PetscCall(PetscSFSetGraph(gsfPoint, numRoots + numGhostCells, numLeaves, glocalPoints, PETSC_OWN_POINTER, gremotePoints, PETSC_OWN_POINTER));
953: }
954: PetscCall(DMDestroy(dmSplit));
955: *dmSplit = sdm;
956: PetscFunctionReturn(PETSC_SUCCESS);
957: }
959: PetscErrorCode CreatePartitionVec(DM dm, DM *dmCell, Vec *partition)
960: {
961: PetscSF sfPoint;
962: PetscSection coordSection;
963: Vec coordinates;
964: PetscSection sectionCell;
965: PetscScalar *part;
966: PetscInt cStart, cEnd, c;
967: PetscMPIInt rank;
969: PetscFunctionBeginUser;
970: PetscCall(DMGetCoordinateSection(dm, &coordSection));
971: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
972: PetscCall(DMClone(dm, dmCell));
973: PetscCall(DMGetPointSF(dm, &sfPoint));
974: PetscCall(DMSetPointSF(*dmCell, sfPoint));
975: PetscCall(DMSetCoordinateSection(*dmCell, PETSC_DETERMINE, coordSection));
976: PetscCall(DMSetCoordinatesLocal(*dmCell, coordinates));
977: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)dm), &rank));
978: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionCell));
979: PetscCall(DMPlexGetHeightStratum(*dmCell, 0, &cStart, &cEnd));
980: PetscCall(PetscSectionSetChart(sectionCell, cStart, cEnd));
981: for (c = cStart; c < cEnd; ++c) PetscCall(PetscSectionSetDof(sectionCell, c, 1));
982: PetscCall(PetscSectionSetUp(sectionCell));
983: PetscCall(DMSetLocalSection(*dmCell, sectionCell));
984: PetscCall(PetscSectionDestroy(§ionCell));
985: PetscCall(DMCreateLocalVector(*dmCell, partition));
986: PetscCall(PetscObjectSetName((PetscObject)*partition, "partition"));
987: PetscCall(VecGetArray(*partition, &part));
988: for (c = cStart; c < cEnd; ++c) {
989: PetscScalar *p;
991: PetscCall(DMPlexPointLocalRef(*dmCell, c, part, &p));
992: p[0] = rank;
993: }
994: PetscCall(VecRestoreArray(*partition, &part));
995: PetscFunctionReturn(PETSC_SUCCESS);
996: }
998: PetscErrorCode CreateMassMatrix(DM dm, Vec *massMatrix, User user)
999: {
1000: DM plex, dmMass, dmFace, dmCell, dmCoord;
1001: PetscSection coordSection;
1002: Vec coordinates, facegeom, cellgeom;
1003: PetscSection sectionMass;
1004: PetscScalar *m;
1005: const PetscScalar *fgeom, *cgeom, *coords;
1006: PetscInt vStart, vEnd, v;
1008: PetscFunctionBeginUser;
1009: PetscCall(DMConvert(dm, DMPLEX, &plex));
1010: PetscCall(DMGetCoordinateSection(dm, &coordSection));
1011: PetscCall(DMGetCoordinatesLocal(dm, &coordinates));
1012: PetscCall(DMClone(dm, &dmMass));
1013: PetscCall(DMSetCoordinateSection(dmMass, PETSC_DETERMINE, coordSection));
1014: PetscCall(DMSetCoordinatesLocal(dmMass, coordinates));
1015: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), §ionMass));
1016: PetscCall(DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd));
1017: PetscCall(PetscSectionSetChart(sectionMass, vStart, vEnd));
1018: for (v = vStart; v < vEnd; ++v) {
1019: PetscInt numFaces;
1021: PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
1022: PetscCall(PetscSectionSetDof(sectionMass, v, numFaces * numFaces));
1023: }
1024: PetscCall(PetscSectionSetUp(sectionMass));
1025: PetscCall(DMSetLocalSection(dmMass, sectionMass));
1026: PetscCall(PetscSectionDestroy(§ionMass));
1027: PetscCall(DMGetLocalVector(dmMass, massMatrix));
1028: PetscCall(VecGetArray(*massMatrix, &m));
1029: PetscCall(DMPlexGetGeometryFVM(plex, &facegeom, &cellgeom, NULL));
1030: PetscCall(VecGetDM(facegeom, &dmFace));
1031: PetscCall(VecGetArrayRead(facegeom, &fgeom));
1032: PetscCall(VecGetDM(cellgeom, &dmCell));
1033: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
1034: PetscCall(DMGetCoordinateDM(dm, &dmCoord));
1035: PetscCall(VecGetArrayRead(coordinates, &coords));
1036: for (v = vStart; v < vEnd; ++v) {
1037: const PetscInt *faces;
1038: const PetscFVFaceGeom *fgA, *fgB, *cg;
1039: const PetscScalar *vertex;
1040: PetscInt numFaces, sides[2], f, g;
1042: PetscCall(DMPlexPointLocalRead(dmCoord, v, coords, &vertex));
1043: PetscCall(DMPlexGetSupportSize(dmMass, v, &numFaces));
1044: PetscCall(DMPlexGetSupport(dmMass, v, &faces));
1045: for (f = 0; f < numFaces; ++f) {
1046: sides[0] = faces[f];
1047: PetscCall(DMPlexPointLocalRead(dmFace, faces[f], fgeom, &fgA));
1048: for (g = 0; g < numFaces; ++g) {
1049: const PetscInt *cells = NULL;
1050: PetscReal area = 0.0;
1051: PetscInt numCells;
1053: sides[1] = faces[g];
1054: PetscCall(DMPlexPointLocalRead(dmFace, faces[g], fgeom, &fgB));
1055: PetscCall(DMPlexGetJoin(dmMass, 2, sides, &numCells, &cells));
1056: PetscCheck(numCells == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Invalid join for faces");
1057: PetscCall(DMPlexPointLocalRead(dmCell, cells[0], cgeom, &cg));
1058: area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgA->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgA->centroid[0] - cg->centroid[0]));
1059: area += PetscAbsScalar((vertex[0] - cg->centroid[0]) * (fgB->centroid[1] - cg->centroid[1]) - (vertex[1] - cg->centroid[1]) * (fgB->centroid[0] - cg->centroid[0]));
1060: m[f * numFaces + g] = Dot2(fgA->normal, fgB->normal) * area * 0.5;
1061: PetscCall(DMPlexRestoreJoin(dmMass, 2, sides, &numCells, &cells));
1062: }
1063: }
1064: }
1065: PetscCall(VecRestoreArrayRead(facegeom, &fgeom));
1066: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
1067: PetscCall(VecRestoreArrayRead(coordinates, &coords));
1068: PetscCall(VecRestoreArray(*massMatrix, &m));
1069: PetscCall(DMDestroy(&dmMass));
1070: PetscCall(DMDestroy(&plex));
1071: PetscFunctionReturn(PETSC_SUCCESS);
1072: }
1074: PetscErrorCode SetUpLocalSpace(DM dm, User user)
1075: {
1076: PetscSection stateSection;
1077: Physics phys;
1078: PetscInt dof = user->model->physics->dof, *cind, d, stateSize, cStart, cEnd, cEndInterior, c, i;
1080: PetscFunctionBeginUser;
1081: PetscCall(DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd));
1082: PetscCall(DMPlexGetGhostCellStratum(dm, &cEndInterior, NULL));
1083: PetscCall(PetscSectionCreate(PetscObjectComm((PetscObject)dm), &stateSection));
1084: phys = user->model->physics;
1085: PetscCall(PetscSectionSetNumFields(stateSection, phys->nfields));
1086: for (i = 0; i < phys->nfields; i++) {
1087: PetscCall(PetscSectionSetFieldName(stateSection, i, phys->field_desc[i].name));
1088: PetscCall(PetscSectionSetFieldComponents(stateSection, i, phys->field_desc[i].dof));
1089: }
1090: PetscCall(PetscSectionSetChart(stateSection, cStart, cEnd));
1091: for (c = cStart; c < cEnd; ++c) {
1092: for (i = 0; i < phys->nfields; i++) PetscCall(PetscSectionSetFieldDof(stateSection, c, i, phys->field_desc[i].dof));
1093: PetscCall(PetscSectionSetDof(stateSection, c, dof));
1094: }
1095: for (c = cEndInterior; c < cEnd; ++c) PetscCall(PetscSectionSetConstraintDof(stateSection, c, dof));
1096: PetscCall(PetscSectionSetUp(stateSection));
1097: PetscCall(PetscMalloc1(dof, &cind));
1098: for (d = 0; d < dof; ++d) cind[d] = d;
1099: #if 0
1100: for (c = cStart; c < cEnd; ++c) {
1101: PetscInt val;
1103: PetscCall(DMGetLabelValue(dm, "vtk", c, &val));
1104: if (val < 0) PetscCall(PetscSectionSetConstraintIndices(stateSection, c, cind));
1105: }
1106: #endif
1107: for (c = cEndInterior; c < cEnd; ++c) PetscCall(PetscSectionSetConstraintIndices(stateSection, c, cind));
1108: PetscCall(PetscFree(cind));
1109: PetscCall(PetscSectionGetStorageSize(stateSection, &stateSize));
1110: PetscCall(DMSetLocalSection(dm, stateSection));
1111: PetscCall(PetscSectionDestroy(&stateSection));
1112: PetscFunctionReturn(PETSC_SUCCESS);
1113: }
1115: #if 0
1116: PetscErrorCode SetUpBoundaries(DM dm, User user)
1117: {
1118: Model mod = user->model;
1119: BoundaryLink b;
1121: PetscFunctionBeginUser;
1122: PetscOptionsBegin(PetscObjectComm((PetscObject)dm),NULL,"Boundary condition options","");
1123: for (b = mod->boundary; b; b=b->next) {
1124: char optname[512];
1125: PetscInt ids[512],len = 512;
1126: PetscBool flg;
1127: PetscCall(PetscSNPrintf(optname,sizeof optname,"-bc_%s",b->name));
1128: PetscCall(PetscMemzero(ids,sizeof(ids)));
1129: PetscCall(PetscOptionsIntArray(optname,"List of boundary IDs","",ids,&len,&flg));
1130: if (flg) {
1131: /* TODO: check all IDs to make sure they exist in the mesh */
1132: PetscCall(PetscFree(b->ids));
1133: b->numids = len;
1134: PetscCall(PetscMalloc1(len,&b->ids));
1135: PetscCall(PetscArraycpy(b->ids,ids,len));
1136: }
1137: }
1138: PetscOptionsEnd();
1139: PetscFunctionReturn(PETSC_SUCCESS);
1140: }
1141: #endif
1143: /* Behavior will be different for multi-physics or when using non-default boundary conditions */
1144: static PetscErrorCode ModelSolutionSetDefault(Model mod, SolutionFunction func, void *ctx)
1145: {
1146: PetscFunctionBeginUser;
1147: mod->solution = func;
1148: mod->solutionctx = ctx;
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: static PetscErrorCode ModelFunctionalRegister(Model mod, const char *name, PetscInt *offset, FunctionalFunction func, void *ctx)
1153: {
1154: FunctionalLink link, *ptr;
1155: PetscInt lastoffset = -1;
1157: PetscFunctionBeginUser;
1158: for (ptr = &mod->functionalRegistry; *ptr; ptr = &(*ptr)->next) lastoffset = (*ptr)->offset;
1159: PetscCall(PetscNew(&link));
1160: PetscCall(PetscStrallocpy(name, &link->name));
1161: link->offset = lastoffset + 1;
1162: link->func = func;
1163: link->ctx = ctx;
1164: link->next = NULL;
1165: *ptr = link;
1166: *offset = link->offset;
1167: PetscFunctionReturn(PETSC_SUCCESS);
1168: }
1170: static PetscErrorCode ModelFunctionalSetFromOptions(Model mod, PetscOptions *PetscOptionsObject)
1171: {
1172: PetscInt i, j;
1173: FunctionalLink link;
1174: char *names[256];
1176: PetscFunctionBeginUser;
1177: mod->numMonitored = PETSC_STATIC_ARRAY_LENGTH(names);
1178: PetscCall(PetscOptionsStringArray("-monitor", "list of functionals to monitor", "", names, &mod->numMonitored, NULL));
1179: /* Create list of functionals that will be computed somehow */
1180: PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalMonitored));
1181: /* Create index of calls that we will have to make to compute these functionals (over-allocation in general). */
1182: PetscCall(PetscMalloc1(mod->numMonitored, &mod->functionalCall));
1183: mod->numCall = 0;
1184: for (i = 0; i < mod->numMonitored; i++) {
1185: for (link = mod->functionalRegistry; link; link = link->next) {
1186: PetscBool match;
1187: PetscCall(PetscStrcasecmp(names[i], link->name, &match));
1188: if (match) break;
1189: }
1190: PetscCheck(link, mod->comm, PETSC_ERR_USER, "No known functional '%s'", names[i]);
1191: mod->functionalMonitored[i] = link;
1192: for (j = 0; j < i; j++) {
1193: if (mod->functionalCall[j]->func == link->func && mod->functionalCall[j]->ctx == link->ctx) goto next_name;
1194: }
1195: mod->functionalCall[mod->numCall++] = link; /* Just points to the first link using the result. There may be more results. */
1196: next_name:
1197: PetscCall(PetscFree(names[i]));
1198: }
1200: /* Find out the maximum index of any functional computed by a function we will be calling (even if we are not using it) */
1201: mod->maxComputed = -1;
1202: for (link = mod->functionalRegistry; link; link = link->next) {
1203: for (i = 0; i < mod->numCall; i++) {
1204: FunctionalLink call = mod->functionalCall[i];
1205: if (link->func == call->func && link->ctx == call->ctx) mod->maxComputed = PetscMax(mod->maxComputed, link->offset);
1206: }
1207: }
1208: PetscFunctionReturn(PETSC_SUCCESS);
1209: }
1211: static PetscErrorCode FunctionalLinkDestroy(FunctionalLink *link)
1212: {
1213: FunctionalLink l, next;
1215: PetscFunctionBeginUser;
1216: if (!link) PetscFunctionReturn(PETSC_SUCCESS);
1217: l = *link;
1218: *link = NULL;
1219: for (; l; l = next) {
1220: next = l->next;
1221: PetscCall(PetscFree(l->name));
1222: PetscCall(PetscFree(l));
1223: }
1224: PetscFunctionReturn(PETSC_SUCCESS);
1225: }
1227: PetscErrorCode SetInitialCondition(DM dm, Vec X, User user)
1228: {
1229: DM plex, dmCell;
1230: Model mod = user->model;
1231: Vec cellgeom;
1232: const PetscScalar *cgeom;
1233: PetscScalar *x;
1234: PetscInt cStart, cEnd, c;
1236: PetscFunctionBeginUser;
1237: PetscCall(DMConvert(dm, DMPLEX, &plex));
1238: PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1239: PetscCall(VecGetDM(cellgeom, &dmCell));
1240: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1241: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
1242: PetscCall(VecGetArray(X, &x));
1243: for (c = cStart; c < cEnd; ++c) {
1244: const PetscFVCellGeom *cg;
1245: PetscScalar *xc;
1247: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
1248: PetscCall(DMPlexPointGlobalRef(dm, c, x, &xc));
1249: if (xc) PetscCall((*mod->solution)(mod, 0.0, cg->centroid, xc, mod->solutionctx));
1250: }
1251: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
1252: PetscCall(VecRestoreArray(X, &x));
1253: PetscCall(DMDestroy(&plex));
1254: PetscFunctionReturn(PETSC_SUCCESS);
1255: }
1257: static PetscErrorCode OutputVTK(DM dm, const char *filename, PetscViewer *viewer)
1258: {
1259: PetscFunctionBeginUser;
1260: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1261: PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERVTK));
1262: PetscCall(PetscViewerFileSetName(*viewer, filename));
1263: PetscFunctionReturn(PETSC_SUCCESS);
1264: }
1266: static PetscErrorCode MonitorVTK(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
1267: {
1268: User user = (User)ctx;
1269: DM dm, plex;
1270: PetscViewer viewer;
1271: char filename[PETSC_MAX_PATH_LEN], *ftable = NULL;
1272: PetscReal xnorm;
1274: PetscFunctionBeginUser;
1275: PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1276: PetscCall(VecGetDM(X, &dm));
1277: PetscCall(VecNorm(X, NORM_INFINITY, &xnorm));
1278: if (stepnum >= 0) { /* No summary for final time */
1279: Model mod = user->model;
1280: Vec cellgeom;
1281: PetscInt c, cStart, cEnd, fcount, i;
1282: size_t ftableused, ftablealloc;
1283: const PetscScalar *cgeom, *x;
1284: DM dmCell;
1285: PetscReal *fmin, *fmax, *fintegral, *ftmp;
1287: PetscCall(DMConvert(dm, DMPLEX, &plex));
1288: PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1289: fcount = mod->maxComputed + 1;
1290: PetscCall(PetscMalloc4(fcount, &fmin, fcount, &fmax, fcount, &fintegral, fcount, &ftmp));
1291: for (i = 0; i < fcount; i++) {
1292: fmin[i] = PETSC_MAX_REAL;
1293: fmax[i] = PETSC_MIN_REAL;
1294: fintegral[i] = 0;
1295: }
1296: PetscCall(DMPlexGetSimplexOrBoxCells(dm, 0, &cStart, &cEnd));
1297: PetscCall(VecGetDM(cellgeom, &dmCell));
1298: PetscCall(VecGetArrayRead(cellgeom, &cgeom));
1299: PetscCall(VecGetArrayRead(X, &x));
1300: for (c = cStart; c < cEnd; ++c) {
1301: const PetscFVCellGeom *cg;
1302: const PetscScalar *cx;
1303: PetscCall(DMPlexPointLocalRead(dmCell, c, cgeom, &cg));
1304: PetscCall(DMPlexPointGlobalRead(dm, c, x, &cx));
1305: if (!cx) continue; /* not a global cell */
1306: for (i = 0; i < mod->numCall; i++) {
1307: FunctionalLink flink = mod->functionalCall[i];
1308: PetscCall((*flink->func)(mod, time, cg->centroid, cx, ftmp, flink->ctx));
1309: }
1310: for (i = 0; i < fcount; i++) {
1311: fmin[i] = PetscMin(fmin[i], ftmp[i]);
1312: fmax[i] = PetscMax(fmax[i], ftmp[i]);
1313: fintegral[i] += cg->volume * ftmp[i];
1314: }
1315: }
1316: PetscCall(VecRestoreArrayRead(cellgeom, &cgeom));
1317: PetscCall(VecRestoreArrayRead(X, &x));
1318: PetscCall(DMDestroy(&plex));
1319: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmin, fcount, MPIU_REAL, MPIU_MIN, PetscObjectComm((PetscObject)ts)));
1320: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fmax, fcount, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)ts)));
1321: PetscCall(MPIU_Allreduce(MPI_IN_PLACE, fintegral, fcount, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)ts)));
1323: ftablealloc = fcount * 100;
1324: ftableused = 0;
1325: PetscCall(PetscMalloc1(ftablealloc, &ftable));
1326: for (i = 0; i < mod->numMonitored; i++) {
1327: size_t countused;
1328: char buffer[256], *p;
1329: FunctionalLink flink = mod->functionalMonitored[i];
1330: PetscInt id = flink->offset;
1331: if (i % 3) {
1332: PetscCall(PetscArraycpy(buffer, " ", 2));
1333: p = buffer + 2;
1334: } else if (i) {
1335: char newline[] = "\n";
1336: PetscCall(PetscArraycpy(buffer, newline, sizeof(newline) - 1));
1337: p = buffer + sizeof(newline) - 1;
1338: } else {
1339: p = buffer;
1340: }
1341: 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]));
1342: countused += p - buffer;
1343: if (countused > ftablealloc - ftableused - 1) { /* reallocate */
1344: char *ftablenew;
1345: ftablealloc = 2 * ftablealloc + countused;
1346: PetscCall(PetscMalloc(ftablealloc, &ftablenew));
1347: PetscCall(PetscArraycpy(ftablenew, ftable, ftableused));
1348: PetscCall(PetscFree(ftable));
1349: ftable = ftablenew;
1350: }
1351: PetscCall(PetscArraycpy(ftable + ftableused, buffer, countused));
1352: ftableused += countused;
1353: ftable[ftableused] = 0;
1354: }
1355: PetscCall(PetscFree4(fmin, fmax, fintegral, ftmp));
1357: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "% 3" PetscInt_FMT " time %8.4g |x| %8.4g %s\n", stepnum, (double)time, (double)xnorm, ftable ? ftable : ""));
1358: PetscCall(PetscFree(ftable));
1359: }
1360: if (user->vtkInterval < 1) PetscFunctionReturn(PETSC_SUCCESS);
1361: if ((stepnum == -1) ^ (stepnum % user->vtkInterval == 0)) {
1362: if (stepnum == -1) { /* Final time is not multiple of normal time interval, write it anyway */
1363: PetscCall(TSGetStepNumber(ts, &stepnum));
1364: }
1365: PetscCall(PetscSNPrintf(filename, sizeof filename, "ex11-%03" PetscInt_FMT ".vtu", stepnum));
1366: PetscCall(OutputVTK(dm, filename, &viewer));
1367: PetscCall(VecView(X, viewer));
1368: PetscCall(PetscViewerDestroy(&viewer));
1369: }
1370: PetscFunctionReturn(PETSC_SUCCESS);
1371: }
1373: static PetscErrorCode OutputBIN(DM dm, const char *filename, PetscViewer *viewer)
1374: {
1375: PetscFunctionBeginUser;
1376: PetscCall(PetscViewerCreate(PetscObjectComm((PetscObject)dm), viewer));
1377: PetscCall(PetscViewerSetType(*viewer, PETSCVIEWERBINARY));
1378: PetscCall(PetscViewerFileSetMode(*viewer, FILE_MODE_WRITE));
1379: PetscCall(PetscViewerFileSetName(*viewer, filename));
1380: PetscFunctionReturn(PETSC_SUCCESS);
1381: }
1383: static PetscErrorCode TestMonitor(DM dm, const char *filename, Vec X, PetscReal time)
1384: {
1385: Vec odesolution;
1386: PetscInt Nr;
1387: PetscBool equal;
1388: PetscReal timeread;
1389: PetscViewer viewer;
1391: PetscFunctionBeginUser;
1392: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &viewer));
1393: PetscCall(VecCreate(PETSC_COMM_WORLD, &odesolution));
1394: PetscCall(VecLoad(odesolution, viewer));
1395: VecEqual(X, odesolution, &equal);
1396: PetscCheck(equal, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the vec data from file");
1397: else
1398: {
1399: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for Vec\n"));
1400: }
1401: /*Nr = 1;
1402: PetscCall(PetscRealLoad(Nr,&Nr,&timeread,viewer));*/
1403: PetscCall(PetscViewerBinaryRead(viewer, &timeread, 1, NULL, PETSC_REAL));
1405: PetscCheck(timeread == time, PETSC_COMM_WORLD, PETSC_ERR_FILE_UNEXPECTED, "Error in reading the scalar data from file");
1406: else
1407: {
1408: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "IO test OK for PetscReal\n"));
1409: }
1411: PetscCall(PetscViewerDestroy(&viewer));
1412: PetscCall(VecDestroy(&odesolution));
1413: PetscFunctionReturn(PETSC_SUCCESS);
1414: }
1416: static PetscErrorCode MonitorBIN(TS ts, PetscInt stepnum, PetscReal time, Vec X, void *ctx)
1417: {
1418: User user = (User)ctx;
1419: DM dm;
1420: PetscViewer viewer;
1421: char filename[PETSC_MAX_PATH_LEN];
1423: PetscFunctionBeginUser;
1424: PetscCall(VecGetDM(X, &dm));
1425: PetscCall(PetscSNPrintf(filename, sizeof filename, "ex11-SA-%06d.bin", stepnum));
1426: PetscCall(OutputBIN(dm, filename, &viewer));
1427: PetscCall(VecView(X, viewer));
1428: PetscCall(PetscRealView(1, &time, viewer));
1429: /* PetscCall(PetscViewerBinaryWrite(viewer,&time,1,PETSC_REAL));*/
1430: PetscCall(PetscViewerDestroy(&viewer));
1431: PetscCall(TestMonitor(dm, filename, X, time));
1432: PetscFunctionReturn(PETSC_SUCCESS);
1433: }
1435: int main(int argc, char **argv)
1436: {
1437: MPI_Comm comm;
1438: PetscFV fvm;
1439: User user;
1440: Model mod;
1441: Physics phys;
1442: DM dm, plex;
1443: PetscReal ftime, cfl, dt, minRadius;
1444: PetscInt dim, nsteps;
1445: TS ts;
1446: TSConvergedReason reason;
1447: Vec X;
1448: PetscViewer viewer;
1449: PetscMPIInt rank;
1450: PetscBool vtkCellGeom, splitFaces;
1451: PetscInt overlap, f;
1452: char filename[PETSC_MAX_PATH_LEN] = "sevenside.exo";
1454: PetscFunctionBeginUser;
1455: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1456: comm = PETSC_COMM_WORLD;
1457: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1459: PetscCall(PetscNew(&user));
1460: PetscCall(PetscNew(&user->model));
1461: PetscCall(PetscNew(&user->model->physics));
1462: mod = user->model;
1463: phys = mod->physics;
1464: mod->comm = comm;
1466: /* Register physical models to be available on the command line */
1467: PetscCall(PetscFunctionListAdd(&PhysicsList, "advect", PhysicsCreate_Advect));
1468: PetscCall(PetscFunctionListAdd(&PhysicsList, "sw", PhysicsCreate_SW));
1469: PetscCall(PetscFunctionListAdd(&PhysicsList, "euler", PhysicsCreate_Euler));
1471: PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Mesh Options", "");
1472: {
1473: cfl = 0.9 * 4; /* default SSPRKS2 with s=5 stages is stable for CFL number s-1 */
1474: PetscCall(PetscOptionsReal("-ufv_cfl", "CFL number per step", "", cfl, &cfl, NULL));
1475: PetscCall(PetscOptionsString("-f", "Exodus.II filename to read", "", filename, filename, sizeof(filename), NULL));
1476: splitFaces = PETSC_FALSE;
1477: PetscCall(PetscOptionsBool("-ufv_split_faces", "Split faces between cell sets", "", splitFaces, &splitFaces, NULL));
1478: overlap = 1;
1479: PetscCall(PetscOptionsInt("-ufv_mesh_overlap", "Number of cells to overlap partitions", "", overlap, &overlap, NULL));
1480: user->vtkInterval = 1;
1481: PetscCall(PetscOptionsInt("-ufv_vtk_interval", "VTK output interval (0 to disable)", "", user->vtkInterval, &user->vtkInterval, NULL));
1482: vtkCellGeom = PETSC_FALSE;
1483: PetscCall(PetscOptionsBool("-ufv_vtk_cellgeom", "Write cell geometry (for debugging)", "", vtkCellGeom, &vtkCellGeom, NULL));
1484: }
1485: PetscOptionsEnd();
1486: PetscCall(DMPlexCreateExodusFromFile(comm, filename, PETSC_TRUE, &dm));
1487: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1488: PetscCall(DMGetDimension(dm, &dim));
1490: PetscOptionsBegin(comm, NULL, "Unstructured Finite Volume Physics Options", "");
1491: {
1492: PetscDS prob;
1493: PetscErrorCode (*physcreate)(PetscDS, Model, Physics);
1494: char physname[256] = "advect";
1496: PetscCall(DMCreateLabel(dm, "Face Sets"));
1497: PetscCall(PetscOptionsFList("-physics", "Physics module to solve", "", PhysicsList, physname, physname, sizeof physname, NULL));
1498: PetscCall(PetscFunctionListFind(PhysicsList, physname, &physcreate));
1499: PetscCall(PetscMemzero(phys, sizeof(struct _n_Physics)));
1500: PetscCall(DMGetDS(dm, &prob));
1501: PetscCall((*physcreate)(prob, mod, phys));
1502: mod->maxspeed = phys->maxspeed;
1503: /* Count number of fields and dofs */
1504: for (phys->nfields = 0, phys->dof = 0; phys->field_desc[phys->nfields].name; phys->nfields++) phys->dof += phys->field_desc[phys->nfields].dof;
1506: PetscCheck(mod->maxspeed > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set maxspeed", physname);
1507: PetscCheck(phys->dof > 0, comm, PETSC_ERR_ARG_WRONGSTATE, "Physics '%s' did not set dof", physname);
1508: PetscCall(ModelFunctionalSetFromOptions(mod, PetscOptionsObject));
1509: }
1510: PetscOptionsEnd();
1511: {
1512: DM dmDist;
1514: PetscCall(DMSetBasicAdjacency(dm, PETSC_TRUE, PETSC_FALSE));
1515: PetscCall(DMPlexDistribute(dm, overlap, NULL, &dmDist));
1516: if (dmDist) {
1517: PetscCall(DMDestroy(&dm));
1518: dm = dmDist;
1519: }
1520: }
1521: PetscCall(DMSetFromOptions(dm));
1522: {
1523: DM gdm;
1525: PetscCall(DMPlexConstructGhostCells(dm, NULL, NULL, &gdm));
1526: PetscCall(DMDestroy(&dm));
1527: dm = gdm;
1528: PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));
1529: }
1530: if (splitFaces) PetscCall(ConstructCellBoundary(dm, user));
1531: PetscCall(SplitFaces(&dm, "split faces", user));
1533: PetscCall(PetscFVCreate(comm, &fvm));
1534: PetscCall(PetscFVSetFromOptions(fvm));
1535: PetscCall(DMSetNumFields(dm, phys->nfields));
1536: for (f = 0; f < phys->nfields; ++f) PetscCall(DMSetField(dm, f, (PetscObject)fvm));
1537: PetscCall(PetscFVSetNumComponents(fvm, phys->dof));
1538: PetscCall(PetscFVSetSpatialDimension(fvm, dim));
1540: /* Set up DM with section describing local vector and configure local vector. */
1541: PetscCall(SetUpLocalSpace(dm, user));
1543: PetscCall(TSCreate(comm, &ts));
1544: PetscCall(TSSetType(ts, TSSSP));
1545: /* PetscCall(TSSetType(ts, TSRK)); */
1546: PetscCall(TSSetDM(ts, dm));
1547: /* PetscCall(TSMonitorSet(ts,MonitorVTK,user,NULL)); */
1548: PetscCall(TSMonitorSet(ts, MonitorBIN, user, NULL));
1549: PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, user));
1550: PetscCall(DMCreateGlobalVector(dm, &X));
1551: PetscCall(PetscObjectSetName((PetscObject)X, "solution"));
1552: PetscCall(SetInitialCondition(dm, X, user));
1553: PetscCall(DMConvert(dm, DMPLEX, &plex));
1554: if (vtkCellGeom) {
1555: DM dmCell;
1556: Vec cellgeom, partition;
1558: PetscCall(DMPlexGetGeometryFVM(plex, NULL, &cellgeom, NULL));
1559: PetscCall(OutputVTK(dm, "ex11-cellgeom.vtk", &viewer));
1560: PetscCall(VecView(cellgeom, viewer));
1561: PetscCall(PetscViewerDestroy(&viewer));
1562: PetscCall(CreatePartitionVec(dm, &dmCell, &partition));
1563: PetscCall(OutputVTK(dmCell, "ex11-partition.vtk", &viewer));
1564: PetscCall(VecView(partition, viewer));
1565: PetscCall(PetscViewerDestroy(&viewer));
1566: PetscCall(VecDestroy(&partition));
1567: PetscCall(DMDestroy(&dmCell));
1568: }
1570: PetscCall(DMPlexGetGeometryFVM(plex, NULL, NULL, &minRadius));
1571: PetscCall(DMDestroy(&plex));
1572: PetscCall(TSSetMaxTime(ts, 2.0));
1573: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1574: dt = cfl * minRadius / user->model->maxspeed;
1575: PetscCall(TSSetTimeStep(ts, dt));
1576: PetscCall(TSSetFromOptions(ts));
1577: PetscCall(TSSolve(ts, X));
1578: PetscCall(TSGetSolveTime(ts, &ftime));
1579: PetscCall(TSGetStepNumber(ts, &nsteps));
1580: PetscCall(TSGetConvergedReason(ts, &reason));
1581: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%s at time %g after %" PetscInt_FMT " steps\n", TSConvergedReasons[reason], (double)ftime, nsteps));
1582: PetscCall(TSDestroy(&ts));
1584: PetscCall(PetscFunctionListDestroy(&PhysicsList));
1585: PetscCall(FunctionalLinkDestroy(&user->model->functionalRegistry));
1586: PetscCall(PetscFree(user->model->functionalMonitored));
1587: PetscCall(PetscFree(user->model->functionalCall));
1588: PetscCall(PetscFree(user->model->physics->data));
1589: PetscCall(PetscFree(user->model->physics));
1590: PetscCall(PetscFree(user->model));
1591: PetscCall(PetscFree(user));
1592: PetscCall(VecDestroy(&X));
1593: PetscCall(PetscFVDestroy(&fvm));
1594: PetscCall(DMDestroy(&dm));
1595: PetscCall(PetscFinalize());
1596: return 0;
1597: }