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