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_dt 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_dt 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_dt 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_dt 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_dt 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, cfl_idt = 0;
554: PetscScalar *x, *f, *slope;
555: Vec Xloc;
556: DM da;
558: PetscFunctionBeginUser;
559: PetscCall(TSGetDM(ts, &da));
560: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
561: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* 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, 0, 0, &xm, 0, 0));
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: cfl_idt = PetscMax(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: PetscCall(MPIU_Allreduce(&cfl_idt, &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) {
718: 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, cfl_idt = 0;
731: PetscScalar *x, *f, *slope;
732: Vec Xloc;
733: DM da;
735: PetscFunctionBeginUser;
736: PetscCall(TSGetDM(ts, &da));
737: PetscCall(DMGetLocalVector(da, &Xloc));
738: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
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, 0, 0, &xm, 0, 0));
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: cfl_idt = PetscMax(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: PetscCall(MPIU_Allreduce(&cfl_idt, &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, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
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, 0, 0, &xm, 0, 0));
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, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
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, 0, 0, &xm, 0, 0));
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 = 0, physics = 0;
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, 0, 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, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
1258: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
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, mass_differenceg;
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: PetscCall(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
1352: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
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_dt 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_dt 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_dt 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_dt 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_dt 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*/