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