Actual source code: ex8.c
1: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
2: " advection - Constant coefficient scalar advection\n"
3: " u_t + (a*u)_x = 0\n"
4: " for this toy problem, we choose different meshsizes for different sub-domains (slow-medium-fast-medium-slow), \n"
5: " the meshsize ratio between two adjacient sub-domains is controlled with -hratio,\n"
6: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
7: " the states across shocks and rarefactions\n"
8: " simulation - use reference solution which is generated by smaller time step size to be true solution,\n"
9: " also the reference solution should be generated by user and stored in a binary file.\n"
10: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
11: "Several initial conditions can be chosen with -initial N\n\n"
12: "The problem size should be set with -da_grid_x M\n\n";
14: #include <petscts.h>
15: #include <petscdm.h>
16: #include <petscdmda.h>
17: #include "finitevolume1d.h"
19: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
20: {
21: PetscReal range = xmax - xmin;
22: return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
23: }
25: /* --------------------------------- Advection ----------------------------------- */
26: typedef struct {
27: PetscReal a; /* advective velocity */
28: } AdvectCtx;
30: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
31: {
32: AdvectCtx *ctx = (AdvectCtx *)vctx;
33: PetscReal speed;
35: PetscFunctionBeginUser;
36: speed = ctx->a;
37: flux[0] = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
38: *maxspeed = speed;
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
43: {
44: AdvectCtx *ctx = (AdvectCtx *)vctx;
46: PetscFunctionBeginUser;
47: X[0] = 1.;
48: Xi[0] = 1.;
49: speeds[0] = ctx->a;
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
54: {
55: AdvectCtx *ctx = (AdvectCtx *)vctx;
56: PetscReal a = ctx->a, x0;
58: PetscFunctionBeginUser;
59: switch (bctype) {
60: case FVBC_OUTFLOW:
61: x0 = x - a * t;
62: break;
63: case FVBC_PERIODIC:
64: x0 = RangeMod(x - a * t, xmin, xmax);
65: break;
66: default:
67: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
68: }
69: switch (initial) {
70: case 0:
71: u[0] = (x0 < 0) ? 1 : -1;
72: break;
73: case 1:
74: u[0] = (x0 < 0) ? -1 : 1;
75: break;
76: case 2:
77: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
78: break;
79: case 3:
80: u[0] = PetscSinReal(2 * PETSC_PI * x0);
81: break;
82: case 4:
83: u[0] = PetscAbs(x0);
84: break;
85: case 5:
86: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
87: break;
88: case 6:
89: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
90: break;
91: case 7:
92: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
93: break;
94: default:
95: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
96: }
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
101: {
102: AdvectCtx *user;
104: PetscFunctionBeginUser;
105: PetscCall(PetscNew(&user));
106: ctx->physics2.sample2 = PhysicsSample_Advect;
107: ctx->physics2.riemann2 = PhysicsRiemann_Advect;
108: ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
109: ctx->physics2.destroy = PhysicsDestroy_SimpleFree;
110: ctx->physics2.user = user;
111: ctx->physics2.dof = 1;
112: PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0]));
113: user->a = 1;
114: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
115: {
116: PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
117: }
118: PetscOptionsEnd();
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: PetscErrorCode FVSample_3WaySplit(FVCtx *ctx, DM da, PetscReal time, Vec U)
123: {
124: PetscScalar *u, *uj, xj, xi;
125: PetscInt i, j, k, dof, xs, xm, Mx;
126: const PetscInt N = 200;
127: PetscReal hs, hm, hf;
129: PetscFunctionBeginUser;
130: PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
131: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
132: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
133: PetscCall(DMDAVecGetArray(da, U, &u));
134: PetscCall(PetscMalloc1(dof, &uj));
136: hs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
137: hm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
138: hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
139: for (i = xs; i < xs + xm; i++) {
140: if (i < ctx->sm) {
141: xi = ctx->xmin + 0.5 * hs + i * hs;
142: /* Integrate over cell i using trapezoid rule with N points. */
143: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
144: for (j = 0; j < N + 1; j++) {
145: xj = xi + hs * (j - N / 2) / (PetscReal)N;
146: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
147: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
148: }
149: } else if (i < ctx->mf) {
150: xi = ctx->xmin + ctx->sm * hs + 0.5 * hm + (i - ctx->sm) * hm;
151: /* Integrate over cell i using trapezoid rule with N points. */
152: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
153: for (j = 0; j < N + 1; j++) {
154: xj = xi + hm * (j - N / 2) / (PetscReal)N;
155: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
156: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
157: }
158: } else if (i < ctx->fm) {
159: xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + 0.5 * hf + (i - ctx->mf) * hf;
160: /* Integrate over cell i using trapezoid rule with N points. */
161: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
162: for (j = 0; j < N + 1; j++) {
163: xj = xi + hf * (j - N / 2) / (PetscReal)N;
164: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
165: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
166: }
167: } else if (i < ctx->ms) {
168: xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + (ctx->fm - ctx->mf) * hf + 0.5 * hm + (i - ctx->fm) * hm;
169: /* Integrate over cell i using trapezoid rule with N points. */
170: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
171: for (j = 0; j < N + 1; j++) {
172: xj = xi + hm * (j - N / 2) / (PetscReal)N;
173: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
174: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
175: }
176: } else {
177: xi = ctx->xmin + ctx->sm * hs + (ctx->mf - ctx->sm) * hm + (ctx->fm - ctx->mf) * hf + (ctx->ms - ctx->fm) * hm + 0.5 * hs + (i - ctx->ms) * hs;
178: /* Integrate over cell i using trapezoid rule with N points. */
179: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
180: for (j = 0; j < N + 1; j++) {
181: xj = xi + hs * (j - N / 2) / (PetscReal)N;
182: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
183: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
184: }
185: }
186: }
187: PetscCall(DMDAVecRestoreArray(da, U, &u));
188: PetscCall(PetscFree(uj));
189: PetscFunctionReturn(PETSC_SUCCESS);
190: }
192: static PetscErrorCode SolutionErrorNorms_3WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
193: {
194: Vec Y;
195: PetscInt i, Mx;
196: const PetscScalar *ptr_X, *ptr_Y;
197: PetscReal hs, hm, hf;
199: PetscFunctionBeginUser;
200: PetscCall(VecGetSize(X, &Mx));
201: hs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
202: hm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
203: hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
204: PetscCall(VecDuplicate(X, &Y));
205: PetscCall(FVSample_3WaySplit(ctx, da, t, Y));
206: PetscCall(VecGetArrayRead(X, &ptr_X));
207: PetscCall(VecGetArrayRead(Y, &ptr_Y));
208: for (i = 0; i < Mx; i++) {
209: if (i < ctx->sm || i > ctx->ms - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
210: else if (i < ctx->mf || i > ctx->fm - 1) *nrm1 += hm * PetscAbs(ptr_X[i] - ptr_Y[i]);
211: else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
212: }
213: PetscCall(VecRestoreArrayRead(X, &ptr_X));
214: PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
215: PetscCall(VecDestroy(&Y));
216: PetscFunctionReturn(PETSC_SUCCESS);
217: }
219: PetscErrorCode FVRHSFunction_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
220: {
221: FVCtx *ctx = (FVCtx *)vctx;
222: PetscInt i, j, k, Mx, dof, xs, xm, sm = ctx->sm, mf = ctx->mf, fm = ctx->fm, ms = ctx->ms;
223: PetscReal hxf, hxm, hxs;
224: PetscScalar *x, *f, *slope;
225: Vec Xloc;
226: DM da;
228: PetscFunctionBeginUser;
229: ctx->cfl_idt = 0;
230: PetscCall(TSGetDM(ts, &da));
231: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
232: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); /* Mx is the number of center points */
233: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
234: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
235: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
236: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
237: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
239: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
241: PetscCall(DMDAVecGetArray(da, Xloc, &x));
242: PetscCall(DMDAVecGetArray(da, F, &f));
243: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
245: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
247: if (ctx->bctype == FVBC_OUTFLOW) {
248: for (i = xs - 2; i < 0; i++) {
249: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
250: }
251: for (i = Mx; i < xs + xm + 2; i++) {
252: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
253: }
254: }
255: for (i = xs - 1; i < xs + xm + 1; i++) {
256: struct _LimitInfo info;
257: PetscScalar *cjmpL, *cjmpR;
258: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
259: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
260: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
261: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
262: cjmpL = &ctx->cjmpLR[0];
263: cjmpR = &ctx->cjmpLR[dof];
264: for (j = 0; j < dof; j++) {
265: PetscScalar jmpL, jmpR;
266: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
267: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
268: for (k = 0; k < dof; k++) {
269: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
270: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
271: }
272: }
273: /* Apply limiter to the left and right characteristic jumps */
274: info.m = dof;
275: info.hxs = hxs;
276: info.hxm = hxm;
277: info.hxf = hxf;
278: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
279: for (j = 0; j < dof; j++) {
280: PetscScalar tmp = 0;
281: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
282: slope[i * dof + j] = tmp;
283: }
284: }
286: for (i = xs; i < xs + xm + 1; i++) {
287: PetscReal maxspeed;
288: PetscScalar *uL, *uR;
289: uL = &ctx->uLR[0];
290: uR = &ctx->uLR[dof];
291: if (i < sm || i > ms) { /* slow region */
292: for (j = 0; j < dof; j++) {
293: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
294: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
295: }
296: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
297: if (i > xs) {
298: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
299: }
300: if (i < xs + xm) {
301: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
302: }
303: } else if (i == sm) { /* interface between slow and medium component */
304: for (j = 0; j < dof; j++) {
305: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
306: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
307: }
308: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
309: if (i > xs) {
310: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
311: }
312: if (i < xs + xm) {
313: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
314: }
315: } else if (i == ms) { /* interface between medium and slow regions */
316: for (j = 0; j < dof; j++) {
317: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
318: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
319: }
320: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
321: if (i > xs) {
322: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
323: }
324: if (i < xs + xm) {
325: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
326: }
327: } else if (i < mf || i > fm) { /* medium region */
328: for (j = 0; j < dof; j++) {
329: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
330: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
331: }
332: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
333: if (i > xs) {
334: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
335: }
336: if (i < xs + xm) {
337: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
338: }
339: } else if (i == mf) { /* interface between medium and fast regions */
340: for (j = 0; j < dof; j++) {
341: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
342: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
343: }
344: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
345: if (i > xs) {
346: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
347: }
348: if (i < xs + xm) {
349: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
350: }
351: } else if (i == fm) { /* interface between fast and medium regions */
352: for (j = 0; j < dof; j++) {
353: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
354: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
355: }
356: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
357: if (i > xs) {
358: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
359: }
360: if (i < xs + xm) {
361: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
362: }
363: } else { /* fast region */
364: for (j = 0; j < dof; j++) {
365: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
366: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
367: }
368: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
369: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
370: if (i > xs) {
371: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
372: }
373: if (i < xs + xm) {
374: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
375: }
376: }
377: }
378: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
379: PetscCall(DMDAVecRestoreArray(da, F, &f));
380: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
381: PetscCall(DMRestoreLocalVector(da, &Xloc));
382: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
383: if (0) {
384: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
385: PetscReal dt, tnow;
386: PetscCall(TSGetTimeStep(ts, &dt));
387: PetscCall(TSGetTime(ts, &tnow));
388: if (dt > 0.5 / ctx->cfl_idt) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(1 / (2 * ctx->cfl_idt))));
389: }
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
394: PetscErrorCode FVRHSFunctionslow_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
395: {
396: FVCtx *ctx = (FVCtx *)vctx;
397: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
398: PetscReal hxs, hxm, hxf;
399: PetscScalar *x, *f, *slope;
400: Vec Xloc;
401: DM da;
403: PetscFunctionBeginUser;
404: ctx->cfl_idt = 0;
405: PetscCall(TSGetDM(ts, &da));
406: PetscCall(DMGetLocalVector(da, &Xloc));
407: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
408: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
409: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
410: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
411: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
412: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
413: PetscCall(VecZeroEntries(F));
414: PetscCall(DMDAVecGetArray(da, Xloc, &x));
415: PetscCall(VecGetArray(F, &f));
416: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
417: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
419: if (ctx->bctype == FVBC_OUTFLOW) {
420: for (i = xs - 2; i < 0; i++) {
421: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
422: }
423: for (i = Mx; i < xs + xm + 2; i++) {
424: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
425: }
426: }
427: for (i = xs - 1; i < xs + xm + 1; i++) {
428: struct _LimitInfo info;
429: PetscScalar *cjmpL, *cjmpR;
430: if (i < sm - lsbwidth + 1 || i > ms + rsbwidth - 2) { /* slow components and the first and last fast components */
431: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
432: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
433: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
434: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
435: cjmpL = &ctx->cjmpLR[0];
436: cjmpR = &ctx->cjmpLR[dof];
437: for (j = 0; j < dof; j++) {
438: PetscScalar jmpL, jmpR;
439: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
440: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
441: for (k = 0; k < dof; k++) {
442: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
443: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
444: }
445: }
446: /* Apply limiter to the left and right characteristic jumps */
447: info.m = dof;
448: info.hxs = hxs;
449: info.hxm = hxm;
450: info.hxf = hxf;
451: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
452: for (j = 0; j < dof; j++) {
453: PetscScalar tmp = 0;
454: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
455: slope[i * dof + j] = tmp;
456: }
457: }
458: }
460: for (i = xs; i < xs + xm + 1; i++) {
461: PetscReal maxspeed;
462: PetscScalar *uL, *uR;
463: uL = &ctx->uLR[0];
464: uR = &ctx->uLR[dof];
465: if (i < sm - lsbwidth) { /* slow region */
466: for (j = 0; j < dof; j++) {
467: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
468: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
469: }
470: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
471: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
472: if (i > xs) {
473: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
474: }
475: if (i < xs + xm) {
476: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
477: islow++;
478: }
479: }
480: if (i == sm - lsbwidth) { /* interface between slow and medium regions */
481: for (j = 0; j < dof; j++) {
482: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
483: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
484: }
485: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
486: if (i > xs) {
487: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
488: }
489: }
490: if (i == ms + rsbwidth) { /* interface between medium and slow regions */
491: for (j = 0; j < dof; j++) {
492: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
493: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
494: }
495: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
496: if (i < xs + xm) {
497: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
498: islow++;
499: }
500: }
501: if (i > ms + rsbwidth) { /* slow region */
502: for (j = 0; j < dof; j++) {
503: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
504: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
505: }
506: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
507: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
508: if (i > xs) {
509: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
510: }
511: if (i < xs + xm) {
512: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
513: islow++;
514: }
515: }
516: }
517: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
518: PetscCall(VecRestoreArray(F, &f));
519: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
520: PetscCall(DMRestoreLocalVector(da, &Xloc));
521: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
522: PetscFunctionReturn(PETSC_SUCCESS);
523: }
525: PetscErrorCode FVRHSFunctionslowbuffer_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
526: {
527: FVCtx *ctx = (FVCtx *)vctx;
528: PetscInt i, j, k, Mx, dof, xs, xm, islowbuffer = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
529: PetscReal hxs, hxm, hxf;
530: PetscScalar *x, *f, *slope;
531: Vec Xloc;
532: DM da;
534: PetscFunctionBeginUser;
535: PetscCall(TSGetDM(ts, &da));
536: PetscCall(DMGetLocalVector(da, &Xloc));
537: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
538: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
539: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
540: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
541: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
542: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
543: PetscCall(VecZeroEntries(F));
544: PetscCall(DMDAVecGetArray(da, Xloc, &x));
545: PetscCall(VecGetArray(F, &f));
546: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
547: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
549: if (ctx->bctype == FVBC_OUTFLOW) {
550: for (i = xs - 2; i < 0; i++) {
551: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
552: }
553: for (i = Mx; i < xs + xm + 2; i++) {
554: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
555: }
556: }
557: for (i = xs - 1; i < xs + xm + 1; i++) {
558: struct _LimitInfo info;
559: PetscScalar *cjmpL, *cjmpR;
560: if ((i > sm - lsbwidth - 2 && i < sm + 1) || (i > ms - 2 && i < ms + rsbwidth + 1)) {
561: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
562: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
563: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
564: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
565: cjmpL = &ctx->cjmpLR[0];
566: cjmpR = &ctx->cjmpLR[dof];
567: for (j = 0; j < dof; j++) {
568: PetscScalar jmpL, jmpR;
569: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
570: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
571: for (k = 0; k < dof; k++) {
572: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
573: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
574: }
575: }
576: /* Apply limiter to the left and right characteristic jumps */
577: info.m = dof;
578: info.hxs = hxs;
579: info.hxm = hxm;
580: info.hxf = hxf;
581: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
582: for (j = 0; j < dof; j++) {
583: PetscScalar tmp = 0;
584: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
585: slope[i * dof + j] = tmp;
586: }
587: }
588: }
590: for (i = xs; i < xs + xm + 1; i++) {
591: PetscReal maxspeed;
592: PetscScalar *uL, *uR;
593: uL = &ctx->uLR[0];
594: uR = &ctx->uLR[dof];
595: if (i == sm - lsbwidth) {
596: for (j = 0; j < dof; j++) {
597: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
598: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
599: }
600: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
601: if (i < xs + xm) {
602: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
603: islowbuffer++;
604: }
605: }
606: if (i > sm - lsbwidth && i < sm) {
607: for (j = 0; j < dof; j++) {
608: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
609: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
610: }
611: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
612: if (i > xs) {
613: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
614: }
615: if (i < xs + xm) {
616: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
617: islowbuffer++;
618: }
619: }
620: if (i == sm) { /* interface between the slow region and the medium region */
621: for (j = 0; j < dof; j++) {
622: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
623: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
624: }
625: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
626: if (i > xs) {
627: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
628: }
629: }
630: if (i == ms) { /* interface between the medium region and the slow region */
631: for (j = 0; j < dof; j++) {
632: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
633: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
634: }
635: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
636: if (i < xs + xm) {
637: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
638: islowbuffer++;
639: }
640: }
641: if (i > ms && i < ms + rsbwidth) {
642: for (j = 0; j < dof; j++) {
643: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
644: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
645: }
646: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
647: if (i > xs) {
648: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
649: }
650: if (i < xs + xm) {
651: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
652: islowbuffer++;
653: }
654: }
655: if (i == ms + rsbwidth) {
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] * hxs / 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[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
663: }
664: }
665: }
666: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
667: PetscCall(VecRestoreArray(F, &f));
668: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
669: PetscCall(DMRestoreLocalVector(da, &Xloc));
670: PetscFunctionReturn(PETSC_SUCCESS);
671: }
673: /* --------------------------------- Finite Volume Solver for medium components ----------------------------------- */
674: PetscErrorCode FVRHSFunctionmedium_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
675: {
676: FVCtx *ctx = (FVCtx *)vctx;
677: PetscInt i, j, k, Mx, dof, xs, xm, imedium = 0, sm = ctx->sm, mf = ctx->mf, fm = ctx->fm, ms = ctx->ms, lmbwidth = ctx->lmbwidth, rmbwidth = ctx->rmbwidth;
678: PetscReal hxs, hxm, hxf;
679: PetscScalar *x, *f, *slope;
680: Vec Xloc;
681: DM da;
683: PetscFunctionBeginUser;
684: PetscCall(TSGetDM(ts, &da));
685: PetscCall(DMGetLocalVector(da, &Xloc));
686: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
687: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
688: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
689: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
690: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
691: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
692: PetscCall(VecZeroEntries(F));
693: PetscCall(DMDAVecGetArray(da, Xloc, &x));
694: PetscCall(VecGetArray(F, &f));
695: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
696: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
698: if (ctx->bctype == FVBC_OUTFLOW) {
699: for (i = xs - 2; i < 0; i++) {
700: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
701: }
702: for (i = Mx; i < xs + xm + 2; i++) {
703: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
704: }
705: }
706: for (i = xs - 1; i < xs + xm + 1; i++) {
707: struct _LimitInfo info;
708: PetscScalar *cjmpL, *cjmpR;
709: if ((i > sm - 2 && i < mf - lmbwidth + 1) || (i > fm + rmbwidth - 2 && i < ms + 1)) { /* slow components and the first and last fast components */
710: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
711: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
712: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
713: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
714: cjmpL = &ctx->cjmpLR[0];
715: cjmpR = &ctx->cjmpLR[dof];
716: for (j = 0; j < dof; j++) {
717: PetscScalar jmpL, jmpR;
718: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
719: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
720: for (k = 0; k < dof; k++) {
721: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
722: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
723: }
724: }
725: /* Apply limiter to the left and right characteristic jumps */
726: info.m = dof;
727: info.hxs = hxs;
728: info.hxm = hxm;
729: info.hxf = hxf;
730: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
731: for (j = 0; j < dof; j++) {
732: PetscScalar tmp = 0;
733: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
734: slope[i * dof + j] = tmp;
735: }
736: }
737: }
739: for (i = xs; i < xs + xm + 1; i++) {
740: PetscReal maxspeed;
741: PetscScalar *uL, *uR;
742: uL = &ctx->uLR[0];
743: uR = &ctx->uLR[dof];
744: if (i == sm) { /* interface between slow and medium regions */
745: for (j = 0; j < dof; j++) {
746: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
747: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
748: }
749: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
750: if (i < xs + xm) {
751: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
752: imedium++;
753: }
754: }
755: if (i > sm && i < mf - lmbwidth) { /* medium region */
756: for (j = 0; j < dof; j++) {
757: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
758: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
759: }
760: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
761: if (i > xs) {
762: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
763: }
764: if (i < xs + xm) {
765: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
766: imedium++;
767: }
768: }
769: if (i == mf - lmbwidth) { /* interface between medium and fast regions */
770: for (j = 0; j < dof; j++) {
771: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
772: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
773: }
774: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
775: if (i > xs) {
776: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
777: }
778: }
779: if (i == fm + rmbwidth) { /* interface between fast and medium regions */
780: for (j = 0; j < dof; j++) {
781: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
782: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
783: }
784: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
785: if (i < xs + xm) {
786: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
787: imedium++;
788: }
789: }
790: if (i > fm + rmbwidth && i < ms) { /* medium region */
791: for (j = 0; j < dof; j++) {
792: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
793: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
794: }
795: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
796: if (i > xs) {
797: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
798: }
799: if (i < xs + xm) {
800: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
801: imedium++;
802: }
803: }
804: if (i == ms) { /* interface between medium and slow regions */
805: for (j = 0; j < dof; j++) {
806: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
807: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
808: }
809: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
810: if (i > xs) {
811: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
812: }
813: }
814: }
815: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
816: PetscCall(VecRestoreArray(F, &f));
817: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
818: PetscCall(DMRestoreLocalVector(da, &Xloc));
819: PetscFunctionReturn(PETSC_SUCCESS);
820: }
822: PetscErrorCode FVRHSFunctionmediumbuffer_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
823: {
824: FVCtx *ctx = (FVCtx *)vctx;
825: PetscInt i, j, k, Mx, dof, xs, xm, imediumbuffer = 0, mf = ctx->mf, fm = ctx->fm, lmbwidth = ctx->lmbwidth, rmbwidth = ctx->rmbwidth;
826: PetscReal hxs, hxm, hxf;
827: PetscScalar *x, *f, *slope;
828: Vec Xloc;
829: DM da;
831: PetscFunctionBeginUser;
832: PetscCall(TSGetDM(ts, &da));
833: PetscCall(DMGetLocalVector(da, &Xloc));
834: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
835: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
836: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
837: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
838: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
839: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
840: PetscCall(VecZeroEntries(F));
841: PetscCall(DMDAVecGetArray(da, Xloc, &x));
842: PetscCall(VecGetArray(F, &f));
843: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
844: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
846: if (ctx->bctype == FVBC_OUTFLOW) {
847: for (i = xs - 2; i < 0; i++) {
848: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
849: }
850: for (i = Mx; i < xs + xm + 2; i++) {
851: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
852: }
853: }
854: for (i = xs - 1; i < xs + xm + 1; i++) {
855: struct _LimitInfo info;
856: PetscScalar *cjmpL, *cjmpR;
857: if ((i > mf - lmbwidth - 2 && i < mf + 1) || (i > fm - 2 && i < fm + rmbwidth + 1)) {
858: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
859: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
860: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
861: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
862: cjmpL = &ctx->cjmpLR[0];
863: cjmpR = &ctx->cjmpLR[dof];
864: for (j = 0; j < dof; j++) {
865: PetscScalar jmpL, jmpR;
866: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
867: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
868: for (k = 0; k < dof; k++) {
869: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
870: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
871: }
872: }
873: /* Apply limiter to the left and right characteristic jumps */
874: info.m = dof;
875: info.hxs = hxs;
876: info.hxm = hxm;
877: info.hxf = hxf;
878: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
879: for (j = 0; j < dof; j++) {
880: PetscScalar tmp = 0;
881: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
882: slope[i * dof + j] = tmp;
883: }
884: }
885: }
887: for (i = xs; i < xs + xm + 1; i++) {
888: PetscReal maxspeed;
889: PetscScalar *uL, *uR;
890: uL = &ctx->uLR[0];
891: uR = &ctx->uLR[dof];
892: if (i == mf - lmbwidth) {
893: for (j = 0; j < dof; j++) {
894: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
895: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
896: }
897: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
898: if (i < xs + xm) {
899: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
900: imediumbuffer++;
901: }
902: }
903: if (i > mf - lmbwidth && i < mf) {
904: for (j = 0; j < dof; j++) {
905: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
906: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
907: }
908: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
909: if (i > xs) {
910: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
911: }
912: if (i < xs + xm) {
913: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
914: imediumbuffer++;
915: }
916: }
917: if (i == mf) { /* interface between the medium region and the fast region */
918: for (j = 0; j < dof; j++) {
919: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
920: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
921: }
922: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
923: if (i > xs) {
924: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
925: }
926: }
927: if (i == fm) { /* interface between the fast region and the medium region */
928: for (j = 0; j < dof; j++) {
929: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
930: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
931: }
932: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
933: if (i < xs + xm) {
934: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
935: imediumbuffer++;
936: }
937: }
938: if (i > fm && i < fm + rmbwidth) {
939: for (j = 0; j < dof; j++) {
940: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
941: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
942: }
943: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
944: if (i > xs) {
945: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
946: }
947: if (i < xs + xm) {
948: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
949: imediumbuffer++;
950: }
951: }
952: if (i == fm + rmbwidth) {
953: for (j = 0; j < dof; j++) {
954: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
955: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
956: }
957: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
958: if (i > xs) {
959: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
960: }
961: }
962: }
963: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
964: PetscCall(VecRestoreArray(F, &f));
965: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
966: PetscCall(DMRestoreLocalVector(da, &Xloc));
967: PetscFunctionReturn(PETSC_SUCCESS);
968: }
970: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
971: PetscErrorCode FVRHSFunctionfast_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
972: {
973: FVCtx *ctx = (FVCtx *)vctx;
974: PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, mf = ctx->mf, fm = ctx->fm;
975: PetscReal hxs, hxm, hxf;
976: PetscScalar *x, *f, *slope;
977: Vec Xloc;
978: DM da;
980: PetscFunctionBeginUser;
981: PetscCall(TSGetDM(ts, &da));
982: PetscCall(DMGetLocalVector(da, &Xloc));
983: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
984: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
985: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
986: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
987: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
988: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
989: PetscCall(VecZeroEntries(F));
990: PetscCall(DMDAVecGetArray(da, Xloc, &x));
991: PetscCall(VecGetArray(F, &f));
992: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
993: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
995: if (ctx->bctype == FVBC_OUTFLOW) {
996: for (i = xs - 2; i < 0; i++) {
997: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
998: }
999: for (i = Mx; i < xs + xm + 2; i++) {
1000: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
1001: }
1002: }
1003: for (i = xs - 1; i < xs + xm + 1; i++) { /* fast components and the last slow components before fast components and the first slow component after fast components */
1004: struct _LimitInfo info;
1005: PetscScalar *cjmpL, *cjmpR;
1006: if (i > mf - 2 && i < fm + 1) {
1007: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
1008: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
1009: cjmpL = &ctx->cjmpLR[0];
1010: cjmpR = &ctx->cjmpLR[dof];
1011: for (j = 0; j < dof; j++) {
1012: PetscScalar jmpL, jmpR;
1013: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
1014: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
1015: for (k = 0; k < dof; k++) {
1016: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
1017: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
1018: }
1019: }
1020: /* Apply limiter to the left and right characteristic jumps */
1021: info.m = dof;
1022: info.hxs = hxs;
1023: info.hxm = hxm;
1024: info.hxf = hxf;
1025: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
1026: for (j = 0; j < dof; j++) {
1027: PetscScalar tmp = 0;
1028: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
1029: slope[i * dof + j] = tmp;
1030: }
1031: }
1032: }
1034: for (i = xs; i < xs + xm + 1; i++) {
1035: PetscReal maxspeed;
1036: PetscScalar *uL, *uR;
1037: uL = &ctx->uLR[0];
1038: uR = &ctx->uLR[dof];
1039: if (i == mf) { /* interface between medium and fast regions */
1040: for (j = 0; j < dof; j++) {
1041: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
1042: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1043: }
1044: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1045: if (i < xs + xm) {
1046: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1047: ifast++;
1048: }
1049: }
1050: if (i > mf && i < fm) { /* fast region */
1051: for (j = 0; j < dof; j++) {
1052: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1053: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1054: }
1055: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1056: if (i > xs) {
1057: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1058: }
1059: if (i < xs + xm) {
1060: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1061: ifast++;
1062: }
1063: }
1064: if (i == fm) { /* interface between fast and medium regions */
1065: for (j = 0; j < dof; j++) {
1066: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1067: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
1068: }
1069: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1070: if (i > xs) {
1071: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1072: }
1073: }
1074: }
1075: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
1076: PetscCall(VecRestoreArray(F, &f));
1077: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
1078: PetscCall(DMRestoreLocalVector(da, &Xloc));
1079: PetscFunctionReturn(PETSC_SUCCESS);
1080: }
1082: int main(int argc, char *argv[])
1083: {
1084: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
1085: PetscFunctionList limiters = NULL, physics = NULL;
1086: MPI_Comm comm;
1087: TS ts;
1088: DM da;
1089: Vec X, X0, R;
1090: FVCtx ctx;
1091: PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_medium, count_fast, islow = 0, imedium = 0, ifast = 0, *index_slow, *index_medium, *index_fast, islowbuffer = 0, *index_slowbuffer, imediumbuffer = 0, *index_mediumbuffer;
1092: PetscBool view_final = PETSC_FALSE;
1093: PetscReal ptime;
1095: PetscFunctionBeginUser;
1096: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1097: comm = PETSC_COMM_WORLD;
1098: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
1100: /* Register limiters to be available on the command line */
1101: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit3_Upwind));
1102: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit3_LaxWendroff));
1103: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit3_BeamWarming));
1104: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit3_Fromm));
1105: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit3_Minmod));
1106: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit3_Superbee));
1107: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit3_MC));
1108: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit3_Koren3));
1110: /* Register physical models to be available on the command line */
1111: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
1113: ctx.comm = comm;
1114: ctx.cfl = 0.9;
1115: ctx.bctype = FVBC_PERIODIC;
1116: ctx.xmin = -1.0;
1117: ctx.xmax = 1.0;
1118: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
1119: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
1120: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
1121: PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
1122: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
1123: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
1124: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
1125: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
1126: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
1127: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
1128: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
1129: PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
1130: PetscOptionsEnd();
1132: /* Choose the limiter from the list of registered limiters */
1133: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit3));
1134: PetscCheck(ctx.limit3, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
1136: /* Choose the physics from the list of registered models */
1137: {
1138: PetscErrorCode (*r)(FVCtx *);
1139: PetscCall(PetscFunctionListFind(physics, physname, &r));
1140: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
1141: /* Create the physics, will set the number of fields and their names */
1142: PetscCall((*r)(&ctx));
1143: }
1145: /* Create a DMDA to manage the parallel grid */
1146: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da));
1147: PetscCall(DMSetFromOptions(da));
1148: PetscCall(DMSetUp(da));
1149: /* Inform the DMDA of the field names provided by the physics. */
1150: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1151: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i]));
1152: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
1153: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
1155: /* Set coordinates of cell centers */
1156: 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));
1158: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1159: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
1160: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
1162: /* Create a vector to store the solution and to save the initial state */
1163: PetscCall(DMCreateGlobalVector(da, &X));
1164: PetscCall(VecDuplicate(X, &X0));
1165: PetscCall(VecDuplicate(X, &R));
1167: /* create index for slow parts and fast parts,
1168: count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
1169: count_slow = Mx / (1 + ctx.hratio) / (1 + ctx.hratio);
1170: count_medium = 2 * ctx.hratio * count_slow;
1171: PetscCheck((count_slow % 2) == 0 && (count_medium % 2) == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio)^2 and Mx*2*hratio/(1+hratio)^2 is even");
1172: count_fast = ctx.hratio * ctx.hratio * count_slow;
1173: ctx.sm = count_slow / 2;
1174: ctx.mf = ctx.sm + count_medium / 2;
1175: ctx.fm = ctx.mf + count_fast;
1176: ctx.ms = ctx.fm + count_medium / 2;
1177: PetscCall(PetscMalloc1(xm * dof, &index_slow));
1178: PetscCall(PetscMalloc1(xm * dof, &index_medium));
1179: PetscCall(PetscMalloc1(xm * dof, &index_fast));
1180: PetscCall(PetscMalloc1(6 * dof, &index_slowbuffer));
1181: PetscCall(PetscMalloc1(6 * dof, &index_mediumbuffer));
1182: if (((AdvectCtx *)ctx.physics2.user)->a > 0) {
1183: ctx.lsbwidth = 2;
1184: ctx.rsbwidth = 4;
1185: ctx.lmbwidth = 2;
1186: ctx.rmbwidth = 4;
1187: } else {
1188: ctx.lsbwidth = 4;
1189: ctx.rsbwidth = 2;
1190: ctx.lmbwidth = 4;
1191: ctx.rmbwidth = 2;
1192: }
1194: for (i = xs; i < xs + xm; i++) {
1195: if (i < ctx.sm - ctx.lsbwidth || i > ctx.ms + ctx.rsbwidth - 1)
1196: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
1197: else if ((i >= ctx.sm - ctx.lsbwidth && i < ctx.sm) || (i > ctx.ms - 1 && i <= ctx.ms + ctx.rsbwidth - 1))
1198: for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
1199: else if (i < ctx.mf - ctx.lmbwidth || i > ctx.fm + ctx.rmbwidth - 1)
1200: for (k = 0; k < dof; k++) index_medium[imedium++] = i * dof + k;
1201: else if ((i >= ctx.mf - ctx.lmbwidth && i < ctx.mf) || (i > ctx.fm - 1 && i <= ctx.fm + ctx.rmbwidth - 1))
1202: for (k = 0; k < dof; k++) index_mediumbuffer[imediumbuffer++] = i * dof + k;
1203: else
1204: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
1205: }
1206: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
1207: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, imedium, index_medium, PETSC_COPY_VALUES, &ctx.ism));
1208: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
1209: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));
1210: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, imediumbuffer, index_mediumbuffer, PETSC_COPY_VALUES, &ctx.ismb));
1212: /* Create a time-stepping object */
1213: PetscCall(TSCreate(comm, &ts));
1214: PetscCall(TSSetDM(ts, da));
1215: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_3WaySplit, &ctx));
1216: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
1217: PetscCall(TSRHSSplitSetIS(ts, "medium", ctx.ism));
1218: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
1219: PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
1220: PetscCall(TSRHSSplitSetIS(ts, "mediumbuffer", ctx.ismb));
1221: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_3WaySplit, &ctx));
1222: PetscCall(TSRHSSplitSetRHSFunction(ts, "medium", NULL, FVRHSFunctionmedium_3WaySplit, &ctx));
1223: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_3WaySplit, &ctx));
1224: PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_3WaySplit, &ctx));
1225: PetscCall(TSRHSSplitSetRHSFunction(ts, "mediumbuffer", NULL, FVRHSFunctionmediumbuffer_3WaySplit, &ctx));
1227: PetscCall(TSSetType(ts, TSSSP));
1228: /*PetscCall(TSSetType(ts,TSMPRK));*/
1229: PetscCall(TSSetMaxTime(ts, 10));
1230: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1232: /* Compute initial conditions and starting time step */
1233: PetscCall(FVSample_3WaySplit(&ctx, da, 0, X0));
1234: PetscCall(FVRHSFunction_3WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
1235: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
1236: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
1237: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
1238: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1239: {
1240: PetscInt steps;
1241: PetscScalar mass_initial, mass_final, mass_difference;
1242: const PetscScalar *ptr_X, *ptr_X0;
1243: const PetscReal hs = (ctx.xmax - ctx.xmin) / 4.0 / count_slow;
1244: const PetscReal hm = (ctx.xmax - ctx.xmin) / 2.0 / count_medium;
1245: const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;
1247: PetscCall(TSSolve(ts, X));
1248: PetscCall(TSGetSolveTime(ts, &ptime));
1249: PetscCall(TSGetStepNumber(ts, &steps));
1250: /* calculate the total mass at initial time and final time */
1251: mass_initial = 0.0;
1252: mass_final = 0.0;
1253: PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
1254: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
1255: for (i = xs; i < xs + xm; i++) {
1256: if (i < ctx.sm || i > ctx.ms - 1)
1257: for (k = 0; k < dof; k++) {
1258: mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
1259: mass_final = mass_final + hs * ptr_X[i * dof + k];
1260: }
1261: else if (i < ctx.mf || i > ctx.fm - 1)
1262: for (k = 0; k < dof; k++) {
1263: mass_initial = mass_initial + hm * ptr_X0[i * dof + k];
1264: mass_final = mass_final + hm * ptr_X[i * dof + k];
1265: }
1266: else {
1267: for (k = 0; k < dof; k++) {
1268: mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
1269: mass_final = mass_final + hf * ptr_X[i * dof + k];
1270: }
1271: }
1272: }
1273: PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
1274: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
1275: mass_difference = mass_final - mass_initial;
1276: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mass_difference, 1, MPIU_SCALAR, MPIU_SUM, comm));
1277: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
1278: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
1279: PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1 / ctx.cfl_idt)));
1280: if (ctx.exact) {
1281: PetscReal nrm1 = 0;
1282: PetscCall(SolutionErrorNorms_3WaySplit(&ctx, da, ptime, X, &nrm1));
1283: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
1284: }
1285: if (ctx.simulation) {
1286: PetscReal nrm1 = 0;
1287: PetscViewer fd;
1288: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
1289: Vec XR;
1290: PetscBool flg;
1291: const PetscScalar *ptr_XR;
1292: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
1293: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
1294: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
1295: PetscCall(VecDuplicate(X0, &XR));
1296: PetscCall(VecLoad(XR, fd));
1297: PetscCall(PetscViewerDestroy(&fd));
1298: PetscCall(VecGetArrayRead(X, &ptr_X));
1299: PetscCall(VecGetArrayRead(XR, &ptr_XR));
1300: for (i = xs; i < xs + xm; i++) {
1301: if (i < ctx.sm || i > ctx.ms - 1)
1302: for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1303: else if (i < ctx.mf || i < ctx.fm - 1)
1304: for (k = 0; k < dof; k++) nrm1 = nrm1 + hm * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1305: else
1306: for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1307: }
1308: PetscCall(VecRestoreArrayRead(X, &ptr_X));
1309: PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
1310: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
1311: PetscCall(VecDestroy(&XR));
1312: }
1313: }
1315: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1316: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
1317: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
1318: if (draw & 0x4) {
1319: Vec Y;
1320: PetscCall(VecDuplicate(X, &Y));
1321: PetscCall(FVSample_3WaySplit(&ctx, da, ptime, Y));
1322: PetscCall(VecAYPX(Y, -1, X));
1323: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
1324: PetscCall(VecDestroy(&Y));
1325: }
1327: if (view_final) {
1328: PetscViewer viewer;
1329: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
1330: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
1331: PetscCall(VecView(X, viewer));
1332: PetscCall(PetscViewerPopFormat(viewer));
1333: PetscCall(PetscViewerDestroy(&viewer));
1334: }
1336: /* Clean up */
1337: PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
1338: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
1339: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
1340: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
1341: PetscCall(VecDestroy(&X));
1342: PetscCall(VecDestroy(&X0));
1343: PetscCall(VecDestroy(&R));
1344: PetscCall(DMDestroy(&da));
1345: PetscCall(TSDestroy(&ts));
1346: PetscCall(ISDestroy(&ctx.iss));
1347: PetscCall(ISDestroy(&ctx.ism));
1348: PetscCall(ISDestroy(&ctx.isf));
1349: PetscCall(ISDestroy(&ctx.issb));
1350: PetscCall(ISDestroy(&ctx.ismb));
1351: PetscCall(PetscFree(index_slow));
1352: PetscCall(PetscFree(index_medium));
1353: PetscCall(PetscFree(index_fast));
1354: PetscCall(PetscFree(index_slowbuffer));
1355: PetscCall(PetscFree(index_mediumbuffer));
1356: PetscCall(PetscFunctionListDestroy(&limiters));
1357: PetscCall(PetscFunctionListDestroy(&physics));
1358: PetscCall(PetscFinalize());
1359: return 0;
1360: }
1362: /*TEST
1364: build:
1365: requires: !complex
1366: depends: finitevolume1d.c
1368: test:
1369: suffix: 1
1370: args: -da_grid_x 90 -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 2a23 -ts_use_splitrhsfunction 0
1372: test:
1373: suffix: 2
1374: args: -da_grid_x 90 -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 2a23 -ts_use_splitrhsfunction 1
1376: TEST*/