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;
225: PetscScalar *x, *f, *slope;
226: Vec Xloc;
227: DM da;
229: PetscFunctionBeginUser;
230: ctx->cfl_idt = 0;
231: PetscCall(TSGetDM(ts, &da));
232: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
233: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */
234: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
235: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
236: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
237: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
238: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
240: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
242: PetscCall(DMDAVecGetArray(da, Xloc, &x));
243: PetscCall(DMDAVecGetArray(da, F, &f));
244: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
246: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
248: if (ctx->bctype == FVBC_OUTFLOW) {
249: for (i = xs - 2; i < 0; i++) {
250: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
251: }
252: for (i = Mx; i < xs + xm + 2; i++) {
253: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
254: }
255: }
256: for (i = xs - 1; i < xs + xm + 1; i++) {
257: struct _LimitInfo info;
258: PetscScalar *cjmpL, *cjmpR;
259: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
260: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
261: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
262: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
263: cjmpL = &ctx->cjmpLR[0];
264: cjmpR = &ctx->cjmpLR[dof];
265: for (j = 0; j < dof; j++) {
266: PetscScalar jmpL, jmpR;
267: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
268: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
269: for (k = 0; k < dof; k++) {
270: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
271: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
272: }
273: }
274: /* Apply limiter to the left and right characteristic jumps */
275: info.m = dof;
276: info.hxs = hxs;
277: info.hxm = hxm;
278: info.hxf = hxf;
279: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
280: for (j = 0; j < dof; j++) {
281: PetscScalar tmp = 0;
282: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
283: slope[i * dof + j] = tmp;
284: }
285: }
287: for (i = xs; i < xs + xm + 1; i++) {
288: PetscReal maxspeed;
289: PetscScalar *uL, *uR;
290: uL = &ctx->uLR[0];
291: uR = &ctx->uLR[dof];
292: if (i < sm || i > ms) { /* slow region */
293: for (j = 0; j < dof; j++) {
294: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
295: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
296: }
297: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
298: if (i > xs) {
299: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
300: }
301: if (i < xs + xm) {
302: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
303: }
304: } else if (i == sm) { /* interface between slow and medium component */
305: for (j = 0; j < dof; j++) {
306: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
307: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
308: }
309: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
310: if (i > xs) {
311: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
312: }
313: if (i < xs + xm) {
314: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
315: }
316: } else if (i == ms) { /* interface between medium and slow regions */
317: for (j = 0; j < dof; j++) {
318: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
319: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
320: }
321: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
322: if (i > xs) {
323: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
324: }
325: if (i < xs + xm) {
326: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
327: }
328: } else if (i < mf || i > fm) { /* medium region */
329: for (j = 0; j < dof; j++) {
330: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
331: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
332: }
333: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
334: if (i > xs) {
335: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
336: }
337: if (i < xs + xm) {
338: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
339: }
340: } else if (i == mf) { /* interface between medium and fast regions */
341: for (j = 0; j < dof; j++) {
342: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
343: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
344: }
345: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
346: if (i > xs) {
347: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxm;
348: }
349: if (i < xs + xm) {
350: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
351: }
352: } else if (i == fm) { /* interface between fast and medium regions */
353: for (j = 0; j < dof; j++) {
354: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
355: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
356: }
357: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
358: if (i > xs) {
359: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
360: }
361: if (i < xs + xm) {
362: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxm;
363: }
364: } else { /* fast region */
365: for (j = 0; j < dof; j++) {
366: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
367: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
368: }
369: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
370: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
371: if (i > xs) {
372: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
373: }
374: if (i < xs + xm) {
375: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
376: }
377: }
378: }
379: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
380: PetscCall(DMDAVecRestoreArray(da, F, &f));
381: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
382: PetscCall(DMRestoreLocalVector(da, &Xloc));
383: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
384: if (0) {
385: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
386: PetscReal dt, tnow;
387: PetscCall(TSGetTimeStep(ts, &dt));
388: PetscCall(TSGetTime(ts, &tnow));
389: 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))));
390: }
391: PetscFunctionReturn(PETSC_SUCCESS);
392: }
394: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
395: PetscErrorCode FVRHSFunctionslow_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
396: {
397: FVCtx *ctx = (FVCtx *)vctx;
398: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
399: PetscReal hxs, hxm, hxf;
400: PetscScalar *x, *f, *slope;
401: Vec Xloc;
402: DM da;
404: PetscFunctionBeginUser;
405: ctx->cfl_idt = 0;
406: PetscCall(TSGetDM(ts, &da));
407: PetscCall(DMGetLocalVector(da, &Xloc));
408: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
409: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
410: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
411: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
412: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
413: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
414: PetscCall(VecZeroEntries(F));
415: PetscCall(DMDAVecGetArray(da, Xloc, &x));
416: PetscCall(VecGetArray(F, &f));
417: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
418: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
420: if (ctx->bctype == FVBC_OUTFLOW) {
421: for (i = xs - 2; i < 0; i++) {
422: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
423: }
424: for (i = Mx; i < xs + xm + 2; i++) {
425: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
426: }
427: }
428: for (i = xs - 1; i < xs + xm + 1; i++) {
429: struct _LimitInfo info;
430: PetscScalar *cjmpL, *cjmpR;
431: if (i < sm - lsbwidth + 1 || i > ms + rsbwidth - 2) { /* slow components and the first and last fast components */
432: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
433: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
434: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
435: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
436: cjmpL = &ctx->cjmpLR[0];
437: cjmpR = &ctx->cjmpLR[dof];
438: for (j = 0; j < dof; j++) {
439: PetscScalar jmpL, jmpR;
440: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
441: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
442: for (k = 0; k < dof; k++) {
443: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
444: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
445: }
446: }
447: /* Apply limiter to the left and right characteristic jumps */
448: info.m = dof;
449: info.hxs = hxs;
450: info.hxm = hxm;
451: info.hxf = hxf;
452: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
453: for (j = 0; j < dof; j++) {
454: PetscScalar tmp = 0;
455: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
456: slope[i * dof + j] = tmp;
457: }
458: }
459: }
461: for (i = xs; i < xs + xm + 1; i++) {
462: PetscReal maxspeed;
463: PetscScalar *uL, *uR;
464: uL = &ctx->uLR[0];
465: uR = &ctx->uLR[dof];
466: if (i < sm - lsbwidth) { /* slow region */
467: for (j = 0; j < dof; j++) {
468: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
469: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
470: }
471: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
472: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
473: if (i > xs) {
474: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
475: }
476: if (i < xs + xm) {
477: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
478: islow++;
479: }
480: }
481: if (i == sm - lsbwidth) { /* interface between slow and medium regions */
482: for (j = 0; j < dof; j++) {
483: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
484: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
485: }
486: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
487: if (i > xs) {
488: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
489: }
490: }
491: if (i == ms + rsbwidth) { /* interface between medium and slow regions */
492: for (j = 0; j < dof; j++) {
493: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
494: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
495: }
496: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
497: if (i < xs + xm) {
498: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
499: islow++;
500: }
501: }
502: if (i > ms + rsbwidth) { /* slow region */
503: for (j = 0; j < dof; j++) {
504: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
505: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
506: }
507: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
508: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
509: if (i > xs) {
510: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
511: }
512: if (i < xs + xm) {
513: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
514: islow++;
515: }
516: }
517: }
518: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
519: PetscCall(VecRestoreArray(F, &f));
520: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
521: PetscCall(DMRestoreLocalVector(da, &Xloc));
522: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
523: PetscFunctionReturn(PETSC_SUCCESS);
524: }
526: PetscErrorCode FVRHSFunctionslowbuffer_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
527: {
528: FVCtx *ctx = (FVCtx *)vctx;
529: PetscInt i, j, k, Mx, dof, xs, xm, islowbuffer = 0, sm = ctx->sm, ms = ctx->ms, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
530: PetscReal hxs, hxm, hxf;
531: PetscScalar *x, *f, *slope;
532: Vec Xloc;
533: DM da;
535: PetscFunctionBeginUser;
536: PetscCall(TSGetDM(ts, &da));
537: PetscCall(DMGetLocalVector(da, &Xloc));
538: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
539: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
540: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
541: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
542: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
543: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
544: PetscCall(VecZeroEntries(F));
545: PetscCall(DMDAVecGetArray(da, Xloc, &x));
546: PetscCall(VecGetArray(F, &f));
547: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
548: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
550: if (ctx->bctype == FVBC_OUTFLOW) {
551: for (i = xs - 2; i < 0; i++) {
552: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
553: }
554: for (i = Mx; i < xs + xm + 2; i++) {
555: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
556: }
557: }
558: for (i = xs - 1; i < xs + xm + 1; i++) {
559: struct _LimitInfo info;
560: PetscScalar *cjmpL, *cjmpR;
561: if ((i > sm - lsbwidth - 2 && i < sm + 1) || (i > ms - 2 && i < ms + rsbwidth + 1)) {
562: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
563: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
564: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
565: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
566: cjmpL = &ctx->cjmpLR[0];
567: cjmpR = &ctx->cjmpLR[dof];
568: for (j = 0; j < dof; j++) {
569: PetscScalar jmpL, jmpR;
570: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
571: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
572: for (k = 0; k < dof; k++) {
573: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
574: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
575: }
576: }
577: /* Apply limiter to the left and right characteristic jumps */
578: info.m = dof;
579: info.hxs = hxs;
580: info.hxm = hxm;
581: info.hxf = hxf;
582: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
583: for (j = 0; j < dof; j++) {
584: PetscScalar tmp = 0;
585: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
586: slope[i * dof + j] = tmp;
587: }
588: }
589: }
591: for (i = xs; i < xs + xm + 1; i++) {
592: PetscReal maxspeed;
593: PetscScalar *uL, *uR;
594: uL = &ctx->uLR[0];
595: uR = &ctx->uLR[dof];
596: if (i == sm - lsbwidth) {
597: for (j = 0; j < dof; j++) {
598: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
599: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
600: }
601: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
602: if (i < xs + xm) {
603: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
604: islowbuffer++;
605: }
606: }
607: if (i > sm - lsbwidth && i < sm) {
608: for (j = 0; j < dof; j++) {
609: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
610: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
611: }
612: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
613: if (i > xs) {
614: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
615: }
616: if (i < xs + xm) {
617: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
618: islowbuffer++;
619: }
620: }
621: if (i == sm) { /* interface between the slow region and the medium region */
622: for (j = 0; j < dof; j++) {
623: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
624: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
625: }
626: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
627: if (i > xs) {
628: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
629: }
630: }
631: if (i == ms) { /* interface between the medium region and the slow region */
632: for (j = 0; j < dof; j++) {
633: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
634: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
635: }
636: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
637: if (i < xs + xm) {
638: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
639: islowbuffer++;
640: }
641: }
642: if (i > ms && i < ms + rsbwidth) {
643: for (j = 0; j < dof; j++) {
644: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
645: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
646: }
647: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
648: if (i > xs) {
649: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
650: }
651: if (i < xs + xm) {
652: for (j = 0; j < dof; j++) f[islowbuffer * dof + j] += ctx->flux[j] / hxs;
653: islowbuffer++;
654: }
655: }
656: if (i == ms + rsbwidth) {
657: for (j = 0; j < dof; j++) {
658: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
659: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
660: }
661: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
662: if (i > xs) {
663: for (j = 0; j < dof; j++) f[(islowbuffer - 1) * dof + j] -= ctx->flux[j] / hxs;
664: }
665: }
666: }
667: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
668: PetscCall(VecRestoreArray(F, &f));
669: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
670: PetscCall(DMRestoreLocalVector(da, &Xloc));
671: PetscFunctionReturn(PETSC_SUCCESS);
672: }
674: /* --------------------------------- Finite Volume Solver for medium components ----------------------------------- */
675: PetscErrorCode FVRHSFunctionmedium_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
676: {
677: FVCtx *ctx = (FVCtx *)vctx;
678: 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;
679: PetscReal hxs, hxm, hxf;
680: PetscScalar *x, *f, *slope;
681: Vec Xloc;
682: DM da;
684: PetscFunctionBeginUser;
685: PetscCall(TSGetDM(ts, &da));
686: PetscCall(DMGetLocalVector(da, &Xloc));
687: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
688: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
689: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
690: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
691: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
692: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
693: PetscCall(VecZeroEntries(F));
694: PetscCall(DMDAVecGetArray(da, Xloc, &x));
695: PetscCall(VecGetArray(F, &f));
696: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
697: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
699: if (ctx->bctype == FVBC_OUTFLOW) {
700: for (i = xs - 2; i < 0; i++) {
701: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
702: }
703: for (i = Mx; i < xs + xm + 2; i++) {
704: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
705: }
706: }
707: for (i = xs - 1; i < xs + xm + 1; i++) {
708: struct _LimitInfo info;
709: PetscScalar *cjmpL, *cjmpR;
710: if ((i > sm - 2 && i < mf - lmbwidth + 1) || (i > fm + rmbwidth - 2 && i < ms + 1)) { /* slow components and the first and last fast components */
711: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
712: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
713: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
714: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
715: cjmpL = &ctx->cjmpLR[0];
716: cjmpR = &ctx->cjmpLR[dof];
717: for (j = 0; j < dof; j++) {
718: PetscScalar jmpL, jmpR;
719: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
720: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
721: for (k = 0; k < dof; k++) {
722: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
723: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
724: }
725: }
726: /* Apply limiter to the left and right characteristic jumps */
727: info.m = dof;
728: info.hxs = hxs;
729: info.hxm = hxm;
730: info.hxf = hxf;
731: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
732: for (j = 0; j < dof; j++) {
733: PetscScalar tmp = 0;
734: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
735: slope[i * dof + j] = tmp;
736: }
737: }
738: }
740: for (i = xs; i < xs + xm + 1; i++) {
741: PetscReal maxspeed;
742: PetscScalar *uL, *uR;
743: uL = &ctx->uLR[0];
744: uR = &ctx->uLR[dof];
745: if (i == sm) { /* interface between slow and medium regions */
746: for (j = 0; j < dof; j++) {
747: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
748: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
749: }
750: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
751: if (i < xs + xm) {
752: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
753: imedium++;
754: }
755: }
756: if (i > sm && i < mf - lmbwidth) { /* medium region */
757: for (j = 0; j < dof; j++) {
758: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
759: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
760: }
761: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
762: if (i > xs) {
763: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
764: }
765: if (i < xs + xm) {
766: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
767: imedium++;
768: }
769: }
770: if (i == mf - lmbwidth) { /* interface between medium and fast regions */
771: for (j = 0; j < dof; j++) {
772: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
773: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
774: }
775: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
776: if (i > xs) {
777: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
778: }
779: }
780: if (i == fm + rmbwidth) { /* interface between fast and medium regions */
781: for (j = 0; j < dof; j++) {
782: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
783: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
784: }
785: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
786: if (i < xs + xm) {
787: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
788: imedium++;
789: }
790: }
791: if (i > fm + rmbwidth && i < ms) { /* medium region */
792: for (j = 0; j < dof; j++) {
793: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
794: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
795: }
796: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
797: if (i > xs) {
798: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
799: }
800: if (i < xs + xm) {
801: for (j = 0; j < dof; j++) f[imedium * dof + j] += ctx->flux[j] / hxm;
802: imedium++;
803: }
804: }
805: if (i == ms) { /* interface between medium and slow regions */
806: for (j = 0; j < dof; j++) {
807: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
808: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
809: }
810: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
811: if (i > xs) {
812: for (j = 0; j < dof; j++) f[(imedium - 1) * dof + j] -= ctx->flux[j] / hxm;
813: }
814: }
815: }
816: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
817: PetscCall(VecRestoreArray(F, &f));
818: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
819: PetscCall(DMRestoreLocalVector(da, &Xloc));
820: PetscFunctionReturn(PETSC_SUCCESS);
821: }
823: PetscErrorCode FVRHSFunctionmediumbuffer_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
824: {
825: FVCtx *ctx = (FVCtx *)vctx;
826: PetscInt i, j, k, Mx, dof, xs, xm, imediumbuffer = 0, mf = ctx->mf, fm = ctx->fm, lmbwidth = ctx->lmbwidth, rmbwidth = ctx->rmbwidth;
827: PetscReal hxs, hxm, hxf;
828: PetscScalar *x, *f, *slope;
829: Vec Xloc;
830: DM da;
832: PetscFunctionBeginUser;
833: PetscCall(TSGetDM(ts, &da));
834: PetscCall(DMGetLocalVector(da, &Xloc));
835: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
836: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
837: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
838: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
839: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
840: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
841: PetscCall(VecZeroEntries(F));
842: PetscCall(DMDAVecGetArray(da, Xloc, &x));
843: PetscCall(VecGetArray(F, &f));
844: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
845: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
847: if (ctx->bctype == FVBC_OUTFLOW) {
848: for (i = xs - 2; i < 0; i++) {
849: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
850: }
851: for (i = Mx; i < xs + xm + 2; i++) {
852: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
853: }
854: }
855: for (i = xs - 1; i < xs + xm + 1; i++) {
856: struct _LimitInfo info;
857: PetscScalar *cjmpL, *cjmpR;
858: if ((i > mf - lmbwidth - 2 && i < mf + 1) || (i > fm - 2 && i < fm + rmbwidth + 1)) {
859: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
860: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
861: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
862: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
863: cjmpL = &ctx->cjmpLR[0];
864: cjmpR = &ctx->cjmpLR[dof];
865: for (j = 0; j < dof; j++) {
866: PetscScalar jmpL, jmpR;
867: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
868: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
869: for (k = 0; k < dof; k++) {
870: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
871: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
872: }
873: }
874: /* Apply limiter to the left and right characteristic jumps */
875: info.m = dof;
876: info.hxs = hxs;
877: info.hxm = hxm;
878: info.hxf = hxf;
879: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
880: for (j = 0; j < dof; j++) {
881: PetscScalar tmp = 0;
882: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
883: slope[i * dof + j] = tmp;
884: }
885: }
886: }
888: for (i = xs; i < xs + xm + 1; i++) {
889: PetscReal maxspeed;
890: PetscScalar *uL, *uR;
891: uL = &ctx->uLR[0];
892: uR = &ctx->uLR[dof];
893: if (i == mf - lmbwidth) {
894: for (j = 0; j < dof; j++) {
895: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
896: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
897: }
898: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
899: if (i < xs + xm) {
900: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
901: imediumbuffer++;
902: }
903: }
904: if (i > mf - lmbwidth && i < mf) {
905: for (j = 0; j < dof; j++) {
906: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
907: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
908: }
909: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
910: if (i > xs) {
911: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
912: }
913: if (i < xs + xm) {
914: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
915: imediumbuffer++;
916: }
917: }
918: if (i == mf) { /* interface between the medium region and the fast region */
919: for (j = 0; j < dof; j++) {
920: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
921: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
922: }
923: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
924: if (i > xs) {
925: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
926: }
927: }
928: if (i == fm) { /* interface between the fast region and the medium region */
929: for (j = 0; j < dof; j++) {
930: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
931: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
932: }
933: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
934: if (i < xs + xm) {
935: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
936: imediumbuffer++;
937: }
938: }
939: if (i > fm && i < fm + rmbwidth) {
940: for (j = 0; j < dof; j++) {
941: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
942: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
943: }
944: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
945: if (i > xs) {
946: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
947: }
948: if (i < xs + xm) {
949: for (j = 0; j < dof; j++) f[imediumbuffer * dof + j] += ctx->flux[j] / hxm;
950: imediumbuffer++;
951: }
952: }
953: if (i == fm + rmbwidth) {
954: for (j = 0; j < dof; j++) {
955: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
956: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
957: }
958: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
959: if (i > xs) {
960: for (j = 0; j < dof; j++) f[(imediumbuffer - 1) * dof + j] -= ctx->flux[j] / hxm;
961: }
962: }
963: }
964: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
965: PetscCall(VecRestoreArray(F, &f));
966: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
967: PetscCall(DMRestoreLocalVector(da, &Xloc));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
972: PetscErrorCode FVRHSFunctionfast_3WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
973: {
974: FVCtx *ctx = (FVCtx *)vctx;
975: PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, mf = ctx->mf, fm = ctx->fm;
976: PetscReal hxs, hxm, hxf;
977: PetscScalar *x, *f, *slope;
978: Vec Xloc;
979: DM da;
981: PetscFunctionBeginUser;
982: PetscCall(TSGetDM(ts, &da));
983: PetscCall(DMGetLocalVector(da, &Xloc));
984: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
985: hxs = (ctx->xmax - ctx->xmin) / 8.0 / ctx->sm;
986: hxm = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->mf - ctx->sm);
987: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fm - ctx->mf);
988: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
989: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
990: PetscCall(VecZeroEntries(F));
991: PetscCall(DMDAVecGetArray(da, Xloc, &x));
992: PetscCall(VecGetArray(F, &f));
993: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
994: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
996: if (ctx->bctype == FVBC_OUTFLOW) {
997: for (i = xs - 2; i < 0; i++) {
998: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
999: }
1000: for (i = Mx; i < xs + xm + 2; i++) {
1001: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
1002: }
1003: }
1004: 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 */
1005: struct _LimitInfo info;
1006: PetscScalar *cjmpL, *cjmpR;
1007: if (i > mf - 2 && i < fm + 1) {
1008: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
1009: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
1010: cjmpL = &ctx->cjmpLR[0];
1011: cjmpR = &ctx->cjmpLR[dof];
1012: for (j = 0; j < dof; j++) {
1013: PetscScalar jmpL, jmpR;
1014: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
1015: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
1016: for (k = 0; k < dof; k++) {
1017: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
1018: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
1019: }
1020: }
1021: /* Apply limiter to the left and right characteristic jumps */
1022: info.m = dof;
1023: info.hxs = hxs;
1024: info.hxm = hxm;
1025: info.hxf = hxf;
1026: (*ctx->limit3)(&info, cjmpL, cjmpR, ctx->sm, ctx->mf, ctx->fm, ctx->ms, i, ctx->cslope);
1027: for (j = 0; j < dof; j++) {
1028: PetscScalar tmp = 0;
1029: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
1030: slope[i * dof + j] = tmp;
1031: }
1032: }
1033: }
1035: for (i = xs; i < xs + xm + 1; i++) {
1036: PetscReal maxspeed;
1037: PetscScalar *uL, *uR;
1038: uL = &ctx->uLR[0];
1039: uR = &ctx->uLR[dof];
1040: if (i == mf) { /* interface between medium and fast regions */
1041: for (j = 0; j < dof; j++) {
1042: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxm / 2;
1043: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1044: }
1045: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1046: if (i < xs + xm) {
1047: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1048: ifast++;
1049: }
1050: }
1051: if (i > mf && i < fm) { /* fast region */
1052: for (j = 0; j < dof; j++) {
1053: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1054: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
1055: }
1056: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1057: if (i > xs) {
1058: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1059: }
1060: if (i < xs + xm) {
1061: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
1062: ifast++;
1063: }
1064: }
1065: if (i == fm) { /* interface between fast and medium regions */
1066: for (j = 0; j < dof; j++) {
1067: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
1068: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxm / 2;
1069: }
1070: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
1071: if (i > xs) {
1072: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
1073: }
1074: }
1075: }
1076: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
1077: PetscCall(VecRestoreArray(F, &f));
1078: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
1079: PetscCall(DMRestoreLocalVector(da, &Xloc));
1080: PetscFunctionReturn(PETSC_SUCCESS);
1081: }
1083: int main(int argc, char *argv[])
1084: {
1085: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
1086: PetscFunctionList limiters = 0, physics = 0;
1087: MPI_Comm comm;
1088: TS ts;
1089: DM da;
1090: Vec X, X0, R;
1091: FVCtx ctx;
1092: 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;
1093: PetscBool view_final = PETSC_FALSE;
1094: PetscReal ptime;
1096: PetscFunctionBeginUser;
1097: PetscCall(PetscInitialize(&argc, &argv, 0, help));
1098: comm = PETSC_COMM_WORLD;
1099: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
1101: /* Register limiters to be available on the command line */
1102: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit3_Upwind));
1103: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit3_LaxWendroff));
1104: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit3_BeamWarming));
1105: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit3_Fromm));
1106: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit3_Minmod));
1107: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit3_Superbee));
1108: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit3_MC));
1109: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit3_Koren3));
1111: /* Register physical models to be available on the command line */
1112: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
1114: ctx.comm = comm;
1115: ctx.cfl = 0.9;
1116: ctx.bctype = FVBC_PERIODIC;
1117: ctx.xmin = -1.0;
1118: ctx.xmax = 1.0;
1119: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
1120: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
1121: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
1122: PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
1123: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
1124: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
1125: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
1126: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
1127: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
1128: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
1129: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
1130: PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
1131: PetscOptionsEnd();
1133: /* Choose the limiter from the list of registered limiters */
1134: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit3));
1135: PetscCheck(ctx.limit3, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
1137: /* Choose the physics from the list of registered models */
1138: {
1139: PetscErrorCode (*r)(FVCtx *);
1140: PetscCall(PetscFunctionListFind(physics, physname, &r));
1141: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
1142: /* Create the physics, will set the number of fields and their names */
1143: PetscCall((*r)(&ctx));
1144: }
1146: /* Create a DMDA to manage the parallel grid */
1147: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da));
1148: PetscCall(DMSetFromOptions(da));
1149: PetscCall(DMSetUp(da));
1150: /* Inform the DMDA of the field names provided by the physics. */
1151: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
1152: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i]));
1153: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
1154: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
1156: /* Set coordinates of cell centers */
1157: 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));
1159: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
1160: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
1161: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
1163: /* Create a vector to store the solution and to save the initial state */
1164: PetscCall(DMCreateGlobalVector(da, &X));
1165: PetscCall(VecDuplicate(X, &X0));
1166: PetscCall(VecDuplicate(X, &R));
1168: /* create index for slow parts and fast parts,
1169: count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
1170: count_slow = Mx / (1 + ctx.hratio) / (1 + ctx.hratio);
1171: count_medium = 2 * ctx.hratio * count_slow;
1172: 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");
1173: count_fast = ctx.hratio * ctx.hratio * count_slow;
1174: ctx.sm = count_slow / 2;
1175: ctx.mf = ctx.sm + count_medium / 2;
1176: ctx.fm = ctx.mf + count_fast;
1177: ctx.ms = ctx.fm + count_medium / 2;
1178: PetscCall(PetscMalloc1(xm * dof, &index_slow));
1179: PetscCall(PetscMalloc1(xm * dof, &index_medium));
1180: PetscCall(PetscMalloc1(xm * dof, &index_fast));
1181: PetscCall(PetscMalloc1(6 * dof, &index_slowbuffer));
1182: PetscCall(PetscMalloc1(6 * dof, &index_mediumbuffer));
1183: if (((AdvectCtx *)ctx.physics2.user)->a > 0) {
1184: ctx.lsbwidth = 2;
1185: ctx.rsbwidth = 4;
1186: ctx.lmbwidth = 2;
1187: ctx.rmbwidth = 4;
1188: } else {
1189: ctx.lsbwidth = 4;
1190: ctx.rsbwidth = 2;
1191: ctx.lmbwidth = 4;
1192: ctx.rmbwidth = 2;
1193: }
1195: for (i = xs; i < xs + xm; i++) {
1196: if (i < ctx.sm - ctx.lsbwidth || i > ctx.ms + ctx.rsbwidth - 1)
1197: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
1198: else if ((i >= ctx.sm - ctx.lsbwidth && i < ctx.sm) || (i > ctx.ms - 1 && i <= ctx.ms + ctx.rsbwidth - 1))
1199: for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
1200: else if (i < ctx.mf - ctx.lmbwidth || i > ctx.fm + ctx.rmbwidth - 1)
1201: for (k = 0; k < dof; k++) index_medium[imedium++] = i * dof + k;
1202: else if ((i >= ctx.mf - ctx.lmbwidth && i < ctx.mf) || (i > ctx.fm - 1 && i <= ctx.fm + ctx.rmbwidth - 1))
1203: for (k = 0; k < dof; k++) index_mediumbuffer[imediumbuffer++] = i * dof + k;
1204: else
1205: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
1206: }
1207: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
1208: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, imedium, index_medium, PETSC_COPY_VALUES, &ctx.ism));
1209: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
1210: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));
1211: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, imediumbuffer, index_mediumbuffer, PETSC_COPY_VALUES, &ctx.ismb));
1213: /* Create a time-stepping object */
1214: PetscCall(TSCreate(comm, &ts));
1215: PetscCall(TSSetDM(ts, da));
1216: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_3WaySplit, &ctx));
1217: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
1218: PetscCall(TSRHSSplitSetIS(ts, "medium", ctx.ism));
1219: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
1220: PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
1221: PetscCall(TSRHSSplitSetIS(ts, "mediumbuffer", ctx.ismb));
1222: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_3WaySplit, &ctx));
1223: PetscCall(TSRHSSplitSetRHSFunction(ts, "medium", NULL, FVRHSFunctionmedium_3WaySplit, &ctx));
1224: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_3WaySplit, &ctx));
1225: PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_3WaySplit, &ctx));
1226: PetscCall(TSRHSSplitSetRHSFunction(ts, "mediumbuffer", NULL, FVRHSFunctionmediumbuffer_3WaySplit, &ctx));
1228: PetscCall(TSSetType(ts, TSSSP));
1229: /*PetscCall(TSSetType(ts,TSMPRK));*/
1230: PetscCall(TSSetMaxTime(ts, 10));
1231: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1233: /* Compute initial conditions and starting time step */
1234: PetscCall(FVSample_3WaySplit(&ctx, da, 0, X0));
1235: PetscCall(FVRHSFunction_3WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
1236: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
1237: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
1238: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
1239: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1240: {
1241: PetscInt steps;
1242: PetscScalar mass_initial, mass_final, mass_difference;
1243: const PetscScalar *ptr_X, *ptr_X0;
1244: const PetscReal hs = (ctx.xmax - ctx.xmin) / 4.0 / count_slow;
1245: const PetscReal hm = (ctx.xmax - ctx.xmin) / 2.0 / count_medium;
1246: const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;
1248: PetscCall(TSSolve(ts, X));
1249: PetscCall(TSGetSolveTime(ts, &ptime));
1250: PetscCall(TSGetStepNumber(ts, &steps));
1251: /* calculate the total mass at initial time and final time */
1252: mass_initial = 0.0;
1253: mass_final = 0.0;
1254: PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
1255: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
1256: for (i = xs; i < xs + xm; i++) {
1257: if (i < ctx.sm || i > ctx.ms - 1)
1258: for (k = 0; k < dof; k++) {
1259: mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
1260: mass_final = mass_final + hs * ptr_X[i * dof + k];
1261: }
1262: else if (i < ctx.mf || i > ctx.fm - 1)
1263: for (k = 0; k < dof; k++) {
1264: mass_initial = mass_initial + hm * ptr_X0[i * dof + k];
1265: mass_final = mass_final + hm * ptr_X[i * dof + k];
1266: }
1267: else {
1268: for (k = 0; k < dof; k++) {
1269: mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
1270: mass_final = mass_final + hf * ptr_X[i * dof + k];
1271: }
1272: }
1273: }
1274: PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
1275: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
1276: mass_difference = mass_final - mass_initial;
1277: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mass_difference, 1, MPIU_SCALAR, MPIU_SUM, comm));
1278: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
1279: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
1280: PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1 / ctx.cfl_idt)));
1281: if (ctx.exact) {
1282: PetscReal nrm1 = 0;
1283: PetscCall(SolutionErrorNorms_3WaySplit(&ctx, da, ptime, X, &nrm1));
1284: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
1285: }
1286: if (ctx.simulation) {
1287: PetscReal nrm1 = 0;
1288: PetscViewer fd;
1289: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
1290: Vec XR;
1291: PetscBool flg;
1292: const PetscScalar *ptr_XR;
1293: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
1294: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
1295: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
1296: PetscCall(VecDuplicate(X0, &XR));
1297: PetscCall(VecLoad(XR, fd));
1298: PetscCall(PetscViewerDestroy(&fd));
1299: PetscCall(VecGetArrayRead(X, &ptr_X));
1300: PetscCall(VecGetArrayRead(XR, &ptr_XR));
1301: for (i = xs; i < xs + xm; i++) {
1302: if (i < ctx.sm || i > ctx.ms - 1)
1303: for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1304: else if (i < ctx.mf || i < ctx.fm - 1)
1305: for (k = 0; k < dof; k++) nrm1 = nrm1 + hm * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1306: else
1307: for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
1308: }
1309: PetscCall(VecRestoreArrayRead(X, &ptr_X));
1310: PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
1311: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
1312: PetscCall(VecDestroy(&XR));
1313: }
1314: }
1316: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
1317: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
1318: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
1319: if (draw & 0x4) {
1320: Vec Y;
1321: PetscCall(VecDuplicate(X, &Y));
1322: PetscCall(FVSample_3WaySplit(&ctx, da, ptime, Y));
1323: PetscCall(VecAYPX(Y, -1, X));
1324: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
1325: PetscCall(VecDestroy(&Y));
1326: }
1328: if (view_final) {
1329: PetscViewer viewer;
1330: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
1331: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
1332: PetscCall(VecView(X, viewer));
1333: PetscCall(PetscViewerPopFormat(viewer));
1334: PetscCall(PetscViewerDestroy(&viewer));
1335: }
1337: /* Clean up */
1338: PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
1339: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
1340: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
1341: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
1342: PetscCall(VecDestroy(&X));
1343: PetscCall(VecDestroy(&X0));
1344: PetscCall(VecDestroy(&R));
1345: PetscCall(DMDestroy(&da));
1346: PetscCall(TSDestroy(&ts));
1347: PetscCall(ISDestroy(&ctx.iss));
1348: PetscCall(ISDestroy(&ctx.ism));
1349: PetscCall(ISDestroy(&ctx.isf));
1350: PetscCall(ISDestroy(&ctx.issb));
1351: PetscCall(ISDestroy(&ctx.ismb));
1352: PetscCall(PetscFree(index_slow));
1353: PetscCall(PetscFree(index_medium));
1354: PetscCall(PetscFree(index_fast));
1355: PetscCall(PetscFree(index_slowbuffer));
1356: PetscCall(PetscFree(index_mediumbuffer));
1357: PetscCall(PetscFunctionListDestroy(&limiters));
1358: PetscCall(PetscFunctionListDestroy(&physics));
1359: PetscCall(PetscFinalize());
1360: return 0;
1361: }
1363: /*TEST
1365: build:
1366: requires: !complex
1367: depends: finitevolume1d.c
1369: test:
1370: suffix: 1
1371: 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
1373: test:
1374: suffix: 2
1375: 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
1377: TEST*/