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