Actual source code: ex9.c
1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
2: "Solves scalar and vector problems, choose the physical model with -physics\n"
3: " advection - Constant coefficient scalar advection\n"
4: " u_t + (a*u)_x = 0\n"
5: " burgers - Burgers equation\n"
6: " u_t + (u^2/2)_x = 0\n"
7: " traffic - Traffic equation\n"
8: " u_t + (u*(1-u))_x = 0\n"
9: " acoustics - Acoustic wave propagation\n"
10: " u_t + (c*z*v)_x = 0\n"
11: " v_t + (c/z*u)_x = 0\n"
12: " isogas - Isothermal gas dynamics\n"
13: " rho_t + (rho*u)_x = 0\n"
14: " (rho*u)_t + (rho*u^2 + c^2*rho)_x = 0\n"
15: " shallow - Shallow water equations\n"
16: " h_t + (h*u)_x = 0\n"
17: " (h*u)_t + (h*u^2 + g*h^2/2)_x = 0\n"
18: "Some of these physical models have multiple Riemann solvers, select these with -physics_xxx_riemann\n"
19: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
20: " the states across shocks and rarefactions\n"
21: " roe - Linearized scheme, usually with an entropy fix inside sonic rarefactions\n"
22: "The systems provide a choice of reconstructions with -physics_xxx_reconstruct\n"
23: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
24: " conservative - Limit the conservative variables directly, can cause undesired interaction of waves\n\n"
25: "A variety of limiters for high-resolution TVD limiters are available with -limit\n"
26: " upwind,minmod,superbee,mc,vanleer,vanalbada,koren,cada-torillhon (last two are nominally third order)\n"
27: " and non-TVD schemes lax-wendroff,beam-warming,fromm\n\n"
28: "To preserve the TVD property, one should time step with a strong stability preserving method.\n"
29: "The optimal high order explicit Runge-Kutta methods in TSSSP are recommended for non-stiff problems.\n\n"
30: "Several initial conditions can be chosen with -initial N\n\n"
31: "The problem size should be set with -da_grid_x M\n\n";
33: #include <petscts.h>
34: #include <petscdm.h>
35: #include <petscdmda.h>
37: #include <petsc/private/kernels/blockinvert.h>
39: static inline PetscReal Sgn(PetscReal a)
40: {
41: return (a < 0) ? -1 : 1;
42: }
43: static inline PetscReal Abs(PetscReal a)
44: {
45: return (a < 0) ? 0 : a;
46: }
47: static inline PetscReal Sqr(PetscReal a)
48: {
49: return a * a;
50: }
51: static inline PetscReal MaxAbs(PetscReal a, PetscReal b)
52: {
53: return (PetscAbs(a) > PetscAbs(b)) ? a : b;
54: }
55: PETSC_UNUSED static inline PetscReal MinAbs(PetscReal a, PetscReal b)
56: {
57: return (PetscAbs(a) < PetscAbs(b)) ? a : b;
58: }
59: static inline PetscReal MinMod2(PetscReal a, PetscReal b)
60: {
61: return (a * b < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscAbs(b));
62: }
63: static inline PetscReal MaxMod2(PetscReal a, PetscReal b)
64: {
65: return (a * b < 0) ? 0 : Sgn(a) * PetscMax(PetscAbs(a), PetscAbs(b));
66: }
67: static inline PetscReal MinMod3(PetscReal a, PetscReal b, PetscReal c)
68: {
69: return (a * b < 0 || a * c < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscMin(PetscAbs(b), PetscAbs(c)));
70: }
72: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
73: {
74: PetscReal range = xmax - xmin;
75: return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
76: }
78: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
79: typedef struct _LimitInfo {
80: PetscReal hx;
81: PetscInt m;
82: } *LimitInfo;
83: static void Limit_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
84: {
85: PetscInt i;
86: for (i = 0; i < info->m; i++) lmt[i] = 0;
87: }
88: static void Limit_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
89: {
90: PetscInt i;
91: for (i = 0; i < info->m; i++) lmt[i] = jR[i];
92: }
93: static void Limit_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
94: {
95: PetscInt i;
96: for (i = 0; i < info->m; i++) lmt[i] = jL[i];
97: }
98: static void Limit_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
99: {
100: PetscInt i;
101: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]);
102: }
103: static void Limit_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
104: {
105: PetscInt i;
106: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]);
107: }
108: static void Limit_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
109: {
110: PetscInt i;
111: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i]));
112: }
113: static void Limit_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
114: {
115: PetscInt i;
116: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]);
117: }
118: static void Limit_VanLeer(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
119: { /* phi = (t + abs(t)) / (1 + abs(t)) */
120: PetscInt i;
121: for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * Abs(jR[i]) + Abs(jL[i]) * jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15);
122: }
123: static void Limit_VanAlbada(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
124: { /* phi = (t + t^2) / (1 + t^2) */
125: PetscInt i;
126: for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
127: }
128: static void Limit_VanAlbadaTVD(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
129: { /* phi = (t + t^2) / (1 + t^2) */
130: PetscInt i;
131: for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * jR[i] < 0) ? 0 : (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
132: }
133: static void Limit_Koren(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
134: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
135: PetscInt i;
136: for (i = 0; i < info->m; i++) lmt[i] = ((jL[i] * Sqr(jR[i]) + 2 * Sqr(jL[i]) * jR[i]) / (2 * Sqr(jL[i]) - jL[i] * jR[i] + 2 * Sqr(jR[i]) + 1e-15));
137: }
138: static void Limit_KorenSym(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
139: { /* Symmetric version of above */
140: PetscInt i;
141: for (i = 0; i < info->m; i++) lmt[i] = (1.5 * (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (2 * Sqr(jL[i]) - jL[i] * jR[i] + 2 * Sqr(jR[i]) + 1e-15));
142: }
143: static void Limit_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
144: { /* Eq 11 of Cada-Torrilhon 2009 */
145: PetscInt i;
146: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]);
147: }
148: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L, PetscReal R)
149: {
150: return PetscMax(0, PetscMin((L + 2 * R) / 3, PetscMax(-0.5 * L, PetscMin(2 * L, PetscMin((L + 2 * R) / 3, 1.6 * R)))));
151: }
152: static void Limit_CadaTorrilhon2(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
153: { /* Cada-Torrilhon 2009, Eq 13 */
154: PetscInt i;
155: for (i = 0; i < info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]);
156: }
157: static void Limit_CadaTorrilhon3R(PetscReal r, LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
158: { /* Cada-Torrilhon 2009, Eq 22 */
159: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
160: const PetscReal eps = 1e-7, hx = info->hx;
161: for (PetscInt i = 0; i < info->m; i++) {
162: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r * hx);
163: lmt[i] = ((eta < 1 - eps) ? (jL[i] + 2 * jR[i]) / 3 : ((eta > 1 + eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]) : 0.5 * ((1 - (eta - 1) / eps) * (jL[i] + 2 * jR[i]) / 3 + (1 + (eta + 1) / eps) * CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]))));
164: }
165: }
166: static void Limit_CadaTorrilhon3R0p1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
167: {
168: Limit_CadaTorrilhon3R(0.1, info, jL, jR, lmt);
169: }
170: static void Limit_CadaTorrilhon3R1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
171: {
172: Limit_CadaTorrilhon3R(1, info, jL, jR, lmt);
173: }
174: static void Limit_CadaTorrilhon3R10(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
175: {
176: Limit_CadaTorrilhon3R(10, info, jL, jR, lmt);
177: }
178: static void Limit_CadaTorrilhon3R100(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
179: {
180: Limit_CadaTorrilhon3R(100, info, jL, jR, lmt);
181: }
183: /* --------------------------------- Finite Volume data structures ----------------------------------- */
185: typedef enum {
186: FVBC_PERIODIC,
187: FVBC_OUTFLOW
188: } FVBCType;
189: static const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "FVBCType", "FVBC_", NULL};
190: typedef PetscErrorCode (*RiemannFunction)(void *, PetscInt, const PetscScalar *, const PetscScalar *, PetscScalar *, PetscReal *);
191: typedef PetscErrorCode (*ReconstructFunction)(void *, PetscInt, const PetscScalar *, PetscScalar *, PetscScalar *, PetscReal *);
193: typedef struct {
194: PetscErrorCode (*sample)(void *, PetscInt, FVBCType, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *);
195: RiemannFunction riemann;
196: ReconstructFunction characteristic;
197: PetscErrorCode (*destroy)(void *);
198: void *user;
199: PetscInt dof;
200: char *fieldname[16];
201: } PhysicsCtx;
203: typedef struct {
204: void (*limit)(LimitInfo, const PetscScalar *, const PetscScalar *, PetscScalar *);
205: PhysicsCtx physics;
206: MPI_Comm comm;
207: char prefix[256];
209: /* Local work arrays */
210: PetscScalar *R, *Rinv; /* Characteristic basis, and it's inverse. COLUMN-MAJOR */
211: PetscScalar *cjmpLR; /* Jumps at left and right edge of cell, in characteristic basis, len=2*dof */
212: PetscScalar *cslope; /* Limited slope, written in characteristic basis */
213: PetscScalar *uLR; /* Solution at left and right of interface, conservative variables, len=2*dof */
214: PetscScalar *flux; /* Flux across interface */
215: PetscReal *speeds; /* Speeds of each wave */
217: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
218: PetscReal cfl;
219: PetscReal xmin, xmax;
220: PetscInt initial;
221: PetscBool exact;
222: FVBCType bctype;
223: } FVCtx;
225: PetscErrorCode RiemannListAdd(PetscFunctionList *flist, const char *name, RiemannFunction rsolve)
226: {
227: PetscFunctionBeginUser;
228: PetscCall(PetscFunctionListAdd(flist, name, rsolve));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: PetscErrorCode RiemannListFind(PetscFunctionList flist, const char *name, RiemannFunction *rsolve)
233: {
234: PetscFunctionBeginUser;
235: PetscCall(PetscFunctionListFind(flist, name, rsolve));
236: PetscCheck(*rsolve, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Riemann solver \"%s\" could not be found", name);
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist, const char *name, ReconstructFunction r)
241: {
242: PetscFunctionBeginUser;
243: PetscCall(PetscFunctionListAdd(flist, name, r));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: PetscErrorCode ReconstructListFind(PetscFunctionList flist, const char *name, ReconstructFunction *r)
248: {
249: PetscFunctionBeginUser;
250: PetscCall(PetscFunctionListFind(flist, name, r));
251: PetscCheck(*r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Reconstruction \"%s\" could not be found", name);
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: /* --------------------------------- Physics ----------------------------------- */
256: /*
257: Each physical model consists of Riemann solver and a function to determine the basis to use for reconstruction. These
258: are set with the PhysicsCreate_XXX function which allocates private storage and sets these methods as well as the
259: number of fields and their names, and a function to deallocate private storage.
260: */
262: /* First a few functions useful to several different physics */
263: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
264: {
265: PetscFunctionBeginUser;
266: for (PetscInt i = 0; i < m; i++) {
267: for (PetscInt j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j);
268: speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
269: }
270: PetscFunctionReturn(PETSC_SUCCESS);
271: }
273: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
274: {
275: PetscFunctionBeginUser;
276: PetscCall(PetscFree(vctx));
277: PetscFunctionReturn(PETSC_SUCCESS);
278: }
280: /* --------------------------------- Advection ----------------------------------- */
282: typedef struct {
283: PetscReal a; /* advective velocity */
284: } AdvectCtx;
286: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
287: {
288: AdvectCtx *ctx = (AdvectCtx *)vctx;
289: PetscReal speed;
291: PetscFunctionBeginUser;
292: speed = ctx->a;
293: flux[0] = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
294: *maxspeed = speed;
295: PetscFunctionReturn(PETSC_SUCCESS);
296: }
298: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
299: {
300: AdvectCtx *ctx = (AdvectCtx *)vctx;
302: PetscFunctionBeginUser;
303: X[0] = 1.;
304: Xi[0] = 1.;
305: speeds[0] = ctx->a;
306: PetscFunctionReturn(PETSC_SUCCESS);
307: }
309: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
310: {
311: AdvectCtx *ctx = (AdvectCtx *)vctx;
312: PetscReal a = ctx->a, x0;
314: PetscFunctionBeginUser;
315: switch (bctype) {
316: case FVBC_OUTFLOW:
317: x0 = x - a * t;
318: break;
319: case FVBC_PERIODIC:
320: x0 = RangeMod(x - a * t, xmin, xmax);
321: break;
322: default:
323: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
324: }
325: switch (initial) {
326: case 0:
327: u[0] = (x0 < 0) ? 1 : -1;
328: break;
329: case 1:
330: u[0] = (x0 < 0) ? -1 : 1;
331: break;
332: case 2:
333: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
334: break;
335: case 3:
336: u[0] = PetscSinReal(2 * PETSC_PI * x0);
337: break;
338: case 4:
339: u[0] = PetscAbs(x0);
340: break;
341: case 5:
342: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
343: break;
344: case 6:
345: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
346: break;
347: case 7:
348: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
349: break;
350: default:
351: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
357: {
358: AdvectCtx *user;
360: PetscFunctionBeginUser;
361: PetscCall(PetscNew(&user));
362: ctx->physics.sample = PhysicsSample_Advect;
363: ctx->physics.riemann = PhysicsRiemann_Advect;
364: ctx->physics.characteristic = PhysicsCharacteristic_Advect;
365: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
366: ctx->physics.user = user;
367: ctx->physics.dof = 1;
368: PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
369: user->a = 1;
370: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
371: {
372: PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
373: }
374: PetscOptionsEnd();
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /* --------------------------------- Burgers ----------------------------------- */
380: typedef struct {
381: PetscReal lxf_speed;
382: } BurgersCtx;
384: static PetscErrorCode PhysicsSample_Burgers(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
385: {
386: PetscFunctionBeginUser;
387: PetscCheck(bctype != FVBC_PERIODIC || t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solution not implemented for periodic");
388: switch (initial) {
389: case 0:
390: u[0] = (x < 0) ? 1 : -1;
391: break;
392: case 1:
393: if (x < -t) u[0] = -1;
394: else if (x < t) u[0] = x / t;
395: else u[0] = 1;
396: break;
397: case 2:
398: if (x <= 0) u[0] = 0;
399: else if (x < t) u[0] = x / t;
400: else if (x < 1 + 0.5 * t) u[0] = 1;
401: else u[0] = 0;
402: break;
403: case 3:
404: if (x < 0.2 * t) u[0] = 0.2;
405: else if (x < t) u[0] = x / t;
406: else u[0] = 1;
407: break;
408: case 4:
409: PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only initial condition available");
410: u[0] = 0.7 + 0.3 * PetscSinReal(2 * PETSC_PI * ((x - xmin) / (xmax - xmin)));
411: break;
412: case 5: /* Pure shock solution */
413: if (x < 0.5 * t) u[0] = 1;
414: else u[0] = 0;
415: break;
416: default:
417: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
418: }
419: PetscFunctionReturn(PETSC_SUCCESS);
420: }
422: static PetscErrorCode PhysicsRiemann_Burgers_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
423: {
424: PetscFunctionBeginUser;
425: if (uL[0] < uR[0]) { /* rarefaction */
426: flux[0] = (uL[0] * uR[0] < 0) ? 0 /* sonic rarefaction */
427: : 0.5 * PetscMin(PetscSqr(uL[0]), PetscSqr(uR[0]));
428: } else { /* shock */
429: flux[0] = 0.5 * PetscMax(PetscSqr(uL[0]), PetscSqr(uR[0]));
430: }
431: *maxspeed = (PetscAbs(uL[0]) > PetscAbs(uR[0])) ? uL[0] : uR[0];
432: PetscFunctionReturn(PETSC_SUCCESS);
433: }
435: static PetscErrorCode PhysicsRiemann_Burgers_Roe(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
436: {
437: PetscReal speed;
439: PetscFunctionBeginUser;
440: speed = 0.5 * (uL[0] + uR[0]);
441: flux[0] = 0.25 * (PetscSqr(uL[0]) + PetscSqr(uR[0])) - 0.5 * PetscAbs(speed) * (uR[0] - uL[0]);
442: if (uL[0] <= 0 && 0 <= uR[0]) flux[0] = 0; /* Entropy fix for sonic rarefaction */
443: *maxspeed = speed;
444: PetscFunctionReturn(PETSC_SUCCESS);
445: }
447: static PetscErrorCode PhysicsRiemann_Burgers_LxF(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
448: {
449: PetscReal c;
450: PetscScalar fL, fR;
452: PetscFunctionBeginUser;
453: c = ((BurgersCtx *)vctx)->lxf_speed;
454: fL = 0.5 * PetscSqr(uL[0]);
455: fR = 0.5 * PetscSqr(uR[0]);
456: flux[0] = 0.5 * (fL + fR) - 0.5 * c * (uR[0] - uL[0]);
457: *maxspeed = c;
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: static PetscErrorCode PhysicsRiemann_Burgers_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
462: {
463: PetscReal c;
464: PetscScalar fL, fR;
466: PetscFunctionBeginUser;
467: c = PetscMax(PetscAbs(uL[0]), PetscAbs(uR[0]));
468: fL = 0.5 * PetscSqr(uL[0]);
469: fR = 0.5 * PetscSqr(uR[0]);
470: flux[0] = 0.5 * (fL + fR) - 0.5 * c * (uR[0] - uL[0]);
471: *maxspeed = c;
472: PetscFunctionReturn(PETSC_SUCCESS);
473: }
475: static PetscErrorCode PhysicsCreate_Burgers(FVCtx *ctx)
476: {
477: BurgersCtx *user;
478: RiemannFunction r;
479: PetscFunctionList rlist = NULL;
480: char rname[256] = "exact";
482: PetscFunctionBeginUser;
483: PetscCall(PetscNew(&user));
485: ctx->physics.sample = PhysicsSample_Burgers;
486: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
487: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
488: ctx->physics.user = user;
489: ctx->physics.dof = 1;
491: PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
492: PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Burgers_Exact));
493: PetscCall(RiemannListAdd(&rlist, "roe", PhysicsRiemann_Burgers_Roe));
494: PetscCall(RiemannListAdd(&rlist, "lxf", PhysicsRiemann_Burgers_LxF));
495: PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_Burgers_Rusanov));
496: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
497: {
498: PetscCall(PetscOptionsFList("-physics_burgers_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
499: }
500: PetscOptionsEnd();
501: PetscCall(RiemannListFind(rlist, rname, &r));
502: PetscCall(PetscFunctionListDestroy(&rlist));
503: ctx->physics.riemann = r;
505: /* *
506: * Hack to deal with LxF in semi-discrete form
507: * max speed is 1 for the basic initial conditions (where |u| <= 1)
508: * */
509: if (r == PhysicsRiemann_Burgers_LxF) user->lxf_speed = 1;
510: PetscFunctionReturn(PETSC_SUCCESS);
511: }
513: /* --------------------------------- Traffic ----------------------------------- */
515: typedef struct {
516: PetscReal lxf_speed;
517: PetscReal a;
518: } TrafficCtx;
520: static inline PetscScalar TrafficFlux(PetscScalar a, PetscScalar u)
521: {
522: return a * u * (1 - u);
523: }
525: static PetscErrorCode PhysicsSample_Traffic(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
526: {
527: PetscReal a = ((TrafficCtx *)vctx)->a;
529: PetscFunctionBeginUser;
530: PetscCheck(bctype != FVBC_PERIODIC || t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solution not implemented for periodic");
531: switch (initial) {
532: case 0:
533: u[0] = (-a * t < x) ? 2 : 0;
534: break;
535: case 1:
536: if (x < PetscMin(2 * a * t, 0.5 + a * t)) u[0] = -1;
537: else if (x < 1) u[0] = 0;
538: else u[0] = 1;
539: break;
540: case 2:
541: PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Only initial condition available");
542: u[0] = 0.7 + 0.3 * PetscSinReal(2 * PETSC_PI * ((x - xmin) / (xmax - xmin)));
543: break;
544: default:
545: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
546: }
547: PetscFunctionReturn(PETSC_SUCCESS);
548: }
550: static PetscErrorCode PhysicsRiemann_Traffic_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
551: {
552: PetscReal a = ((TrafficCtx *)vctx)->a;
554: PetscFunctionBeginUser;
555: if (uL[0] < uR[0]) {
556: flux[0] = PetscMin(TrafficFlux(a, uL[0]), TrafficFlux(a, uR[0]));
557: } else {
558: flux[0] = (uR[0] < 0.5 && 0.5 < uL[0]) ? TrafficFlux(a, 0.5) : PetscMax(TrafficFlux(a, uL[0]), TrafficFlux(a, uR[0]));
559: }
560: *maxspeed = a * MaxAbs(1 - 2 * uL[0], 1 - 2 * uR[0]);
561: PetscFunctionReturn(PETSC_SUCCESS);
562: }
564: static PetscErrorCode PhysicsRiemann_Traffic_Roe(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
565: {
566: PetscReal a = ((TrafficCtx *)vctx)->a;
567: PetscReal speed;
569: PetscFunctionBeginUser;
570: speed = a * (1 - (uL[0] + uR[0]));
571: flux[0] = 0.5 * (TrafficFlux(a, uL[0]) + TrafficFlux(a, uR[0])) - 0.5 * PetscAbs(speed) * (uR[0] - uL[0]);
572: *maxspeed = speed;
573: PetscFunctionReturn(PETSC_SUCCESS);
574: }
576: static PetscErrorCode PhysicsRiemann_Traffic_LxF(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
577: {
578: TrafficCtx *phys = (TrafficCtx *)vctx;
579: PetscReal a = phys->a;
580: PetscReal speed;
582: PetscFunctionBeginUser;
583: speed = a * (1 - (uL[0] + uR[0]));
584: flux[0] = 0.5 * (TrafficFlux(a, uL[0]) + TrafficFlux(a, uR[0])) - 0.5 * phys->lxf_speed * (uR[0] - uL[0]);
585: *maxspeed = speed;
586: PetscFunctionReturn(PETSC_SUCCESS);
587: }
589: static PetscErrorCode PhysicsRiemann_Traffic_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
590: {
591: PetscReal a = ((TrafficCtx *)vctx)->a;
592: PetscReal speed;
594: PetscFunctionBeginUser;
595: speed = a * PetscMax(PetscAbs(1 - 2 * uL[0]), PetscAbs(1 - 2 * uR[0]));
596: flux[0] = 0.5 * (TrafficFlux(a, uL[0]) + TrafficFlux(a, uR[0])) - 0.5 * speed * (uR[0] - uL[0]);
597: *maxspeed = speed;
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: static PetscErrorCode PhysicsCreate_Traffic(FVCtx *ctx)
602: {
603: TrafficCtx *user;
604: RiemannFunction r;
605: PetscFunctionList rlist = NULL;
606: char rname[256] = "exact";
608: PetscFunctionBeginUser;
609: PetscCall(PetscNew(&user));
610: ctx->physics.sample = PhysicsSample_Traffic;
611: ctx->physics.characteristic = PhysicsCharacteristic_Conservative;
612: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
613: ctx->physics.user = user;
614: ctx->physics.dof = 1;
616: PetscCall(PetscStrallocpy("density", &ctx->physics.fieldname[0]));
617: user->a = 0.5;
618: PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Traffic_Exact));
619: PetscCall(RiemannListAdd(&rlist, "roe", PhysicsRiemann_Traffic_Roe));
620: PetscCall(RiemannListAdd(&rlist, "lxf", PhysicsRiemann_Traffic_LxF));
621: PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_Traffic_Rusanov));
622: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Traffic", "");
623: PetscCall(PetscOptionsReal("-physics_traffic_a", "Flux = a*u*(1-u)", "", user->a, &user->a, NULL));
624: PetscCall(PetscOptionsFList("-physics_traffic_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
625: PetscOptionsEnd();
627: PetscCall(RiemannListFind(rlist, rname, &r));
628: PetscCall(PetscFunctionListDestroy(&rlist));
630: ctx->physics.riemann = r;
632: /* *
633: * Hack to deal with LxF in semi-discrete form
634: * max speed is 3*a for the basic initial conditions (-1 <= u <= 2)
635: * */
636: if (r == PhysicsRiemann_Traffic_LxF) user->lxf_speed = 3 * user->a;
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: /* --------------------------------- Linear Acoustics ----------------------------------- */
642: /* Flux: u_t + (A u)_x
643: * z = sqrt(rho*bulk), c = sqrt(rho/bulk)
644: * Spectral decomposition: A = R * D * Rinv
645: * [ cz] = [-z z] [-c ] [-1/2z 1/2]
646: * [c/z ] = [ 1 1] [ c] [ 1/2z 1/2]
647: *
648: * We decompose this into the left-traveling waves Al = R * D^- Rinv
649: * and the right-traveling waves Ar = R * D^+ * Rinv
650: * Multiplying out these expressions produces the following two matrices
651: */
653: typedef struct {
654: PetscReal c; /* speed of sound: c = sqrt(bulk/rho) */
655: PetscReal z; /* impedance: z = sqrt(rho*bulk) */
656: } AcousticsCtx;
658: PETSC_UNUSED static inline void AcousticsFlux(AcousticsCtx *ctx, const PetscScalar *u, PetscScalar *f)
659: {
660: f[0] = ctx->c * ctx->z * u[1];
661: f[1] = ctx->c / ctx->z * u[0];
662: }
664: static PetscErrorCode PhysicsCharacteristic_Acoustics(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
665: {
666: AcousticsCtx *phys = (AcousticsCtx *)vctx;
667: PetscReal z = phys->z, c = phys->c;
669: PetscFunctionBeginUser;
670: X[0 * 2 + 0] = -z;
671: X[0 * 2 + 1] = z;
672: X[1 * 2 + 0] = 1;
673: X[1 * 2 + 1] = 1;
674: Xi[0 * 2 + 0] = -1. / (2 * z);
675: Xi[0 * 2 + 1] = 1. / 2;
676: Xi[1 * 2 + 0] = 1. / (2 * z);
677: Xi[1 * 2 + 1] = 1. / 2;
678: speeds[0] = -c;
679: speeds[1] = c;
680: PetscFunctionReturn(PETSC_SUCCESS);
681: }
683: static PetscErrorCode PhysicsSample_Acoustics_Initial(AcousticsCtx *phys, PetscInt initial, PetscReal xmin, PetscReal xmax, PetscReal x, PetscReal *u)
684: {
685: PetscFunctionBeginUser;
686: switch (initial) {
687: case 0:
688: u[0] = (PetscAbs((x - xmin) / (xmax - xmin) - 0.2) < 0.1) ? 1 : 0.5;
689: u[1] = (PetscAbs((x - xmin) / (xmax - xmin) - 0.7) < 0.1) ? 1 : -0.5;
690: break;
691: case 1:
692: u[0] = PetscCosReal(3 * 2 * PETSC_PI * x / (xmax - xmin));
693: u[1] = PetscExpReal(-PetscSqr(x - (xmax + xmin) / 2) / (2 * PetscSqr(0.2 * (xmax - xmin)))) - 0.5;
694: break;
695: default:
696: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
697: }
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode PhysicsSample_Acoustics(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
702: {
703: AcousticsCtx *phys = (AcousticsCtx *)vctx;
704: PetscReal c = phys->c;
705: PetscReal x0a, x0b, u0a[2], u0b[2], tmp[2];
706: PetscReal X[2][2], Xi[2][2], dummy[2];
708: PetscFunctionBeginUser;
709: switch (bctype) {
710: case FVBC_OUTFLOW:
711: x0a = x + c * t;
712: x0b = x - c * t;
713: break;
714: case FVBC_PERIODIC:
715: x0a = RangeMod(x + c * t, xmin, xmax);
716: x0b = RangeMod(x - c * t, xmin, xmax);
717: break;
718: default:
719: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
720: }
721: PetscCall(PhysicsSample_Acoustics_Initial(phys, initial, xmin, xmax, x0a, u0a));
722: PetscCall(PhysicsSample_Acoustics_Initial(phys, initial, xmin, xmax, x0b, u0b));
723: PetscCall(PhysicsCharacteristic_Acoustics(vctx, 2, u, &X[0][0], &Xi[0][0], dummy));
724: tmp[0] = Xi[0][0] * u0a[0] + Xi[0][1] * u0a[1];
725: tmp[1] = Xi[1][0] * u0b[0] + Xi[1][1] * u0b[1];
726: u[0] = X[0][0] * tmp[0] + X[0][1] * tmp[1];
727: u[1] = X[1][0] * tmp[0] + X[1][1] * tmp[1];
728: PetscFunctionReturn(PETSC_SUCCESS);
729: }
731: static PetscErrorCode PhysicsRiemann_Acoustics_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
732: {
733: AcousticsCtx *phys = (AcousticsCtx *)vctx;
734: PetscReal c = phys->c, z = phys->z;
735: PetscReal Al[2][2] =
736: {
737: {-c / 2, c * z / 2},
738: {c / (2 * z), -c / 2 }
739: }, /* Left traveling waves */
740: Ar[2][2] = {{c / 2, c * z / 2}, {c / (2 * z), c / 2}}; /* Right traveling waves */
742: PetscFunctionBeginUser;
743: flux[0] = Al[0][0] * uR[0] + Al[0][1] * uR[1] + Ar[0][0] * uL[0] + Ar[0][1] * uL[1];
744: flux[1] = Al[1][0] * uR[0] + Al[1][1] * uR[1] + Ar[1][0] * uL[0] + Ar[1][1] * uL[1];
745: *maxspeed = c;
746: PetscFunctionReturn(PETSC_SUCCESS);
747: }
749: static PetscErrorCode PhysicsCreate_Acoustics(FVCtx *ctx)
750: {
751: AcousticsCtx *user;
752: PetscFunctionList rlist = NULL, rclist = NULL;
753: char rname[256] = "exact", rcname[256] = "characteristic";
755: PetscFunctionBeginUser;
756: PetscCall(PetscNew(&user));
757: ctx->physics.sample = PhysicsSample_Acoustics;
758: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
759: ctx->physics.user = user;
760: ctx->physics.dof = 2;
762: PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
763: PetscCall(PetscStrallocpy("v", &ctx->physics.fieldname[1]));
765: user->c = 1;
766: user->z = 1;
768: PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Acoustics_Exact));
769: PetscCall(ReconstructListAdd(&rclist, "characteristic", PhysicsCharacteristic_Acoustics));
770: PetscCall(ReconstructListAdd(&rclist, "conservative", PhysicsCharacteristic_Conservative));
771: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for linear Acoustics", "");
772: {
773: PetscCall(PetscOptionsReal("-physics_acoustics_c", "c = sqrt(bulk/rho)", "", user->c, &user->c, NULL));
774: PetscCall(PetscOptionsReal("-physics_acoustics_z", "z = sqrt(bulk*rho)", "", user->z, &user->z, NULL));
775: PetscCall(PetscOptionsFList("-physics_acoustics_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
776: PetscCall(PetscOptionsFList("-physics_acoustics_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
777: }
778: PetscOptionsEnd();
779: PetscCall(RiemannListFind(rlist, rname, &ctx->physics.riemann));
780: PetscCall(ReconstructListFind(rclist, rcname, &ctx->physics.characteristic));
781: PetscCall(PetscFunctionListDestroy(&rlist));
782: PetscCall(PetscFunctionListDestroy(&rclist));
783: PetscFunctionReturn(PETSC_SUCCESS);
784: }
786: /* --------------------------------- Isothermal Gas Dynamics ----------------------------------- */
788: typedef struct {
789: PetscReal acoustic_speed;
790: } IsoGasCtx;
792: static inline void IsoGasFlux(PetscReal c, const PetscScalar *u, PetscScalar *f)
793: {
794: f[0] = u[1];
795: f[1] = PetscSqr(u[1]) / u[0] + c * c * u[0];
796: }
798: static PetscErrorCode PhysicsSample_IsoGas(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
799: {
800: PetscFunctionBeginUser;
801: PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solutions not implemented for t > 0");
802: switch (initial) {
803: case 0:
804: u[0] = (x < 0) ? 1 : 0.5;
805: u[1] = (x < 0) ? 1 : 0.7;
806: break;
807: case 1:
808: u[0] = 1 + 0.5 * PetscSinReal(2 * PETSC_PI * x);
809: u[1] = 1 * u[0];
810: break;
811: default:
812: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
813: }
814: PetscFunctionReturn(PETSC_SUCCESS);
815: }
817: static PetscErrorCode PhysicsRiemann_IsoGas_Roe(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
818: {
819: IsoGasCtx *phys = (IsoGasCtx *)vctx;
820: PetscReal c = phys->acoustic_speed;
821: PetscScalar ubar, du[2], a[2], fL[2], fR[2], lam[2], ustar[2], R[2][2];
822: PetscInt i;
824: PetscFunctionBeginUser;
825: ubar = (uL[1] / PetscSqrtScalar(uL[0]) + uR[1] / PetscSqrtScalar(uR[0])) / (PetscSqrtScalar(uL[0]) + PetscSqrtScalar(uR[0]));
826: /* write fluxuations in characteristic basis */
827: du[0] = uR[0] - uL[0];
828: du[1] = uR[1] - uL[1];
829: a[0] = (1 / (2 * c)) * ((ubar + c) * du[0] - du[1]);
830: a[1] = (1 / (2 * c)) * ((-ubar + c) * du[0] + du[1]);
831: /* wave speeds */
832: lam[0] = ubar - c;
833: lam[1] = ubar + c;
834: /* Right eigenvectors */
835: R[0][0] = 1;
836: R[0][1] = ubar - c;
837: R[1][0] = 1;
838: R[1][1] = ubar + c;
839: /* Compute state in star region (between the 1-wave and 2-wave) */
840: for (i = 0; i < 2; i++) ustar[i] = uL[i] + a[0] * R[0][i];
841: if (uL[1] / uL[0] < c && c < ustar[1] / ustar[0]) { /* 1-wave is sonic rarefaction */
842: PetscScalar ufan[2];
843: ufan[0] = uL[0] * PetscExpScalar(uL[1] / (uL[0] * c) - 1);
844: ufan[1] = c * ufan[0];
845: IsoGasFlux(c, ufan, flux);
846: } else if (ustar[1] / ustar[0] < -c && -c < uR[1] / uR[0]) { /* 2-wave is sonic rarefaction */
847: PetscScalar ufan[2];
848: ufan[0] = uR[0] * PetscExpScalar(-uR[1] / (uR[0] * c) - 1);
849: ufan[1] = -c * ufan[0];
850: IsoGasFlux(c, ufan, flux);
851: } else { /* Centered form */
852: IsoGasFlux(c, uL, fL);
853: IsoGasFlux(c, uR, fR);
854: for (i = 0; i < 2; i++) {
855: PetscScalar absdu = PetscAbsScalar(lam[0]) * a[0] * R[0][i] + PetscAbsScalar(lam[1]) * a[1] * R[1][i];
856: flux[i] = 0.5 * (fL[i] + fR[i]) - 0.5 * absdu;
857: }
858: }
859: *maxspeed = MaxAbs(lam[0], lam[1]);
860: PetscFunctionReturn(PETSC_SUCCESS);
861: }
863: static PetscErrorCode PhysicsRiemann_IsoGas_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
864: {
865: IsoGasCtx *phys = (IsoGasCtx *)vctx;
866: PetscReal c = phys->acoustic_speed;
867: PetscScalar ustar[2];
868: struct {
869: PetscScalar rho, u;
870: } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}, star;
871: PetscInt i;
873: PetscFunctionBeginUser;
874: PetscCheck(L.rho > 0 && R.rho > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
875: {
876: /* Solve for star state */
877: PetscScalar res, tmp, rho = 0.5 * (L.rho + R.rho); /* initial guess */
878: for (i = 0; i < 20; i++) {
879: PetscScalar fr, fl, dfr, dfl;
880: fl = (L.rho < rho) ? (rho - L.rho) / PetscSqrtScalar(L.rho * rho) /* shock */
881: : PetscLogScalar(rho) - PetscLogScalar(L.rho); /* rarefaction */
882: fr = (R.rho < rho) ? (rho - R.rho) / PetscSqrtScalar(R.rho * rho) /* shock */
883: : PetscLogScalar(rho) - PetscLogScalar(R.rho); /* rarefaction */
884: res = R.u - L.u + c * (fr + fl);
885: PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
886: if (PetscAbsScalar(res) < 1e-10) {
887: star.rho = rho;
888: star.u = L.u - c * fl;
889: goto converged;
890: }
891: dfl = (L.rho < rho) ? 1 / PetscSqrtScalar(L.rho * rho) * (1 - 0.5 * (rho - L.rho) / rho) : 1 / rho;
892: dfr = (R.rho < rho) ? 1 / PetscSqrtScalar(R.rho * rho) * (1 - 0.5 * (rho - R.rho) / rho) : 1 / rho;
893: tmp = rho - res / (c * (dfr + dfl));
894: if (tmp <= 0) rho /= 2; /* Guard against Newton shooting off to a negative density */
895: else rho = tmp;
896: PetscCheck((rho > 0) && PetscIsNormalScalar(rho), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate rho=%g", (double)PetscRealPart(rho));
897: }
898: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Newton iteration for star.rho diverged after %" PetscInt_FMT " iterations", i);
899: }
900: converged:
901: if (L.u - c < 0 && 0 < star.u - c) { /* 1-wave is sonic rarefaction */
902: PetscScalar ufan[2];
903: ufan[0] = L.rho * PetscExpScalar(L.u / c - 1);
904: ufan[1] = c * ufan[0];
905: IsoGasFlux(c, ufan, flux);
906: } else if (star.u + c < 0 && 0 < R.u + c) { /* 2-wave is sonic rarefaction */
907: PetscScalar ufan[2];
908: ufan[0] = R.rho * PetscExpScalar(-R.u / c - 1);
909: ufan[1] = -c * ufan[0];
910: IsoGasFlux(c, ufan, flux);
911: } else if ((L.rho >= star.rho && L.u - c >= 0) || (L.rho < star.rho && (star.rho * star.u - L.rho * L.u) / (star.rho - L.rho) > 0)) {
912: /* 1-wave is supersonic rarefaction, or supersonic shock */
913: IsoGasFlux(c, uL, flux);
914: } else if ((star.rho <= R.rho && R.u + c <= 0) || (star.rho > R.rho && (R.rho * R.u - star.rho * star.u) / (R.rho - star.rho) < 0)) {
915: /* 2-wave is supersonic rarefaction or supersonic shock */
916: IsoGasFlux(c, uR, flux);
917: } else {
918: ustar[0] = star.rho;
919: ustar[1] = star.rho * star.u;
920: IsoGasFlux(c, ustar, flux);
921: }
922: *maxspeed = MaxAbs(MaxAbs(star.u - c, star.u + c), MaxAbs(L.u - c, R.u + c));
923: PetscFunctionReturn(PETSC_SUCCESS);
924: }
926: static PetscErrorCode PhysicsRiemann_IsoGas_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
927: {
928: IsoGasCtx *phys = (IsoGasCtx *)vctx;
929: PetscScalar c = phys->acoustic_speed, fL[2], fR[2], s;
930: struct {
931: PetscScalar rho, u;
932: } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]};
934: PetscFunctionBeginUser;
935: PetscCheck(L.rho > 0 && R.rho > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed density is negative");
936: IsoGasFlux(c, uL, fL);
937: IsoGasFlux(c, uR, fR);
938: s = PetscMax(PetscAbs(L.u), PetscAbs(R.u)) + c;
939: flux[0] = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (uL[0] - uR[0]);
940: flux[1] = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]);
941: *maxspeed = s;
942: PetscFunctionReturn(PETSC_SUCCESS);
943: }
945: static PetscErrorCode PhysicsCharacteristic_IsoGas(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
946: {
947: IsoGasCtx *phys = (IsoGasCtx *)vctx;
948: PetscReal c = phys->acoustic_speed;
950: PetscFunctionBeginUser;
951: speeds[0] = u[1] / u[0] - c;
952: speeds[1] = u[1] / u[0] + c;
953: X[0 * 2 + 0] = 1;
954: X[0 * 2 + 1] = speeds[0];
955: X[1 * 2 + 0] = 1;
956: X[1 * 2 + 1] = speeds[1];
957: PetscCall(PetscArraycpy(Xi, X, 4));
958: PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: static PetscErrorCode PhysicsCreate_IsoGas(FVCtx *ctx)
963: {
964: IsoGasCtx *user;
965: PetscFunctionList rlist = NULL, rclist = NULL;
966: char rname[256] = "exact", rcname[256] = "characteristic";
968: PetscFunctionBeginUser;
969: PetscCall(PetscNew(&user));
970: ctx->physics.sample = PhysicsSample_IsoGas;
971: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
972: ctx->physics.user = user;
973: ctx->physics.dof = 2;
975: PetscCall(PetscStrallocpy("density", &ctx->physics.fieldname[0]));
976: PetscCall(PetscStrallocpy("momentum", &ctx->physics.fieldname[1]));
978: user->acoustic_speed = 1;
980: PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_IsoGas_Exact));
981: PetscCall(RiemannListAdd(&rlist, "roe", PhysicsRiemann_IsoGas_Roe));
982: PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_IsoGas_Rusanov));
983: PetscCall(ReconstructListAdd(&rclist, "characteristic", PhysicsCharacteristic_IsoGas));
984: PetscCall(ReconstructListAdd(&rclist, "conservative", PhysicsCharacteristic_Conservative));
985: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for IsoGas", "");
986: PetscCall(PetscOptionsReal("-physics_isogas_acoustic_speed", "Acoustic speed", "", user->acoustic_speed, &user->acoustic_speed, NULL));
987: PetscCall(PetscOptionsFList("-physics_isogas_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
988: PetscCall(PetscOptionsFList("-physics_isogas_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
989: PetscOptionsEnd();
990: PetscCall(RiemannListFind(rlist, rname, &ctx->physics.riemann));
991: PetscCall(ReconstructListFind(rclist, rcname, &ctx->physics.characteristic));
992: PetscCall(PetscFunctionListDestroy(&rlist));
993: PetscCall(PetscFunctionListDestroy(&rclist));
994: PetscFunctionReturn(PETSC_SUCCESS);
995: }
997: /* --------------------------------- Shallow Water ----------------------------------- */
998: typedef struct {
999: PetscReal gravity;
1000: } ShallowCtx;
1002: static inline void ShallowFlux(ShallowCtx *phys, const PetscScalar *u, PetscScalar *f)
1003: {
1004: f[0] = u[1];
1005: f[1] = PetscSqr(u[1]) / u[0] + 0.5 * phys->gravity * PetscSqr(u[0]);
1006: }
1008: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
1009: {
1010: ShallowCtx *phys = (ShallowCtx *)vctx;
1011: PetscScalar g = phys->gravity, ustar[2], cL, cR, c, cstar;
1012: struct {
1013: PetscScalar h, u;
1014: } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}, star;
1015: PetscInt i;
1017: PetscFunctionBeginUser;
1018: PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
1019: cL = PetscSqrtScalar(g * L.h);
1020: cR = PetscSqrtScalar(g * R.h);
1021: c = PetscMax(cL, cR);
1022: {
1023: /* Solve for star state */
1024: const PetscInt maxits = 50;
1025: PetscScalar tmp, res, res0 = 0, h0, h = 0.5 * (L.h + R.h); /* initial guess */
1026: h0 = h;
1027: for (i = 0; i < maxits; i++) {
1028: PetscScalar fr, fl, dfr, dfl;
1029: fl = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */
1030: : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h); /* rarefaction */
1031: fr = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */
1032: : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h); /* rarefaction */
1033: res = R.u - L.u + fr + fl;
1034: PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
1035: if (PetscAbsScalar(res) < 1e-8 || (i > 0 && PetscAbsScalar(h - h0) < 1e-8)) {
1036: star.h = h;
1037: star.u = L.u - fl;
1038: goto converged;
1039: } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
1040: h = 0.8 * h0 + 0.2 * h;
1041: continue;
1042: }
1043: /* Accept the last step and take another */
1044: res0 = res;
1045: h0 = h;
1046: dfl = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar(g / h);
1047: dfr = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar(g / h);
1048: tmp = h - res / (dfr + dfl);
1049: if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
1050: else h = tmp;
1051: PetscCheck((h > 0) && PetscIsNormalScalar(h), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h=%g", (double)h);
1052: }
1053: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Newton iteration for star.h diverged after %" PetscInt_FMT " iterations", i);
1054: }
1055: converged:
1056: cstar = PetscSqrtScalar(g * star.h);
1057: if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */
1058: PetscScalar ufan[2];
1059: ufan[0] = 1 / g * PetscSqr(L.u / 3 + 2. / 3 * cL);
1060: ufan[1] = PetscSqrtScalar(g * ufan[0]) * ufan[0];
1061: ShallowFlux(phys, ufan, flux);
1062: } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */
1063: PetscScalar ufan[2];
1064: ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR);
1065: ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0];
1066: ShallowFlux(phys, ufan, flux);
1067: } else if ((L.h >= star.h && L.u - c >= 0) || (L.h < star.h && (star.h * star.u - L.h * L.u) / (star.h - L.h) > 0)) {
1068: /* 1-wave is right-travelling shock (supersonic) */
1069: ShallowFlux(phys, uL, flux);
1070: } else if ((star.h <= R.h && R.u + c <= 0) || (star.h > R.h && (R.h * R.u - star.h * star.h) / (R.h - star.h) < 0)) {
1071: /* 2-wave is left-travelling shock (supersonic) */
1072: ShallowFlux(phys, uR, flux);
1073: } else {
1074: ustar[0] = star.h;
1075: ustar[1] = star.h * star.u;
1076: ShallowFlux(phys, ustar, flux);
1077: }
1078: *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR));
1079: PetscFunctionReturn(PETSC_SUCCESS);
1080: }
1082: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
1083: {
1084: ShallowCtx *phys = (ShallowCtx *)vctx;
1085: PetscScalar g = phys->gravity, fL[2], fR[2], s;
1086: struct {
1087: PetscScalar h, u;
1088: } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]};
1090: PetscFunctionBeginUser;
1091: PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
1092: ShallowFlux(phys, uL, fL);
1093: ShallowFlux(phys, uR, fR);
1094: s = PetscMax(PetscAbs(L.u) + PetscSqrtScalar(g * L.h), PetscAbs(R.u) + PetscSqrtScalar(g * R.h));
1095: flux[0] = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (uL[0] - uR[0]);
1096: flux[1] = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]);
1097: *maxspeed = s;
1098: PetscFunctionReturn(PETSC_SUCCESS);
1099: }
1101: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
1102: {
1103: ShallowCtx *phys = (ShallowCtx *)vctx;
1104: PetscReal c;
1106: PetscFunctionBeginUser;
1107: c = PetscSqrtScalar(u[0] * phys->gravity);
1108: speeds[0] = u[1] / u[0] - c;
1109: speeds[1] = u[1] / u[0] + c;
1110: X[0 * 2 + 0] = 1;
1111: X[0 * 2 + 1] = speeds[0];
1112: X[1 * 2 + 0] = 1;
1113: X[1 * 2 + 1] = speeds[1];
1114: PetscCall(PetscArraycpy(Xi, X, 4));
1115: PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
1116: PetscFunctionReturn(PETSC_SUCCESS);
1117: }
1119: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
1120: {
1121: ShallowCtx *user;
1122: PetscFunctionList rlist = NULL, rclist = NULL;
1123: char rname[256] = "exact", rcname[256] = "characteristic";
1125: PetscFunctionBeginUser;
1126: PetscCall(PetscNew(&user));
1127: /* Shallow water and Isothermal Gas dynamics are similar so we reuse initial conditions for now */
1128: ctx->physics.sample = PhysicsSample_IsoGas;
1129: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
1130: ctx->physics.user = user;
1131: ctx->physics.dof = 2;
1133: PetscCall(PetscStrallocpy("density", &ctx->physics.fieldname[0]));
1134: PetscCall(PetscStrallocpy("momentum", &ctx->physics.fieldname[1]));
1136: user->gravity = 1;
1138: PetscCall(RiemannListAdd(&rlist, "exact", PhysicsRiemann_Shallow_Exact));
1139: PetscCall(RiemannListAdd(&rlist, "rusanov", PhysicsRiemann_Shallow_Rusanov));
1140: PetscCall(ReconstructListAdd(&rclist, "characteristic", PhysicsCharacteristic_Shallow));
1141: PetscCall(ReconstructListAdd(&rclist, "conservative", PhysicsCharacteristic_Conservative));
1142: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", "");
1143: PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravity, NULL));
1144: PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
1145: PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
1146: PetscOptionsEnd();
1147: PetscCall(RiemannListFind(rlist, rname, &ctx->physics.riemann));
1148: PetscCall(ReconstructListFind(rclist, rcname, &ctx->physics.characteristic));
1149: PetscCall(PetscFunctionListDestroy(&rlist));
1150: PetscCall(PetscFunctionListDestroy(&rclist));
1151: PetscFunctionReturn(PETSC_SUCCESS);
1152: }
1154: /* --------------------------------- Finite Volume Solver ----------------------------------- */
1156: static PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
1157: {
1158: FVCtx *ctx = (FVCtx *)vctx;
1159: PetscInt i, j, k, Mx, dof, xs, xm;
1160: PetscReal hx;
1161: PetscScalar *x, *f, *slope;
1162: Vec Xloc;
1163: DM da;
1165: PetscFunctionBeginUser;
1166: ctx->cfl_idt = 0;
1167: PetscCall(TSGetDM(ts, &da));
1168: PetscCall(DMGetLocalVector(da, &Xloc));
1169: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1170: hx = (ctx->xmax - ctx->xmin) / Mx;
1171: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
1172: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
1174: PetscCall(VecZeroEntries(F));
1176: PetscCall(DMDAVecGetArray(da, Xloc, &x));
1177: PetscCall(DMDAVecGetArray(da, F, &f));
1178: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
1180: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1182: if (ctx->bctype == FVBC_OUTFLOW) {
1183: for (i = xs - 2; i < 0; i++) {
1184: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
1185: }
1186: for (i = Mx; i < xs + xm + 2; i++) {
1187: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
1188: }
1189: }
1190: for (i = xs - 1; i < xs + xm + 1; i++) {
1191: struct _LimitInfo info;
1192: PetscScalar *cjmpL, *cjmpR;
1193: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
1194: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
1195: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
1196: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
1197: cjmpL = &ctx->cjmpLR[0];
1198: cjmpR = &ctx->cjmpLR[dof];
1199: for (j = 0; j < dof; j++) {
1200: PetscScalar jmpL, jmpR;
1201: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
1202: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
1203: for (k = 0; k < dof; k++) {
1204: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
1205: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
1206: }
1207: }
1208: /* Apply limiter to the left and right characteristic jumps */
1209: info.m = dof;
1210: info.hx = hx;
1211: (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
1212: for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
1213: for (j = 0; j < dof; j++) {
1214: PetscScalar tmp = 0;
1215: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
1216: slope[i * dof + j] = tmp;
1217: }
1218: }
1220: for (i = xs; i < xs + xm + 1; i++) {
1221: PetscReal maxspeed;
1222: PetscScalar *uL, *uR;
1223: uL = &ctx->uLR[0];
1224: uR = &ctx->uLR[dof];
1225: for (j = 0; j < dof; j++) {
1226: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
1227: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
1228: }
1229: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed));
1230: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
1232: if (i > xs) {
1233: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
1234: }
1235: if (i < xs + xm) {
1236: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx;
1237: }
1238: }
1240: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
1241: PetscCall(DMDAVecRestoreArray(da, F, &f));
1242: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
1243: PetscCall(DMRestoreLocalVector(da, &Xloc));
1245: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)da)));
1246: if (0) {
1247: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
1248: PetscReal dt, tnow;
1249: PetscCall(TSGetTimeStep(ts, &dt));
1250: PetscCall(TSGetTime(ts, &tnow));
1251: if (dt > 0.5 / ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
1252: }
1253: PetscFunctionReturn(PETSC_SUCCESS);
1254: }
1256: static PetscErrorCode SmallMatMultADB(PetscScalar *C, PetscInt bs, const PetscScalar *A, const PetscReal *D, const PetscScalar *B)
1257: {
1258: PetscInt i, j, k;
1260: PetscFunctionBeginUser;
1261: for (i = 0; i < bs; i++) {
1262: for (j = 0; j < bs; j++) {
1263: PetscScalar tmp = 0;
1264: for (k = 0; k < bs; k++) tmp += A[i * bs + k] * D[k] * B[k * bs + j];
1265: C[i * bs + j] = tmp;
1266: }
1267: }
1268: PetscFunctionReturn(PETSC_SUCCESS);
1269: }
1271: static PetscErrorCode FVIJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal shift, Mat A, Mat B, void *vctx)
1272: {
1273: FVCtx *ctx = (FVCtx *)vctx;
1274: PetscInt i, j, dof = ctx->physics.dof;
1275: PetscScalar *J;
1276: const PetscScalar *x;
1277: PetscReal hx;
1278: DM da;
1279: DMDALocalInfo dainfo;
1281: PetscFunctionBeginUser;
1282: PetscCall(TSGetDM(ts, &da));
1283: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&x));
1284: PetscCall(DMDAGetLocalInfo(da, &dainfo));
1285: hx = (ctx->xmax - ctx->xmin) / dainfo.mx;
1286: PetscCall(PetscMalloc1(dof * dof, &J));
1287: for (i = dainfo.xs; i < dainfo.xs + dainfo.xm; i++) {
1288: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
1289: for (j = 0; j < dof; j++) ctx->speeds[j] = PetscAbs(ctx->speeds[j]);
1290: PetscCall(SmallMatMultADB(J, dof, ctx->R, ctx->speeds, ctx->Rinv));
1291: for (j = 0; j < dof * dof; j++) J[j] = J[j] / hx + shift * (j / dof == j % dof);
1292: PetscCall(MatSetValuesBlocked(B, 1, &i, 1, &i, J, INSERT_VALUES));
1293: }
1294: PetscCall(PetscFree(J));
1295: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&x));
1297: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1298: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
1299: if (A != B) {
1300: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1301: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
1302: }
1303: PetscFunctionReturn(PETSC_SUCCESS);
1304: }
1306: static PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
1307: {
1308: PetscScalar *u, *uj;
1309: PetscInt i, j, k, dof, xs, xm, Mx;
1311: PetscFunctionBeginUser;
1312: PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
1313: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1314: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1315: PetscCall(DMDAVecGetArray(da, U, &u));
1316: PetscCall(PetscMalloc1(dof, &uj));
1317: for (i = xs; i < xs + xm; i++) {
1318: const PetscReal h = (ctx->xmax - ctx->xmin) / Mx, xi = ctx->xmin + h / 2 + i * h;
1319: const PetscInt N = 200;
1320: /* Integrate over cell i using trapezoid rule with N points. */
1321: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
1322: for (j = 0; j < N + 1; j++) {
1323: PetscScalar xj = xi + h * (j - N / 2) / (PetscReal)N;
1324: PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
1325: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
1326: }
1327: }
1328: PetscCall(DMDAVecRestoreArray(da, U, &u));
1329: PetscCall(PetscFree(uj));
1330: PetscFunctionReturn(PETSC_SUCCESS);
1331: }
1333: static PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
1334: {
1335: PetscReal xmin, xmax;
1336: PetscScalar sum, tvsum;
1337: const PetscScalar *x;
1338: PetscInt imin, imax, Mx, i, j, xs, xm, dof;
1339: Vec Xloc;
1340: PetscBool isascii;
1342: PetscFunctionBeginUser;
1343: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1344: PetscCheck(isascii, PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
1345: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
1346: PetscCall(DMGetLocalVector(da, &Xloc));
1347: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
1348: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
1349: PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
1350: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1351: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1352: tvsum = 0;
1353: for (i = xs; i < xs + xm; i++) {
1354: for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
1355: }
1356: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &tvsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)da)));
1357: PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
1358: PetscCall(DMRestoreLocalVector(da, &Xloc));
1360: PetscCall(VecMin(X, &imin, &xmin));
1361: PetscCall(VecMax(X, &imax, &xmax));
1362: PetscCall(VecSum(X, &sum));
1363: PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%8.5f,%8.5f] with extrema at %" PetscInt_FMT " and %" PetscInt_FMT ", mean %8.5f, ||x||_TV %8.5f\n", (double)xmin, (double)xmax, imin, imax, (double)(sum / Mx), (double)(tvsum / Mx)));
1364: PetscFunctionReturn(PETSC_SUCCESS);
1365: }
1367: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1, PetscReal *nrmsup)
1368: {
1369: Vec Y;
1370: PetscInt Mx;
1372: PetscFunctionBeginUser;
1373: PetscCall(VecGetSize(X, &Mx));
1374: PetscCall(VecDuplicate(X, &Y));
1375: PetscCall(FVSample(ctx, da, t, Y));
1376: PetscCall(VecAYPX(Y, -1, X));
1377: PetscCall(VecNorm(Y, NORM_1, nrm1));
1378: PetscCall(VecNorm(Y, NORM_INFINITY, nrmsup));
1379: *nrm1 /= Mx;
1380: PetscCall(VecDestroy(&Y));
1381: PetscFunctionReturn(PETSC_SUCCESS);
1382: }
1384: int main(int argc, char *argv[])
1385: {
1386: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
1387: PetscFunctionList limiters = NULL, physics = NULL;
1388: MPI_Comm comm;
1389: TS ts;
1390: DM da;
1391: Vec X, X0, R;
1392: Mat B;
1393: FVCtx ctx;
1394: PetscInt i, dof, xs, xm, Mx, draw = 0;
1395: PetscBool view_final = PETSC_FALSE;
1396: PetscReal ptime;
1398: PetscFunctionBeginUser;
1399: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1400: comm = PETSC_COMM_WORLD;
1401: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
1403: /* Register limiters to be available on the command line */
1404: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit_Upwind));
1405: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit_LaxWendroff));
1406: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit_BeamWarming));
1407: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit_Fromm));
1408: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit_Minmod));
1409: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit_Superbee));
1410: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit_MC));
1411: PetscCall(PetscFunctionListAdd(&limiters, "vanleer", Limit_VanLeer));
1412: PetscCall(PetscFunctionListAdd(&limiters, "vanalbada", Limit_VanAlbada));
1413: PetscCall(PetscFunctionListAdd(&limiters, "vanalbadatvd", Limit_VanAlbadaTVD));
1414: PetscCall(PetscFunctionListAdd(&limiters, "koren", Limit_Koren));
1415: PetscCall(PetscFunctionListAdd(&limiters, "korensym", Limit_KorenSym));
1416: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit_Koren3));
1417: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon2", Limit_CadaTorrilhon2));
1418: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r0p1", Limit_CadaTorrilhon3R0p1));
1419: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r1", Limit_CadaTorrilhon3R1));
1420: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r10", Limit_CadaTorrilhon3R10));
1421: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r100", Limit_CadaTorrilhon3R100));
1423: /* Register physical models to be available on the command line */
1424: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
1425: PetscCall(PetscFunctionListAdd(&physics, "burgers", PhysicsCreate_Burgers));
1426: PetscCall(PetscFunctionListAdd(&physics, "traffic", PhysicsCreate_Traffic));
1427: PetscCall(PetscFunctionListAdd(&physics, "acoustics", PhysicsCreate_Acoustics));
1428: PetscCall(PetscFunctionListAdd(&physics, "isogas", PhysicsCreate_IsoGas));
1429: PetscCall(PetscFunctionListAdd(&physics, "shallow", PhysicsCreate_Shallow));
1431: ctx.comm = comm;
1432: ctx.cfl = 0.9;
1433: ctx.bctype = FVBC_PERIODIC;
1434: ctx.xmin = -1;
1435: ctx.xmax = 1;
1436: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
1437: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
1438: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
1439: PetscCall(PetscOptionsFList("-limit", "Name of flux limiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
1440: PetscCall(PetscOptionsFList("-physics", "Name of physics (Riemann solver and characteristics) to use", "", physics, physname, physname, sizeof(physname), NULL));
1441: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
1442: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
1443: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
1444: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
1445: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
1446: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
1447: PetscOptionsEnd();
1449: /* Choose the limiter from the list of registered limiters */
1450: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit));
1451: PetscCheck(ctx.limit, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
1453: /* Choose the physics from the list of registered models */
1454: {
1455: PetscErrorCode (*r)(FVCtx *);
1456: PetscCall(PetscFunctionListFind(physics, physname, &r));
1457: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
1458: /* Create the physics, will set the number of fields and their names */
1459: PetscCall((*r)(&ctx));
1460: }
1462: /* Create a DMDA to manage the parallel grid */
1463: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
1464: PetscCall(DMSetFromOptions(da));
1465: PetscCall(DMSetUp(da));
1466: /* Inform the DMDA of the field names provided by the physics. */
1467: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1468: for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
1469: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1470: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1472: /* Set coordinates of cell centers */
1473: PetscCall(DMDASetUniformCoordinates(da, ctx.xmin + 0.5 * (ctx.xmax - ctx.xmin) / Mx, ctx.xmax + 0.5 * (ctx.xmax - ctx.xmin) / Mx, 0, 0, 0, 0));
1475: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1476: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
1477: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
1479: /* Create a vector to store the solution and to save the initial state */
1480: PetscCall(DMCreateGlobalVector(da, &X));
1481: PetscCall(VecDuplicate(X, &X0));
1482: PetscCall(VecDuplicate(X, &R));
1484: PetscCall(DMCreateMatrix(da, &B));
1486: /* Create a time-stepping object */
1487: PetscCall(TSCreate(comm, &ts));
1488: PetscCall(TSSetDM(ts, da));
1489: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
1490: PetscCall(TSSetIJacobian(ts, B, B, FVIJacobian, &ctx));
1491: PetscCall(TSSetType(ts, TSSSP));
1492: PetscCall(TSSetMaxTime(ts, 10));
1493: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1495: /* Compute initial conditions and starting time step */
1496: PetscCall(FVSample(&ctx, da, 0, X0));
1497: PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
1498: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
1499: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
1500: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
1501: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1502: {
1503: PetscReal nrm1, nrmsup;
1504: PetscInt steps;
1506: PetscCall(TSSolve(ts, X));
1507: PetscCall(TSGetSolveTime(ts, &ptime));
1508: PetscCall(TSGetStepNumber(ts, &steps));
1510: PetscCall(PetscPrintf(comm, "Final time %8.5f, steps %" PetscInt_FMT "\n", (double)ptime, steps));
1511: if (ctx.exact) {
1512: PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1, &nrmsup));
1513: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %8.4e ||x-x_e||_sup %8.4e\n", (double)nrm1, (double)nrmsup));
1514: }
1515: }
1517: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1518: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
1519: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
1520: if (draw & 0x4) {
1521: Vec Y;
1522: PetscCall(VecDuplicate(X, &Y));
1523: PetscCall(FVSample(&ctx, da, ptime, Y));
1524: PetscCall(VecAYPX(Y, -1, X));
1525: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
1526: PetscCall(VecDestroy(&Y));
1527: }
1529: if (view_final) {
1530: PetscViewer viewer;
1531: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
1532: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
1533: PetscCall(VecView(X, viewer));
1534: PetscCall(PetscViewerPopFormat(viewer));
1535: PetscCall(PetscViewerDestroy(&viewer));
1536: }
1538: /* Clean up */
1539: PetscCall((*ctx.physics.destroy)(ctx.physics.user));
1540: for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
1541: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
1542: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
1543: PetscCall(VecDestroy(&X));
1544: PetscCall(VecDestroy(&X0));
1545: PetscCall(VecDestroy(&R));
1546: PetscCall(MatDestroy(&B));
1547: PetscCall(DMDestroy(&da));
1548: PetscCall(TSDestroy(&ts));
1549: PetscCall(PetscFunctionListDestroy(&limiters));
1550: PetscCall(PetscFunctionListDestroy(&physics));
1551: PetscCall(PetscFinalize());
1552: return 0;
1553: }
1555: /*TEST
1557: build:
1558: requires: !complex
1560: test:
1561: args: -da_grid_x 100 -initial 1 -xmin -2 -xmax 5 -exact -limit mc
1562: requires: !complex !single
1564: test:
1565: suffix: 2
1566: args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
1567: filter: sed "s/at 48/at 0/g"
1568: requires: !complex !single
1570: test:
1571: suffix: 3
1572: args: -da_grid_x 100 -initial 2 -xmin -2 -xmax 2 -exact -limit mc -physics burgers -bc_type outflow -ts_max_time 1
1573: nsize: 3
1574: filter: sed "s/at 48/at 0/g"
1575: requires: !complex !single
1577: TEST*/