Actual source code: ex4.c
1: /*
2: Note:
3: -hratio is the ratio between mesh size of coarse grids and fine grids
4: */
6: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
7: " advect - Constant coefficient scalar advection\n"
8: " u_t + (a*u)_x = 0\n"
9: " shallow - 1D Shallow water equations (Saint Venant System)\n"
10: " h_t + (q)_x = 0 \n"
11: " q_t + (\frac{q^2}{h} + g/2*h^2)_x = 0 \n"
12: " where, h(x,t) denotes the height of the water and q(x,t) the momentum.\n"
13: " for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
14: " hxs = hratio*hxf \n"
15: " where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
16: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
17: " the states across shocks and rarefactions\n"
18: " simulation - use reference solution which is generated by smaller time step size to be true solution,\n"
19: " also the reference solution should be generated by user and stored in a binary file.\n"
20: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
21: " bc_type - Boundary condition for the problem, options are: periodic, outflow, inflow "
22: "Several problem descriptions (initial data, physics specific features, boundary data) can be chosen with -initial N\n\n"
23: "The problem size should be set with -da_grid_x M\n\n";
25: /*
26: Example:
27: ./ex4 -da_grid_x 40 -initial 1 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
28: ./ex4 -da_grid_x 40 -initial 2 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 2.5 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
29: ./ex4 -da_grid_x 40 -initial 3 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
30: ./ex4 -da_grid_x 40 -initial 4 -hratio 1 -limit koren3 -ts_time_step 0.01 -ts_max_time 4.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
31: ./ex4 -da_grid_x 40 -initial 5 -hratio 1 -limit mc -ts_time_step 0.01 -ts_max_time 5.0 -ts_type mprk -ts_mprk_type 2a22 -ts_monitor_draw_solution -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
33: Contributed by: Aidan Hamilton <aidan@udel.edu>
34: */
36: #include <petscts.h>
37: #include <petscdm.h>
38: #include <petscdmda.h>
39: #include "finitevolume1d.h"
40: #include <petsc/private/kernels/blockinvert.h>
42: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
43: {
44: PetscReal range = xmax - xmin;
45: return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
46: }
47: static inline PetscReal MaxAbs(PetscReal a, PetscReal b)
48: {
49: return (PetscAbs(a) > PetscAbs(b)) ? a : b;
50: }
51: /* --------------------------------- Advection ----------------------------------- */
52: typedef struct {
53: PetscReal a; /* advective velocity */
54: } AdvectCtx;
56: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
57: {
58: AdvectCtx *ctx = (AdvectCtx *)vctx;
59: PetscReal speed;
61: PetscFunctionBeginUser;
62: speed = ctx->a;
63: flux[0] = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
64: *maxspeed = speed;
65: PetscFunctionReturn(PETSC_SUCCESS);
66: }
68: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
69: {
70: AdvectCtx *ctx = (AdvectCtx *)vctx;
72: PetscFunctionBeginUser;
73: X[0] = 1.;
74: Xi[0] = 1.;
75: speeds[0] = ctx->a;
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
80: {
81: AdvectCtx *ctx = (AdvectCtx *)vctx;
82: PetscReal a = ctx->a, x0;
84: PetscFunctionBeginUser;
85: switch (bctype) {
86: case FVBC_OUTFLOW:
87: x0 = x - a * t;
88: break;
89: case FVBC_PERIODIC:
90: x0 = RangeMod(x - a * t, xmin, xmax);
91: break;
92: default:
93: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
94: }
95: switch (initial) {
96: case 0:
97: u[0] = (x0 < 0) ? 1 : -1;
98: break;
99: case 1:
100: u[0] = (x0 < 0) ? -1 : 1;
101: break;
102: case 2:
103: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
104: break;
105: case 3:
106: u[0] = PetscSinReal(2 * PETSC_PI * x0);
107: break;
108: case 4:
109: u[0] = PetscAbs(x0);
110: break;
111: case 5:
112: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
113: break;
114: case 6:
115: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
116: break;
117: case 7:
118: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
119: break;
120: default:
121: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
122: }
123: PetscFunctionReturn(PETSC_SUCCESS);
124: }
126: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
127: {
128: AdvectCtx *user;
130: PetscFunctionBeginUser;
131: PetscCall(PetscNew(&user));
132: ctx->physics2.sample2 = PhysicsSample_Advect;
133: ctx->physics2.riemann2 = PhysicsRiemann_Advect;
134: ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
135: ctx->physics2.destroy = PhysicsDestroy_SimpleFree;
136: ctx->physics2.user = user;
137: ctx->physics2.dof = 1;
139: PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0]));
140: user->a = 1;
141: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
142: PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
143: PetscOptionsEnd();
144: PetscFunctionReturn(PETSC_SUCCESS);
145: }
146: /* --------------------------------- Shallow Water ----------------------------------- */
148: typedef struct {
149: PetscReal gravity;
150: } ShallowCtx;
152: static inline void ShallowFlux(ShallowCtx *phys, const PetscScalar *u, PetscScalar *f)
153: {
154: f[0] = u[1];
155: f[1] = PetscSqr(u[1]) / u[0] + 0.5 * phys->gravity * PetscSqr(u[0]);
156: }
158: static PetscErrorCode PhysicsRiemann_Shallow_Exact(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
159: {
160: ShallowCtx *phys = (ShallowCtx *)vctx;
161: PetscScalar g = phys->gravity, ustar[2], cL, cR, c, cstar;
162: struct {
163: PetscScalar h, u;
164: } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]}, star;
165: PetscInt i;
167: PetscFunctionBeginUser;
168: PetscCheck(L.h > 0 && R.h > 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Reconstructed thickness is negative");
169: cL = PetscSqrtScalar(g * L.h);
170: cR = PetscSqrtScalar(g * R.h);
171: c = PetscMax(cL, cR);
172: {
173: /* Solve for star state */
174: const PetscInt maxits = 50;
175: PetscScalar tmp, res, res0 = 0, h0, h = 0.5 * (L.h + R.h); /* initial guess */
176: h0 = h;
177: for (i = 0; i < maxits; i++) {
178: PetscScalar fr, fl, dfr, dfl;
179: fl = (L.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - L.h * L.h) * (1 / L.h - 1 / h)) /* shock */
180: : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * L.h); /* rarefaction */
181: fr = (R.h < h) ? PetscSqrtScalar(0.5 * g * (h * h - R.h * R.h) * (1 / R.h - 1 / h)) /* shock */
182: : 2 * PetscSqrtScalar(g * h) - 2 * PetscSqrtScalar(g * R.h); /* rarefaction */
183: res = R.u - L.u + fr + fl;
184: PetscCheck(!PetscIsInfOrNanScalar(res), PETSC_COMM_SELF, PETSC_ERR_FP, "Infinity or Not-a-Number generated in computation");
185: if (PetscAbsScalar(res) < PETSC_SQRT_MACHINE_EPSILON || (i > 0 && PetscAbsScalar(h - h0) < PETSC_SQRT_MACHINE_EPSILON)) {
186: star.h = h;
187: star.u = L.u - fl;
188: goto converged;
189: } else if (i > 0 && PetscAbsScalar(res) >= PetscAbsScalar(res0)) { /* Line search */
190: h = 0.8 * h0 + 0.2 * h;
191: continue;
192: }
193: /* Accept the last step and take another */
194: res0 = res;
195: h0 = h;
196: dfl = (L.h < h) ? 0.5 / fl * 0.5 * g * (-L.h * L.h / (h * h) - 1 + 2 * h / L.h) : PetscSqrtScalar(g / h);
197: dfr = (R.h < h) ? 0.5 / fr * 0.5 * g * (-R.h * R.h / (h * h) - 1 + 2 * h / R.h) : PetscSqrtScalar(g / h);
198: tmp = h - res / (dfr + dfl);
199: if (tmp <= 0) h /= 2; /* Guard against Newton shooting off to a negative thickness */
200: else h = tmp;
201: PetscCheck((h > 0) && PetscIsNormalScalar(h), PETSC_COMM_SELF, PETSC_ERR_FP, "non-normal iterate h=%g", (double)h);
202: }
203: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_NOT_CONVERGED, "Newton iteration for star.h diverged after %" PetscInt_FMT " iterations", i);
204: }
205: converged:
206: cstar = PetscSqrtScalar(g * star.h);
207: if (L.u - cL < 0 && 0 < star.u - cstar) { /* 1-wave is sonic rarefaction */
208: PetscScalar ufan[2];
209: ufan[0] = 1 / g * PetscSqr(L.u / 3 + 2. / 3 * cL);
210: ufan[1] = PetscSqrtScalar(g * ufan[0]) * ufan[0];
211: ShallowFlux(phys, ufan, flux);
212: } else if (star.u + cstar < 0 && 0 < R.u + cR) { /* 2-wave is sonic rarefaction */
213: PetscScalar ufan[2];
214: ufan[0] = 1 / g * PetscSqr(R.u / 3 - 2. / 3 * cR);
215: ufan[1] = -PetscSqrtScalar(g * ufan[0]) * ufan[0];
216: ShallowFlux(phys, ufan, flux);
217: } 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)) {
218: /* 1-wave is right-travelling shock (supersonic) */
219: ShallowFlux(phys, uL, flux);
220: } 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)) {
221: /* 2-wave is left-travelling shock (supersonic) */
222: ShallowFlux(phys, uR, flux);
223: } else {
224: ustar[0] = star.h;
225: ustar[1] = star.h * star.u;
226: ShallowFlux(phys, ustar, flux);
227: }
228: *maxspeed = MaxAbs(MaxAbs(star.u - cstar, star.u + cstar), MaxAbs(L.u - cL, R.u + cR));
229: PetscFunctionReturn(PETSC_SUCCESS);
230: }
232: static PetscErrorCode PhysicsRiemann_Shallow_Rusanov(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
233: {
234: ShallowCtx *phys = (ShallowCtx *)vctx;
235: PetscScalar g = phys->gravity, fL[2], fR[2], s;
236: struct {
237: PetscScalar h, u;
238: } L = {uL[0], uL[1] / uL[0]}, R = {uR[0], uR[1] / uR[0]};
239: PetscReal tol = 1e-6;
241: PetscFunctionBeginUser;
242: /* Positivity preserving modification*/
243: if (L.h < tol) L.u = 0.0;
244: if (R.h < tol) R.u = 0.0;
246: /*simple pos preserve limiter*/
247: if (L.h < 0) L.h = 0;
248: if (R.h < 0) R.h = 0;
250: ShallowFlux(phys, uL, fL);
251: ShallowFlux(phys, uR, fR);
253: s = PetscMax(PetscAbs(L.u) + PetscSqrtScalar(g * L.h), PetscAbs(R.u) + PetscSqrtScalar(g * R.h));
254: flux[0] = 0.5 * (fL[0] + fR[0]) + 0.5 * s * (L.h - R.h);
255: flux[1] = 0.5 * (fL[1] + fR[1]) + 0.5 * s * (uL[1] - uR[1]);
256: *maxspeed = s;
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
260: static PetscErrorCode PhysicsCharacteristic_Conservative(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
261: {
262: PetscInt i, j;
264: PetscFunctionBeginUser;
265: for (i = 0; i < m; i++) {
266: for (j = 0; j < m; j++) Xi[i * m + j] = X[i * m + j] = (PetscScalar)(i == j);
267: speeds[i] = PETSC_MAX_REAL; /* Indicates invalid */
268: }
269: PetscFunctionReturn(PETSC_SUCCESS);
270: }
272: static PetscErrorCode PhysicsCharacteristic_Shallow(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
273: {
274: ShallowCtx *phys = (ShallowCtx *)vctx;
275: PetscReal c;
276: PetscReal tol = 1e-6;
278: PetscFunctionBeginUser;
279: c = PetscSqrtScalar(u[0] * phys->gravity);
281: if (u[0] < tol) { /*Use conservative variables*/
282: X[0 * 2 + 0] = 1;
283: X[0 * 2 + 1] = 0;
284: X[1 * 2 + 0] = 0;
285: X[1 * 2 + 1] = 1;
286: } else {
287: speeds[0] = u[1] / u[0] - c;
288: speeds[1] = u[1] / u[0] + c;
289: X[0 * 2 + 0] = 1;
290: X[0 * 2 + 1] = speeds[0];
291: X[1 * 2 + 0] = 1;
292: X[1 * 2 + 1] = speeds[1];
293: }
295: PetscCall(PetscArraycpy(Xi, X, 4));
296: PetscCall(PetscKernel_A_gets_inverse_A_2(Xi, 0, PETSC_FALSE, NULL));
297: PetscFunctionReturn(PETSC_SUCCESS);
298: }
300: static PetscErrorCode PhysicsSample_Shallow(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
301: {
302: PetscFunctionBeginUser;
303: PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_SUP, "Exact solutions not implemented for t > 0");
304: switch (initial) {
305: case 0:
306: u[0] = (x < 0) ? 2 : 1; /* Standard Dam Break Problem */
307: u[1] = (x < 0) ? 0 : 0;
308: break;
309: case 1:
310: u[0] = (x < 10) ? 1 : 0.1; /*The Next 5 problems are standard Riemann problem tests */
311: u[1] = (x < 10) ? 2.5 : 0;
312: break;
313: case 2:
314: u[0] = (x < 25) ? 1 : 1;
315: u[1] = (x < 25) ? -5 : 5;
316: break;
317: case 3:
318: u[0] = (x < 20) ? 1 : 1e-12;
319: u[1] = (x < 20) ? 0 : 0;
320: break;
321: case 4:
322: u[0] = (x < 30) ? 1e-12 : 1;
323: u[1] = (x < 30) ? 0 : 0;
324: break;
325: case 5:
326: u[0] = (x < 25) ? 0.1 : 0.1;
327: u[1] = (x < 25) ? -0.3 : 0.3;
328: break;
329: case 6:
330: u[0] = 1 + 0.5 * PetscSinReal(2 * PETSC_PI * x);
331: u[1] = 1 * u[0];
332: break;
333: case 7:
334: if (x < -0.1) {
335: u[0] = 1e-9;
336: u[1] = 0.0;
337: } else if (x < 0.1) {
338: u[0] = 1.0;
339: u[1] = 0.0;
340: } else {
341: u[0] = 1e-9;
342: u[1] = 0.0;
343: }
344: break;
345: case 8:
346: if (x < -0.1) {
347: u[0] = 2;
348: u[1] = 0.0;
349: } else if (x < 0.1) {
350: u[0] = 3.0;
351: u[1] = 0.0;
352: } else {
353: u[0] = 2;
354: u[1] = 0.0;
355: }
356: break;
357: default:
358: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
359: }
360: PetscFunctionReturn(PETSC_SUCCESS);
361: }
363: /* Implements inflow conditions for the given initial conditions. Which conditions are actually enforced depends on
364: on the results of PhysicsSetInflowType_Shallow. */
365: static PetscErrorCode PhysicsInflow_Shallow(void *vctx, PetscReal t, PetscReal x, PetscReal *u)
366: {
367: FVCtx *ctx = (FVCtx *)vctx;
369: PetscFunctionBeginUser;
370: if (ctx->bctype == FVBC_INFLOW) {
371: switch (ctx->initial) {
372: case 0:
373: case 1:
374: case 2:
375: case 3:
376: case 4:
377: case 5:
378: u[0] = 0;
379: u[1] = 0.0; /* Left boundary conditions */
380: u[2] = 0;
381: u[3] = 0.0; /* Right boundary conditions */
382: break;
383: case 6:
384: u[0] = 0;
385: u[1] = 0.0; /* Left boundary conditions */
386: u[2] = 0;
387: u[3] = 0.0; /* Right boundary conditions */
388: break;
389: case 7:
390: u[0] = 0;
391: u[1] = 0.0; /* Left boundary conditions */
392: u[2] = 0;
393: u[3] = 0.0; /* Right boundary conditions */
394: break;
395: case 8:
396: u[0] = 0;
397: u[1] = 1.0; /* Left boundary conditions */
398: u[2] = 0;
399: u[3] = -1.0; /* Right boundary conditions */
400: break;
401: default:
402: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
403: }
404: }
405: PetscFunctionReturn(PETSC_SUCCESS);
406: }
408: /* Selects which boundary conditions are marked as inflow and which as outflow when FVBC_INFLOW is selected. */
409: static PetscErrorCode PhysicsSetInflowType_Shallow(FVCtx *ctx)
410: {
411: PetscFunctionBeginUser;
412: switch (ctx->initial) {
413: case 0:
414: case 1:
415: case 2:
416: case 3:
417: case 4:
418: case 5:
419: case 6:
420: case 7:
421: case 8: /* Fix left and right momentum, height is outflow */
422: ctx->physics2.bcinflowindex[0] = PETSC_FALSE;
423: ctx->physics2.bcinflowindex[1] = PETSC_TRUE;
424: ctx->physics2.bcinflowindex[2] = PETSC_FALSE;
425: ctx->physics2.bcinflowindex[3] = PETSC_TRUE;
426: break;
427: default:
428: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
429: }
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: static PetscErrorCode PhysicsCreate_Shallow(FVCtx *ctx)
434: {
435: ShallowCtx *user;
436: PetscFunctionList rlist = NULL, rclist = NULL;
437: char rname[256] = "rusanov", rcname[256] = "characteristic";
439: PetscFunctionBeginUser;
440: PetscCall(PetscNew(&user));
441: ctx->physics2.sample2 = PhysicsSample_Shallow;
442: ctx->physics2.inflow = PhysicsInflow_Shallow;
443: ctx->physics2.destroy = PhysicsDestroy_SimpleFree;
444: ctx->physics2.riemann2 = PhysicsRiemann_Shallow_Rusanov;
445: ctx->physics2.characteristic2 = PhysicsCharacteristic_Shallow;
446: ctx->physics2.user = user;
447: ctx->physics2.dof = 2;
449: PetscCall(PetscMalloc1(2 * (ctx->physics2.dof), &ctx->physics2.bcinflowindex));
450: PetscCall(PhysicsSetInflowType_Shallow(ctx));
452: PetscCall(PetscStrallocpy("height", &ctx->physics2.fieldname[0]));
453: PetscCall(PetscStrallocpy("momentum", &ctx->physics2.fieldname[1]));
455: user->gravity = 9.81;
457: PetscCall(RiemannListAdd_2WaySplit(&rlist, "exact", PhysicsRiemann_Shallow_Exact));
458: PetscCall(RiemannListAdd_2WaySplit(&rlist, "rusanov", PhysicsRiemann_Shallow_Rusanov));
459: PetscCall(ReconstructListAdd_2WaySplit(&rclist, "characteristic", PhysicsCharacteristic_Shallow));
460: PetscCall(ReconstructListAdd_2WaySplit(&rclist, "conservative", PhysicsCharacteristic_Conservative));
461: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for Shallow", "");
462: PetscCall(PetscOptionsReal("-physics_shallow_gravity", "Gravity", "", user->gravity, &user->gravity, NULL));
463: PetscCall(PetscOptionsFList("-physics_shallow_riemann", "Riemann solver", "", rlist, rname, rname, sizeof(rname), NULL));
464: PetscCall(PetscOptionsFList("-physics_shallow_reconstruct", "Reconstruction", "", rclist, rcname, rcname, sizeof(rcname), NULL));
465: PetscOptionsEnd();
466: PetscCall(RiemannListFind_2WaySplit(rlist, rname, &ctx->physics2.riemann2));
467: PetscCall(ReconstructListFind_2WaySplit(rclist, rcname, &ctx->physics2.characteristic2));
468: PetscCall(PetscFunctionListDestroy(&rlist));
469: PetscCall(PetscFunctionListDestroy(&rclist));
470: PetscFunctionReturn(PETSC_SUCCESS);
471: }
473: PetscErrorCode FVSample_2WaySplit(FVCtx *ctx, DM da, PetscReal time, Vec U)
474: {
475: PetscScalar *u, *uj, xj, xi;
476: PetscInt i, j, k, dof, xs, xm, Mx;
477: const PetscInt N = 200;
478: PetscReal hs, hf;
480: PetscFunctionBeginUser;
481: PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
482: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
483: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
484: PetscCall(DMDAVecGetArray(da, U, &u));
485: PetscCall(PetscMalloc1(dof, &uj));
486: hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
487: hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
488: for (i = xs; i < xs + xm; i++) {
489: if (i < ctx->sf) {
490: xi = ctx->xmin + 0.5 * hs + i * hs;
491: /* Integrate over cell i using trapezoid rule with N points. */
492: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
493: for (j = 0; j < N + 1; j++) {
494: xj = xi + hs * (j - N / 2) / (PetscReal)N;
495: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
496: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
497: }
498: } else if (i < ctx->fs) {
499: xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf;
500: /* Integrate over cell i using trapezoid rule with N points. */
501: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
502: for (j = 0; j < N + 1; j++) {
503: xj = xi + hf * (j - N / 2) / (PetscReal)N;
504: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
505: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
506: }
507: } else {
508: xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs;
509: /* Integrate over cell i using trapezoid rule with N points. */
510: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
511: for (j = 0; j < N + 1; j++) {
512: xj = xi + hs * (j - N / 2) / (PetscReal)N;
513: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
514: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
515: }
516: }
517: }
518: PetscCall(DMDAVecRestoreArray(da, U, &u));
519: PetscCall(PetscFree(uj));
520: PetscFunctionReturn(PETSC_SUCCESS);
521: }
523: static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
524: {
525: Vec Y;
526: PetscInt i, Mx;
527: const PetscScalar *ptr_X, *ptr_Y;
528: PetscReal hs, hf;
530: PetscFunctionBeginUser;
531: PetscCall(VecGetSize(X, &Mx));
532: PetscCall(VecDuplicate(X, &Y));
533: PetscCall(FVSample_2WaySplit(ctx, da, t, Y));
534: hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
535: hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
536: PetscCall(VecGetArrayRead(X, &ptr_X));
537: PetscCall(VecGetArrayRead(Y, &ptr_Y));
538: for (i = 0; i < Mx; i++) {
539: if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
540: else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
541: }
542: PetscCall(VecRestoreArrayRead(X, &ptr_X));
543: PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
544: PetscCall(VecDestroy(&Y));
545: PetscFunctionReturn(PETSC_SUCCESS);
546: }
548: PetscErrorCode FVRHSFunction_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
549: {
550: FVCtx *ctx = (FVCtx *)vctx;
551: PetscInt i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
552: PetscReal hxf, hxs;
553: PetscScalar *x, *f, *slope;
554: Vec Xloc;
555: DM da;
557: PetscFunctionBeginUser;
558: ctx->cfl_idt = 0;
559: PetscCall(TSGetDM(ts, &da));
560: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
561: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); /* Mx is the number of center points */
562: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
563: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
564: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
565: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
567: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
569: PetscCall(DMDAVecGetArray(da, Xloc, &x));
570: PetscCall(DMDAVecGetArray(da, F, &f));
571: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
572: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
574: if (ctx->bctype == FVBC_OUTFLOW) {
575: for (i = xs - 2; i < 0; i++) {
576: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
577: }
578: for (i = Mx; i < xs + xm + 2; i++) {
579: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
580: }
581: }
583: if (ctx->bctype == FVBC_INFLOW) {
584: /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
585: pages 137-138 for the scheme. */
586: if (xs == 0) { /* Left Boundary */
587: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
588: for (j = 0; j < dof; j++) {
589: if (ctx->physics2.bcinflowindex[j]) {
590: for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
591: } else {
592: for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
593: }
594: }
595: }
596: if (xs + xm == Mx) { /* Right Boundary */
597: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
598: for (j = 0; j < dof; j++) {
599: if (ctx->physics2.bcinflowindex[dof + j]) {
600: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
601: } else {
602: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
603: }
604: }
605: }
606: }
608: for (i = xs - 1; i < xs + xm + 1; i++) {
609: struct _LimitInfo info;
610: PetscScalar *cjmpL, *cjmpR;
611: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
612: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
613: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
614: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
615: cjmpL = &ctx->cjmpLR[0];
616: cjmpR = &ctx->cjmpLR[dof];
617: for (j = 0; j < dof; j++) {
618: PetscScalar jmpL, jmpR;
619: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
620: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
621: for (k = 0; k < dof; k++) {
622: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
623: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
624: }
625: }
626: /* Apply limiter to the left and right characteristic jumps */
627: info.m = dof;
628: info.hxs = hxs;
629: info.hxf = hxf;
630: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
631: for (j = 0; j < dof; j++) {
632: PetscScalar tmp = 0;
633: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
634: slope[i * dof + j] = tmp;
635: }
636: }
638: for (i = xs; i < xs + xm + 1; i++) {
639: PetscReal maxspeed;
640: PetscScalar *uL, *uR;
641: uL = &ctx->uLR[0];
642: uR = &ctx->uLR[dof];
643: if (i < sf) { /* slow region */
644: for (j = 0; j < dof; j++) {
645: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
646: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
647: }
648: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
649: if (i > xs) {
650: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
651: }
652: if (i < xs + xm) {
653: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
654: }
655: } else if (i == sf) { /* interface between the slow region and the fast region */
656: for (j = 0; j < dof; j++) {
657: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
658: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
659: }
660: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
661: if (i > xs) {
662: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
663: }
664: if (i < xs + xm) {
665: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
666: }
667: } else if (i > sf && i < fs) { /* fast region */
668: for (j = 0; j < dof; j++) {
669: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
670: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
671: }
672: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
673: if (i > xs) {
674: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
675: }
676: if (i < xs + xm) {
677: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
678: }
679: } else if (i == fs) { /* interface between the fast region and the slow region */
680: for (j = 0; j < dof; j++) {
681: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
682: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
683: }
684: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
685: if (i > xs) {
686: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
687: }
688: if (i < xs + xm) {
689: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
690: }
691: } else { /* slow region */
692: for (j = 0; j < dof; j++) {
693: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
694: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
695: }
696: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
697: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
698: if (i > xs) {
699: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
700: }
701: if (i < xs + xm) {
702: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
703: }
704: }
705: }
706: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
707: PetscCall(DMDAVecRestoreArray(da, F, &f));
708: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
709: PetscCall(DMRestoreLocalVector(da, &Xloc));
710: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
711: if (0) {
712: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
713: PetscReal dt, tnow;
714: PetscCall(TSGetTimeStep(ts, &dt));
715: PetscCall(TSGetTime(ts, &tnow));
716: if (dt > 0.5 / ctx->cfl_idt) {
717: if (1) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
718: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_PLIB, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
719: }
720: }
721: PetscFunctionReturn(PETSC_SUCCESS);
722: }
724: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
725: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
726: {
727: FVCtx *ctx = (FVCtx *)vctx;
728: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
729: PetscReal hxs, hxf;
730: PetscScalar *x, *f, *slope;
731: Vec Xloc;
732: DM da;
734: PetscFunctionBeginUser;
735: ctx->cfl_idt = 0;
736: PetscCall(TSGetDM(ts, &da));
737: PetscCall(DMGetLocalVector(da, &Xloc));
738: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
739: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
740: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
741: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
742: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
743: PetscCall(VecZeroEntries(F));
744: PetscCall(DMDAVecGetArray(da, Xloc, &x));
745: PetscCall(VecGetArray(F, &f));
746: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
747: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
749: if (ctx->bctype == FVBC_OUTFLOW) {
750: for (i = xs - 2; i < 0; i++) {
751: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
752: }
753: for (i = Mx; i < xs + xm + 2; i++) {
754: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
755: }
756: }
757: if (ctx->bctype == FVBC_INFLOW) {
758: /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
759: pages 137-138 for the scheme. */
760: if (xs == 0) { /* Left Boundary */
761: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
762: for (j = 0; j < dof; j++) {
763: if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
764: for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
765: } else {
766: for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
767: }
768: }
769: }
770: if (xs + xm == Mx) { /* Right Boundary */
771: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
772: for (j = 0; j < dof; j++) {
773: if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
774: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
775: } else {
776: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
777: }
778: }
779: }
780: }
781: for (i = xs - 1; i < xs + xm + 1; i++) {
782: struct _LimitInfo info;
783: PetscScalar *cjmpL, *cjmpR;
784: if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */
785: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
786: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
787: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
788: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
789: cjmpL = &ctx->cjmpLR[0];
790: cjmpR = &ctx->cjmpLR[dof];
791: for (j = 0; j < dof; j++) {
792: PetscScalar jmpL, jmpR;
793: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
794: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
795: for (k = 0; k < dof; k++) {
796: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
797: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
798: }
799: }
800: /* Apply limiter to the left and right characteristic jumps */
801: info.m = dof;
802: info.hxs = hxs;
803: info.hxf = hxf;
804: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
805: for (j = 0; j < dof; j++) {
806: PetscScalar tmp = 0;
807: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
808: slope[i * dof + j] = tmp;
809: }
810: }
811: }
813: for (i = xs; i < xs + xm + 1; i++) {
814: PetscReal maxspeed;
815: PetscScalar *uL, *uR;
816: uL = &ctx->uLR[0];
817: uR = &ctx->uLR[dof];
818: if (i < sf - lsbwidth) { /* slow region */
819: for (j = 0; j < dof; j++) {
820: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
821: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
822: }
823: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
824: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
825: if (i > xs) {
826: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
827: }
828: if (i < xs + xm) {
829: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
830: islow++;
831: }
832: }
833: if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */
834: for (j = 0; j < dof; j++) {
835: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
836: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
837: }
838: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
839: if (i > xs) {
840: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
841: }
842: }
843: if (i == fs + rsbwidth) { /* slow region */
844: for (j = 0; j < dof; j++) {
845: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
846: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
847: }
848: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
849: if (i < xs + xm) {
850: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
851: islow++;
852: }
853: }
854: if (i > fs + rsbwidth) { /* slow region */
855: for (j = 0; j < dof; j++) {
856: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
857: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
858: }
859: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
860: if (i > xs) {
861: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
862: }
863: if (i < xs + xm) {
864: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
865: islow++;
866: }
867: }
868: }
869: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
870: PetscCall(VecRestoreArray(F, &f));
871: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
872: PetscCall(DMRestoreLocalVector(da, &Xloc));
873: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
874: PetscFunctionReturn(PETSC_SUCCESS);
875: }
877: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
878: {
879: FVCtx *ctx = (FVCtx *)vctx;
880: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
881: PetscReal hxs, hxf;
882: PetscScalar *x, *f, *slope;
883: Vec Xloc;
884: DM da;
886: PetscFunctionBeginUser;
887: PetscCall(TSGetDM(ts, &da));
888: PetscCall(DMGetLocalVector(da, &Xloc));
889: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
890: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
891: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
892: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
893: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
894: PetscCall(VecZeroEntries(F));
895: PetscCall(DMDAVecGetArray(da, Xloc, &x));
896: PetscCall(VecGetArray(F, &f));
897: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
898: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
900: if (ctx->bctype == FVBC_OUTFLOW) {
901: for (i = xs - 2; i < 0; i++) {
902: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
903: }
904: for (i = Mx; i < xs + xm + 2; i++) {
905: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
906: }
907: }
908: if (ctx->bctype == FVBC_INFLOW) {
909: /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
910: pages 137-138 for the scheme. */
911: if (xs == 0) { /* Left Boundary */
912: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
913: for (j = 0; j < dof; j++) {
914: if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
915: for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
916: } else {
917: for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
918: }
919: }
920: }
921: if (xs + xm == Mx) { /* Right Boundary */
922: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
923: for (j = 0; j < dof; j++) {
924: if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
925: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
926: } else {
927: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
928: }
929: }
930: }
931: }
932: for (i = xs - 1; i < xs + xm + 1; i++) {
933: struct _LimitInfo info;
934: PetscScalar *cjmpL, *cjmpR;
935: if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) {
936: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
937: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
938: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
939: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
940: cjmpL = &ctx->cjmpLR[0];
941: cjmpR = &ctx->cjmpLR[dof];
942: for (j = 0; j < dof; j++) {
943: PetscScalar jmpL, jmpR;
944: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
945: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
946: for (k = 0; k < dof; k++) {
947: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
948: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
949: }
950: }
951: /* Apply limiter to the left and right characteristic jumps */
952: info.m = dof;
953: info.hxs = hxs;
954: info.hxf = hxf;
955: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
956: for (j = 0; j < dof; j++) {
957: PetscScalar tmp = 0;
958: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
959: slope[i * dof + j] = tmp;
960: }
961: }
962: }
964: for (i = xs; i < xs + xm + 1; i++) {
965: PetscReal maxspeed;
966: PetscScalar *uL, *uR;
967: uL = &ctx->uLR[0];
968: uR = &ctx->uLR[dof];
969: if (i == sf - lsbwidth) {
970: for (j = 0; j < dof; j++) {
971: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
972: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
973: }
974: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
975: if (i < xs + xm) {
976: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
977: islow++;
978: }
979: }
980: if (i > sf - lsbwidth && i < sf) {
981: for (j = 0; j < dof; j++) {
982: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
983: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
984: }
985: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
986: if (i > xs) {
987: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
988: }
989: if (i < xs + xm) {
990: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
991: islow++;
992: }
993: }
994: if (i == sf) { /* interface between the slow region and the fast region */
995: for (j = 0; j < dof; j++) {
996: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
997: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
998: }
999: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1000: if (i > xs) {
1001: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
1002: }
1003: }
1004: if (i == fs) { /* interface between the fast region and the slow region */
1005: for (j = 0; j < dof; j++) {
1006: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1007: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
1008: }
1009: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1010: if (i < xs + xm) {
1011: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
1012: islow++;
1013: }
1014: }
1015: if (i > fs && i < fs + rsbwidth) {
1016: for (j = 0; j < dof; j++) {
1017: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
1018: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
1019: }
1020: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1021: if (i > xs) {
1022: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
1023: }
1024: if (i < xs + xm) {
1025: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
1026: islow++;
1027: }
1028: }
1029: if (i == fs + rsbwidth) {
1030: for (j = 0; j < dof; j++) {
1031: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
1032: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
1033: }
1034: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1035: if (i > xs) {
1036: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
1037: }
1038: }
1039: }
1040: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
1041: PetscCall(VecRestoreArray(F, &f));
1042: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
1043: PetscCall(DMRestoreLocalVector(da, &Xloc));
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
1048: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
1049: {
1050: FVCtx *ctx = (FVCtx *)vctx;
1051: PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
1052: PetscReal hxs, hxf;
1053: PetscScalar *x, *f, *slope;
1054: Vec Xloc;
1055: DM da;
1057: PetscFunctionBeginUser;
1058: PetscCall(TSGetDM(ts, &da));
1059: PetscCall(DMGetLocalVector(da, &Xloc));
1060: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1061: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
1062: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
1063: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
1064: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
1065: PetscCall(VecZeroEntries(F));
1066: PetscCall(DMDAVecGetArray(da, Xloc, &x));
1067: PetscCall(VecGetArray(F, &f));
1068: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
1069: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1071: if (ctx->bctype == FVBC_OUTFLOW) {
1072: for (i = xs - 2; i < 0; i++) {
1073: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
1074: }
1075: for (i = Mx; i < xs + xm + 2; i++) {
1076: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
1077: }
1078: }
1079: if (ctx->bctype == FVBC_INFLOW) {
1080: /* See LeVeque, R. (2002). Finite Volume Methods for Hyperbolic Problems. doi:10.1017/CBO9780511791253
1081: pages 137-138 for the scheme. */
1082: if (xs == 0) { /* Left Boundary */
1083: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmin, ctx->ub));
1084: for (j = 0; j < dof; j++) {
1085: if (ctx->physics2.bcinflowindex[j] == PETSC_TRUE) {
1086: for (i = -2; i < 0; i++) x[i * dof + j] = 2.0 * ctx->ub[j] - x[-(i + 1) * dof + j];
1087: } else {
1088: for (i = -2; i < 0; i++) x[i * dof + j] = x[j]; /* Outflow */
1089: }
1090: }
1091: }
1092: if (xs + xm == Mx) { /* Right Boundary */
1093: PetscCall(ctx->physics2.inflow(ctx, time, ctx->xmax, ctx->ub));
1094: for (j = 0; j < dof; j++) {
1095: if (ctx->physics2.bcinflowindex[dof + j] == PETSC_TRUE) {
1096: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = 2.0 * ctx->ub[dof + j] - x[(2 * Mx - (i + 1)) * dof + j];
1097: } else {
1098: for (i = Mx; i < Mx + 2; i++) x[i * dof + j] = x[(Mx - 1) * dof + j]; /* Outflow */
1099: }
1100: }
1101: }
1102: }
1103: for (i = xs - 1; i < xs + xm + 1; i++) {
1104: struct _LimitInfo info;
1105: PetscScalar *cjmpL, *cjmpR;
1106: if (i > sf - 2 && i < fs + 1) {
1107: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
1108: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
1109: cjmpL = &ctx->cjmpLR[0];
1110: cjmpR = &ctx->cjmpLR[dof];
1111: for (j = 0; j < dof; j++) {
1112: PetscScalar jmpL, jmpR;
1113: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
1114: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
1115: for (k = 0; k < dof; k++) {
1116: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
1117: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
1118: }
1119: }
1120: /* Apply limiter to the left and right characteristic jumps */
1121: info.m = dof;
1122: info.hxs = hxs;
1123: info.hxf = hxf;
1124: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
1125: for (j = 0; j < dof; j++) {
1126: PetscScalar tmp = 0;
1127: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
1128: slope[i * dof + j] = tmp;
1129: }
1130: }
1131: }
1133: for (i = xs; i < xs + xm + 1; i++) {
1134: PetscReal maxspeed;
1135: PetscScalar *uL, *uR;
1136: uL = &ctx->uLR[0];
1137: uR = &ctx->uLR[dof];
1138: if (i == sf) { /* interface between the slow region and the fast region */
1139: for (j = 0; j < dof; j++) {
1140: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
1141: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1142: }
1143: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1144: if (i < xs + xm) {
1145: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1146: ifast++;
1147: }
1148: }
1149: if (i > sf && i < fs) { /* fast region */
1150: for (j = 0; j < dof; j++) {
1151: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1152: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1153: }
1154: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1155: if (i > xs) {
1156: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1157: }
1158: if (i < xs + xm) {
1159: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1160: ifast++;
1161: }
1162: }
1163: if (i == fs) { /* interface between the fast region and the slow region */
1164: for (j = 0; j < dof; j++) {
1165: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1166: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
1167: }
1168: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1169: if (i > xs) {
1170: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1171: }
1172: }
1173: }
1174: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
1175: PetscCall(VecRestoreArray(F, &f));
1176: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
1177: PetscCall(DMRestoreLocalVector(da, &Xloc));
1178: PetscFunctionReturn(PETSC_SUCCESS);
1179: }
1181: int main(int argc, char *argv[])
1182: {
1183: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
1184: PetscFunctionList limiters = NULL, physics = NULL;
1185: MPI_Comm comm;
1186: TS ts;
1187: DM da;
1188: Vec X, X0, R;
1189: FVCtx ctx;
1190: PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, islowbuffer = 0, *index_slow, *index_fast, *index_slowbuffer;
1191: PetscBool view_final = PETSC_FALSE;
1192: PetscReal ptime, maxtime;
1194: PetscFunctionBeginUser;
1195: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1196: comm = PETSC_COMM_WORLD;
1197: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
1199: /* Register limiters to be available on the command line */
1200: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind));
1201: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff));
1202: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming));
1203: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm));
1204: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod));
1205: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee));
1206: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC));
1207: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3));
1209: /* Register physical models to be available on the command line */
1210: PetscCall(PetscFunctionListAdd(&physics, "shallow", PhysicsCreate_Shallow));
1211: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
1213: ctx.comm = comm;
1214: ctx.cfl = 0.9;
1215: ctx.bctype = FVBC_PERIODIC;
1216: ctx.xmin = -1.0;
1217: ctx.xmax = 1.0;
1218: ctx.initial = 1;
1219: ctx.hratio = 2;
1220: maxtime = 10.0;
1221: ctx.simulation = PETSC_FALSE;
1222: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
1223: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
1224: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
1225: PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
1226: PetscCall(PetscOptionsFList("-physics", "Name of physics model to use", "", physics, physname, physname, sizeof(physname), NULL));
1227: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
1228: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
1229: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
1230: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
1231: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
1232: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
1233: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
1234: PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
1235: PetscOptionsEnd();
1237: /* Choose the limiter from the list of registered limiters */
1238: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit2));
1239: PetscCheck(ctx.limit2, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
1241: /* Choose the physics from the list of registered models */
1242: {
1243: PetscErrorCode (*r)(FVCtx *);
1244: PetscCall(PetscFunctionListFind(physics, physname, &r));
1245: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
1246: /* Create the physics, will set the number of fields and their names */
1247: PetscCall((*r)(&ctx));
1248: }
1250: /* Create a DMDA to manage the parallel grid */
1251: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da));
1252: PetscCall(DMSetFromOptions(da));
1253: PetscCall(DMSetUp(da));
1254: /* Inform the DMDA of the field names provided by the physics. */
1255: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1256: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i]));
1257: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1258: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1260: /* Set coordinates of cell centers */
1261: 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));
1263: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1264: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
1265: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
1266: PetscCall(PetscMalloc1(2 * dof, &ctx.ub));
1268: /* Create a vector to store the solution and to save the initial state */
1269: PetscCall(DMCreateGlobalVector(da, &X));
1270: PetscCall(VecDuplicate(X, &X0));
1271: PetscCall(VecDuplicate(X, &R));
1273: /* create index for slow parts and fast parts,
1274: count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
1275: count_slow = Mx * 3 / (3 + ctx.hratio); // compute Mx / (1.0 + ctx.hratio / 3.0);
1276: PetscCheck(count_slow % 2 == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hratio/3) is even");
1277: count_fast = Mx - count_slow;
1278: ctx.sf = count_slow / 2;
1279: ctx.fs = ctx.sf + count_fast;
1280: PetscCall(PetscMalloc1(xm * dof, &index_slow));
1281: PetscCall(PetscMalloc1(xm * dof, &index_fast));
1282: PetscCall(PetscMalloc1(8 * dof, &index_slowbuffer));
1283: ctx.lsbwidth = 4;
1284: ctx.rsbwidth = 4;
1286: for (i = xs; i < xs + xm; i++) {
1287: if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1)
1288: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
1289: else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1))
1290: for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
1291: else
1292: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
1293: }
1294: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
1295: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
1296: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));
1298: /* Create a time-stepping object */
1299: PetscCall(TSCreate(comm, &ts));
1300: PetscCall(TSSetDM(ts, da));
1301: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx));
1302: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
1303: PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
1304: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
1305: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx));
1306: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx));
1307: PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx));
1309: PetscCall(TSSetType(ts, TSMPRK));
1310: PetscCall(TSSetMaxTime(ts, maxtime));
1311: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1313: /* Compute initial conditions and starting time step */
1314: PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0));
1315: PetscCall(FVRHSFunction_2WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
1316: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
1317: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
1318: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
1319: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1320: {
1321: PetscInt steps;
1322: PetscScalar mass_initial, mass_final, mass_difference;
1323: const PetscScalar *ptr_X, *ptr_X0;
1324: const PetscReal hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow;
1325: const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;
1327: PetscCall(TSSolve(ts, X));
1328: PetscCall(TSGetSolveTime(ts, &ptime));
1329: PetscCall(TSGetStepNumber(ts, &steps));
1330: /* calculate the total mass at initial time and final time */
1331: mass_initial = 0.0;
1332: mass_final = 0.0;
1333: PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
1334: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
1335: for (i = xs; i < xs + xm; i++) {
1336: if (i < ctx.sf || i > ctx.fs - 1) {
1337: for (k = 0; k < dof; k++) {
1338: mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
1339: mass_final = mass_final + hs * ptr_X[i * dof + k];
1340: }
1341: } else {
1342: for (k = 0; k < dof; k++) {
1343: mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
1344: mass_final = mass_final + hf * ptr_X[i * dof + k];
1345: }
1346: }
1347: }
1348: PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
1349: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
1350: mass_difference = mass_final - mass_initial;
1351: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mass_difference, 1, MPIU_SCALAR, MPIU_SUM, comm));
1352: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
1353: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
1354: PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1.0 / ctx.cfl_idt)));
1355: if (ctx.exact) {
1356: PetscReal nrm1 = 0;
1357: PetscCall(SolutionErrorNorms_2WaySplit(&ctx, da, ptime, X, &nrm1));
1358: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
1359: }
1360: if (ctx.simulation) {
1361: PetscReal nrm1 = 0;
1362: PetscViewer fd;
1363: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
1364: Vec XR;
1365: PetscBool flg;
1366: const PetscScalar *ptr_XR;
1367: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
1368: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
1369: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
1370: PetscCall(VecDuplicate(X0, &XR));
1371: PetscCall(VecLoad(XR, fd));
1372: PetscCall(PetscViewerDestroy(&fd));
1373: PetscCall(VecGetArrayRead(X, &ptr_X));
1374: PetscCall(VecGetArrayRead(XR, &ptr_XR));
1375: for (i = xs; i < xs + xm; i++) {
1376: if (i < ctx.sf || i > ctx.fs - 1)
1377: for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1378: else
1379: for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1380: }
1381: PetscCall(VecRestoreArrayRead(X, &ptr_X));
1382: PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
1383: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
1384: PetscCall(VecDestroy(&XR));
1385: }
1386: }
1388: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1389: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
1390: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
1391: if (draw & 0x4) {
1392: Vec Y;
1393: PetscCall(VecDuplicate(X, &Y));
1394: PetscCall(FVSample_2WaySplit(&ctx, da, ptime, Y));
1395: PetscCall(VecAYPX(Y, -1, X));
1396: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
1397: PetscCall(VecDestroy(&Y));
1398: }
1400: if (view_final) {
1401: PetscViewer viewer;
1402: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
1403: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
1404: PetscCall(VecView(X, viewer));
1405: PetscCall(PetscViewerPopFormat(viewer));
1406: PetscCall(PetscViewerDestroy(&viewer));
1407: }
1409: /* Clean up */
1410: PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
1411: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
1412: PetscCall(PetscFree(ctx.physics2.bcinflowindex));
1413: PetscCall(PetscFree(ctx.ub));
1414: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
1415: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
1416: PetscCall(VecDestroy(&X));
1417: PetscCall(VecDestroy(&X0));
1418: PetscCall(VecDestroy(&R));
1419: PetscCall(DMDestroy(&da));
1420: PetscCall(TSDestroy(&ts));
1421: PetscCall(ISDestroy(&ctx.iss));
1422: PetscCall(ISDestroy(&ctx.isf));
1423: PetscCall(ISDestroy(&ctx.issb));
1424: PetscCall(PetscFree(index_slow));
1425: PetscCall(PetscFree(index_fast));
1426: PetscCall(PetscFree(index_slowbuffer));
1427: PetscCall(PetscFunctionListDestroy(&limiters));
1428: PetscCall(PetscFunctionListDestroy(&physics));
1429: PetscCall(PetscFinalize());
1430: return 0;
1431: }
1433: /*TEST
1435: build:
1436: requires: !complex !single
1437: depends: finitevolume1d.c
1439: test:
1440: suffix: 1
1441: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22
1442: output_file: output/ex4_1.out
1444: test:
1445: suffix: 2
1446: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
1447: output_file: output/ex4_1.out
1449: test:
1450: suffix: 3
1451: args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_time_step 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 0
1452: output_file: output/ex4_3.out
1454: test:
1455: suffix: 4
1456: nsize: 2
1457: args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_time_step 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1
1458: output_file: output/ex4_3.out
1460: test:
1461: suffix: 5
1462: nsize: 4
1463: args: -da_grid_x 40 -initial 1 -hratio 2 -limit mc -ts_time_step 0.1 -ts_max_steps 24 -ts_max_time 7.0 -ts_type mprk -ts_mprk_type 2a22 -physics shallow -bc_type outflow -xmin 0 -xmax 50 -ts_use_splitrhsfunction 1
1464: output_file: output/ex4_3.out
1465: TEST*/