Actual source code: ex6.c
1: /*
2: Note:
3: -hratio is the ratio between mesh size of coarse grids and fine grids
4: -ts_rk_dtratio is the ratio between the large stepsize and the small stepsize
5: */
7: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
8: " advection - Constant coefficient scalar advection\n"
9: " u_t + (a*u)_x = 0\n"
10: " for this toy problem, we choose different meshsizes for different sub-domains (slow-fast-slow), say\n"
11: " hxs = hratio*hxf \n"
12: " where hxs and hxf are the grid spacings for coarse and fine grids respectively.\n"
13: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
14: " the states across shocks and rarefactions\n"
15: " simulation - use reference solution which is generated by smaller time step size to be true solution,\n"
16: " also the reference solution should be generated by user and stored in a binary file.\n"
17: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
18: "Several initial conditions can be chosen with -initial N\n\n"
19: "The problem size should be set with -da_grid_x M\n\n";
21: #include <petscts.h>
22: #include <petscdm.h>
23: #include <petscdmda.h>
24: #include <petscdraw.h>
25: #include "finitevolume1d.h"
27: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
28: {
29: PetscReal range = xmax - xmin;
30: return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
31: }
33: /* --------------------------------- Advection ----------------------------------- */
34: typedef struct {
35: PetscReal a; /* advective velocity */
36: } AdvectCtx;
38: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed)
39: {
40: AdvectCtx *ctx = (AdvectCtx *)vctx;
41: PetscReal speed;
43: PetscFunctionBeginUser;
44: speed = ctx->a;
45: flux[0] = PetscMax(0, speed) * uL[0] + PetscMin(0, speed) * uR[0];
46: *maxspeed = speed;
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds)
51: {
52: AdvectCtx *ctx = (AdvectCtx *)vctx;
54: PetscFunctionBeginUser;
55: X[0] = 1.;
56: Xi[0] = 1.;
57: speeds[0] = ctx->a;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
62: {
63: AdvectCtx *ctx = (AdvectCtx *)vctx;
64: PetscReal a = ctx->a, x0;
66: PetscFunctionBeginUser;
67: switch (bctype) {
68: case FVBC_OUTFLOW:
69: x0 = x - a * t;
70: break;
71: case FVBC_PERIODIC:
72: x0 = RangeMod(x - a * t, xmin, xmax);
73: break;
74: default:
75: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
76: }
77: switch (initial) {
78: case 0:
79: u[0] = (x0 < 0) ? 1 : -1;
80: break;
81: case 1:
82: u[0] = (x0 < 0) ? -1 : 1;
83: break;
84: case 2:
85: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
86: break;
87: case 3:
88: u[0] = PetscSinReal(2 * PETSC_PI * x0);
89: break;
90: case 4:
91: u[0] = PetscAbs(x0);
92: break;
93: case 5:
94: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
95: break;
96: case 6:
97: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
98: break;
99: case 7:
100: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
101: break;
102: default:
103: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
104: }
105: PetscFunctionReturn(PETSC_SUCCESS);
106: }
108: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
109: {
110: AdvectCtx *user;
112: PetscFunctionBeginUser;
113: PetscCall(PetscNew(&user));
114: ctx->physics2.sample2 = PhysicsSample_Advect;
115: ctx->physics2.riemann2 = PhysicsRiemann_Advect;
116: ctx->physics2.characteristic2 = PhysicsCharacteristic_Advect;
117: ctx->physics2.destroy = PhysicsDestroy_SimpleFree;
118: ctx->physics2.user = user;
119: ctx->physics2.dof = 1;
120: PetscCall(PetscStrallocpy("u", &ctx->physics2.fieldname[0]));
121: user->a = 1;
122: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
123: {
124: PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
125: }
126: PetscOptionsEnd();
127: PetscFunctionReturn(PETSC_SUCCESS);
128: }
130: PetscErrorCode FVSample_2WaySplit(FVCtx *ctx, DM da, PetscReal time, Vec U)
131: {
132: PetscScalar *u, *uj, xj, xi;
133: PetscInt i, j, k, dof, xs, xm, Mx;
134: const PetscInt N = 200;
135: PetscReal hs, hf;
137: PetscFunctionBeginUser;
138: PetscCheck(ctx->physics2.sample2, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
139: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
140: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
141: PetscCall(DMDAVecGetArray(da, U, &u));
142: PetscCall(PetscMalloc1(dof, &uj));
143: hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
144: hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
145: for (i = xs; i < xs + xm; i++) {
146: if (i < ctx->sf) {
147: xi = ctx->xmin + 0.5 * hs + i * hs;
148: /* Integrate over cell i using trapezoid rule with N points. */
149: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
150: for (j = 0; j < N + 1; j++) {
151: xj = xi + hs * (j - N / 2) / (PetscReal)N;
152: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
153: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
154: }
155: } else if (i < ctx->fs) {
156: xi = ctx->xmin + ctx->sf * hs + 0.5 * hf + (i - ctx->sf) * hf;
157: /* Integrate over cell i using trapezoid rule with N points. */
158: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
159: for (j = 0; j < N + 1; j++) {
160: xj = xi + hf * (j - N / 2) / (PetscReal)N;
161: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
162: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
163: }
164: } else {
165: xi = ctx->xmin + ctx->sf * hs + (ctx->fs - ctx->sf) * hf + 0.5 * hs + (i - ctx->fs) * hs;
166: /* Integrate over cell i using trapezoid rule with N points. */
167: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
168: for (j = 0; j < N + 1; j++) {
169: xj = xi + hs * (j - N / 2) / (PetscReal)N;
170: PetscCall((*ctx->physics2.sample2)(ctx->physics2.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
171: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
172: }
173: }
174: }
175: PetscCall(DMDAVecRestoreArray(da, U, &u));
176: PetscCall(PetscFree(uj));
177: PetscFunctionReturn(PETSC_SUCCESS);
178: }
180: static PetscErrorCode SolutionErrorNorms_2WaySplit(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
181: {
182: Vec Y;
183: PetscInt i, Mx;
184: const PetscScalar *ptr_X, *ptr_Y;
185: PetscReal hs, hf;
187: PetscFunctionBeginUser;
188: PetscCall(VecGetSize(X, &Mx));
189: PetscCall(VecDuplicate(X, &Y));
190: PetscCall(FVSample_2WaySplit(ctx, da, t, Y));
191: hs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
192: hf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
193: PetscCall(VecGetArrayRead(X, &ptr_X));
194: PetscCall(VecGetArrayRead(Y, &ptr_Y));
195: for (i = 0; i < Mx; i++) {
196: if (i < ctx->sf || i > ctx->fs - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
197: else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
198: }
199: PetscCall(VecRestoreArrayRead(X, &ptr_X));
200: PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
201: PetscCall(VecDestroy(&Y));
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: PetscErrorCode FVRHSFunction_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
206: {
207: FVCtx *ctx = (FVCtx *)vctx;
208: PetscInt i, j, k, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
209: PetscReal hxf, hxs, cfl_idt = 0;
210: PetscScalar *x, *f, *slope;
211: Vec Xloc;
212: DM da;
214: PetscFunctionBeginUser;
215: PetscCall(TSGetDM(ts, &da));
216: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
217: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */
218: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
219: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
220: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
221: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
223: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
225: PetscCall(DMDAVecGetArray(da, Xloc, &x));
226: PetscCall(DMDAVecGetArray(da, F, &f));
227: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
229: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
231: if (ctx->bctype == FVBC_OUTFLOW) {
232: for (i = xs - 2; i < 0; i++) {
233: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
234: }
235: for (i = Mx; i < xs + xm + 2; i++) {
236: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
237: }
238: }
239: for (i = xs - 1; i < xs + xm + 1; i++) {
240: struct _LimitInfo info;
241: PetscScalar *cjmpL, *cjmpR;
242: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
243: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
244: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
245: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
246: cjmpL = &ctx->cjmpLR[0];
247: cjmpR = &ctx->cjmpLR[dof];
248: for (j = 0; j < dof; j++) {
249: PetscScalar jmpL, jmpR;
250: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
251: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
252: for (k = 0; k < dof; k++) {
253: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
254: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
255: }
256: }
257: /* Apply limiter to the left and right characteristic jumps */
258: info.m = dof;
259: info.hxs = hxs;
260: info.hxf = hxf;
261: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
262: for (j = 0; j < dof; j++) {
263: PetscScalar tmp = 0;
264: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
265: slope[i * dof + j] = tmp;
266: }
267: }
269: for (i = xs; i < xs + xm + 1; i++) {
270: PetscReal maxspeed;
271: PetscScalar *uL, *uR;
272: uL = &ctx->uLR[0];
273: uR = &ctx->uLR[dof];
274: if (i < sf) { /* slow region */
275: for (j = 0; j < dof; j++) {
276: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
277: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
278: }
279: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
280: if (i > xs) {
281: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
282: }
283: if (i < xs + xm) {
284: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
285: }
286: } else if (i == sf) { /* interface between the slow region and the fast region */
287: for (j = 0; j < dof; j++) {
288: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
289: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
290: }
291: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
292: if (i > xs) {
293: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
294: }
295: if (i < xs + xm) {
296: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
297: }
298: } else if (i > sf && i < fs) { /* fast region */
299: for (j = 0; j < dof; j++) {
300: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
301: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
302: }
303: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
304: if (i > xs) {
305: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
306: }
307: if (i < xs + xm) {
308: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
309: }
310: } else if (i == fs) { /* interface between the fast region and the slow region */
311: for (j = 0; j < dof; j++) {
312: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
313: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
314: }
315: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
316: if (i > xs) {
317: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
318: }
319: if (i < xs + xm) {
320: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
321: }
322: } else { /* slow region */
323: for (j = 0; j < dof; j++) {
324: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
325: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
326: }
327: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
328: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
329: if (i > xs) {
330: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
331: }
332: if (i < xs + xm) {
333: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
334: }
335: }
336: }
337: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
338: PetscCall(DMDAVecRestoreArray(da, F, &f));
339: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
340: PetscCall(DMRestoreLocalVector(da, &Xloc));
341: PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
342: if (0) {
343: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
344: PetscReal dt, tnow;
345: PetscCall(TSGetTimeStep(ts, &dt));
346: PetscCall(TSGetTime(ts, &tnow));
347: if (dt > 0.5 / ctx->cfl_idt) {
348: if (1) PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(1 / (2 * ctx->cfl_idt))));
349: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
350: }
351: }
352: PetscFunctionReturn(PETSC_SUCCESS);
353: }
355: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
356: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
357: {
358: FVCtx *ctx = (FVCtx *)vctx;
359: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
360: PetscReal hxs, hxf, cfl_idt = 0;
361: PetscScalar *x, *f, *slope;
362: Vec Xloc;
363: DM da;
365: PetscFunctionBeginUser;
366: PetscCall(TSGetDM(ts, &da));
367: PetscCall(DMGetLocalVector(da, &Xloc));
368: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
369: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
370: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
371: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
372: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
373: PetscCall(VecZeroEntries(F));
374: PetscCall(DMDAVecGetArray(da, Xloc, &x));
375: PetscCall(VecGetArray(F, &f));
376: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
377: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
379: if (ctx->bctype == FVBC_OUTFLOW) {
380: for (i = xs - 2; i < 0; i++) {
381: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
382: }
383: for (i = Mx; i < xs + xm + 2; i++) {
384: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
385: }
386: }
387: for (i = xs - 1; i < xs + xm + 1; i++) {
388: struct _LimitInfo info;
389: PetscScalar *cjmpL, *cjmpR;
390: if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */
391: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
392: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
393: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
394: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
395: cjmpL = &ctx->cjmpLR[0];
396: cjmpR = &ctx->cjmpLR[dof];
397: for (j = 0; j < dof; j++) {
398: PetscScalar jmpL, jmpR;
399: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
400: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
401: for (k = 0; k < dof; k++) {
402: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
403: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
404: }
405: }
406: /* Apply limiter to the left and right characteristic jumps */
407: info.m = dof;
408: info.hxs = hxs;
409: info.hxf = hxf;
410: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
411: for (j = 0; j < dof; j++) {
412: PetscScalar tmp = 0;
413: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
414: slope[i * dof + j] = tmp;
415: }
416: }
417: }
419: for (i = xs; i < xs + xm + 1; i++) {
420: PetscReal maxspeed;
421: PetscScalar *uL, *uR;
422: uL = &ctx->uLR[0];
423: uR = &ctx->uLR[dof];
424: if (i < sf - lsbwidth) { /* slow region */
425: for (j = 0; j < dof; j++) {
426: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
427: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
428: }
429: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
430: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
431: if (i > xs) {
432: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
433: }
434: if (i < xs + xm) {
435: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
436: islow++;
437: }
438: }
439: if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */
440: for (j = 0; j < dof; j++) {
441: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
442: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
443: }
444: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
445: if (i > xs) {
446: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
447: }
448: }
449: if (i == fs + rsbwidth) { /* slow region */
450: for (j = 0; j < dof; j++) {
451: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
452: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
453: }
454: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
455: if (i < xs + xm) {
456: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
457: islow++;
458: }
459: }
460: if (i > fs + rsbwidth) { /* slow region */
461: for (j = 0; j < dof; j++) {
462: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
463: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
464: }
465: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
466: if (i > xs) {
467: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
468: }
469: if (i < xs + xm) {
470: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
471: islow++;
472: }
473: }
474: }
475: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
476: PetscCall(VecRestoreArray(F, &f));
477: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
478: PetscCall(DMRestoreLocalVector(da, &Xloc));
479: PetscCallMPI(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
480: PetscFunctionReturn(PETSC_SUCCESS);
481: }
483: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
484: {
485: FVCtx *ctx = (FVCtx *)vctx;
486: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
487: PetscReal hxs, hxf;
488: PetscScalar *x, *f, *slope;
489: Vec Xloc;
490: DM da;
492: PetscFunctionBeginUser;
493: PetscCall(TSGetDM(ts, &da));
494: PetscCall(DMGetLocalVector(da, &Xloc));
495: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
496: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
497: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
498: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
499: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
500: PetscCall(VecZeroEntries(F));
501: PetscCall(DMDAVecGetArray(da, Xloc, &x));
502: PetscCall(VecGetArray(F, &f));
503: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
504: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
506: if (ctx->bctype == FVBC_OUTFLOW) {
507: for (i = xs - 2; i < 0; i++) {
508: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
509: }
510: for (i = Mx; i < xs + xm + 2; i++) {
511: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
512: }
513: }
514: for (i = xs - 1; i < xs + xm + 1; i++) {
515: struct _LimitInfo info;
516: PetscScalar *cjmpL, *cjmpR;
517: if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) {
518: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
519: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
520: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
521: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
522: cjmpL = &ctx->cjmpLR[0];
523: cjmpR = &ctx->cjmpLR[dof];
524: for (j = 0; j < dof; j++) {
525: PetscScalar jmpL, jmpR;
526: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
527: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
528: for (k = 0; k < dof; k++) {
529: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
530: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
531: }
532: }
533: /* Apply limiter to the left and right characteristic jumps */
534: info.m = dof;
535: info.hxs = hxs;
536: info.hxf = hxf;
537: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
538: for (j = 0; j < dof; j++) {
539: PetscScalar tmp = 0;
540: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
541: slope[i * dof + j] = tmp;
542: }
543: }
544: }
546: for (i = xs; i < xs + xm + 1; i++) {
547: PetscReal maxspeed;
548: PetscScalar *uL, *uR;
549: uL = &ctx->uLR[0];
550: uR = &ctx->uLR[dof];
551: if (i == sf - lsbwidth) {
552: for (j = 0; j < dof; j++) {
553: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
554: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
555: }
556: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
557: if (i < xs + xm) {
558: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
559: islow++;
560: }
561: }
562: if (i > sf - lsbwidth && i < sf) {
563: for (j = 0; j < dof; j++) {
564: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
565: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
566: }
567: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
568: if (i > xs) {
569: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
570: }
571: if (i < xs + xm) {
572: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
573: islow++;
574: }
575: }
576: if (i == sf) { /* interface between the slow region and the fast region */
577: for (j = 0; j < dof; j++) {
578: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
579: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
580: }
581: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
582: if (i > xs) {
583: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
584: }
585: }
586: if (i == fs) { /* interface between the fast region and the slow region */
587: for (j = 0; j < dof; j++) {
588: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
589: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
590: }
591: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
592: if (i < xs + xm) {
593: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
594: islow++;
595: }
596: }
597: if (i > fs && i < fs + rsbwidth) {
598: for (j = 0; j < dof; j++) {
599: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
600: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
601: }
602: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
603: if (i > xs) {
604: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
605: }
606: if (i < xs + xm) {
607: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
608: islow++;
609: }
610: }
611: if (i == fs + rsbwidth) {
612: for (j = 0; j < dof; j++) {
613: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
614: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
615: }
616: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
617: if (i > xs) {
618: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
619: }
620: }
621: }
622: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
623: PetscCall(VecRestoreArray(F, &f));
624: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
625: PetscCall(DMRestoreLocalVector(da, &Xloc));
626: PetscFunctionReturn(PETSC_SUCCESS);
627: }
629: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
630: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
631: {
632: FVCtx *ctx = (FVCtx *)vctx;
633: PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
634: PetscReal hxs, hxf;
635: PetscScalar *x, *f, *slope;
636: Vec Xloc;
637: DM da;
639: PetscFunctionBeginUser;
640: PetscCall(TSGetDM(ts, &da));
641: PetscCall(DMGetLocalVector(da, &Xloc));
642: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
643: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
644: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
645: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
646: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
647: PetscCall(VecZeroEntries(F));
648: PetscCall(DMDAVecGetArray(da, Xloc, &x));
649: PetscCall(VecGetArray(F, &f));
650: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
651: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
653: if (ctx->bctype == FVBC_OUTFLOW) {
654: for (i = xs - 2; i < 0; i++) {
655: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
656: }
657: for (i = Mx; i < xs + xm + 2; i++) {
658: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
659: }
660: }
661: for (i = xs - 1; i < xs + xm + 1; i++) {
662: struct _LimitInfo info;
663: PetscScalar *cjmpL, *cjmpR;
664: if (i > sf - 2 && i < fs + 1) {
665: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
666: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
667: cjmpL = &ctx->cjmpLR[0];
668: cjmpR = &ctx->cjmpLR[dof];
669: for (j = 0; j < dof; j++) {
670: PetscScalar jmpL, jmpR;
671: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
672: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
673: for (k = 0; k < dof; k++) {
674: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
675: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
676: }
677: }
678: /* Apply limiter to the left and right characteristic jumps */
679: info.m = dof;
680: info.hxs = hxs;
681: info.hxf = hxf;
682: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
683: for (j = 0; j < dof; j++) {
684: PetscScalar tmp = 0;
685: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
686: slope[i * dof + j] = tmp;
687: }
688: }
689: }
691: for (i = xs; i < xs + xm + 1; i++) {
692: PetscReal maxspeed;
693: PetscScalar *uL, *uR;
694: uL = &ctx->uLR[0];
695: uR = &ctx->uLR[dof];
696: if (i == sf) { /* interface between the slow region and the fast region */
697: for (j = 0; j < dof; j++) {
698: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
699: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
700: }
701: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
702: if (i < xs + xm) {
703: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
704: ifast++;
705: }
706: }
707: if (i > sf && i < fs) { /* fast region */
708: for (j = 0; j < dof; j++) {
709: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
710: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
711: }
712: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
713: if (i > xs) {
714: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
715: }
716: if (i < xs + xm) {
717: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
718: ifast++;
719: }
720: }
721: if (i == fs) { /* interface between the fast region and the slow region */
722: for (j = 0; j < dof; j++) {
723: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
724: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
725: }
726: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
727: if (i > xs) {
728: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
729: }
730: }
731: }
732: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
733: PetscCall(VecRestoreArray(F, &f));
734: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
735: PetscCall(DMRestoreLocalVector(da, &Xloc));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: int main(int argc, char *argv[])
740: {
741: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
742: PetscFunctionList limiters = 0, physics = 0;
743: MPI_Comm comm;
744: TS ts;
745: DM da;
746: Vec X, X0, R;
747: FVCtx ctx;
748: PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, islowbuffer = 0, *index_slow, *index_fast, *index_slowbuffer;
749: PetscBool view_final = PETSC_FALSE;
750: PetscReal ptime;
752: PetscFunctionBeginUser;
753: PetscCall(PetscInitialize(&argc, &argv, 0, help));
754: comm = PETSC_COMM_WORLD;
755: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
757: /* Register limiters to be available on the command line */
758: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind));
759: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff));
760: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming));
761: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm));
762: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod));
763: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee));
764: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC));
765: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3));
767: /* Register physical models to be available on the command line */
768: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
770: ctx.comm = comm;
771: ctx.cfl = 0.9;
772: ctx.bctype = FVBC_PERIODIC;
773: ctx.xmin = -1.0;
774: ctx.xmax = 1.0;
775: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
776: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
777: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
778: PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
779: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
780: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
781: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
782: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
783: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
784: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
785: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
786: PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
787: PetscOptionsEnd();
789: /* Choose the limiter from the list of registered limiters */
790: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit2));
791: PetscCheck(ctx.limit2, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
793: /* Choose the physics from the list of registered models */
794: {
795: PetscErrorCode (*r)(FVCtx *);
796: PetscCall(PetscFunctionListFind(physics, physname, &r));
797: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
798: /* Create the physics, will set the number of fields and their names */
799: PetscCall((*r)(&ctx));
800: }
802: /* Create a DMDA to manage the parallel grid */
803: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da));
804: PetscCall(DMSetFromOptions(da));
805: PetscCall(DMSetUp(da));
806: /* Inform the DMDA of the field names provided by the physics. */
807: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
808: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i]));
809: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
810: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
812: /* Set coordinates of cell centers */
813: 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));
815: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
816: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
817: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
819: /* Create a vector to store the solution and to save the initial state */
820: PetscCall(DMCreateGlobalVector(da, &X));
821: PetscCall(VecDuplicate(X, &X0));
822: PetscCall(VecDuplicate(X, &R));
824: /* create index for slow parts and fast parts,
825: count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
826: count_slow = Mx / (1.0 + ctx.hratio / 3.0);
827: PetscCheck(count_slow % 2 == 0, PETSC_COMM_WORLD, PETSC_ERR_USER, "Please adjust grid size Mx (-da_grid_x) and hratio (-hratio) so that Mx/(1+hartio/3) is even");
828: count_fast = Mx - count_slow;
829: ctx.sf = count_slow / 2;
830: ctx.fs = ctx.sf + count_fast;
831: PetscCall(PetscMalloc1(xm * dof, &index_slow));
832: PetscCall(PetscMalloc1(xm * dof, &index_fast));
833: PetscCall(PetscMalloc1(6 * dof, &index_slowbuffer));
834: if (((AdvectCtx *)ctx.physics2.user)->a > 0) {
835: ctx.lsbwidth = 2;
836: ctx.rsbwidth = 4;
837: } else {
838: ctx.lsbwidth = 4;
839: ctx.rsbwidth = 2;
840: }
841: for (i = xs; i < xs + xm; i++) {
842: if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1)
843: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
844: else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1))
845: for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
846: else
847: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
848: }
849: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
850: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
851: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));
853: /* Create a time-stepping object */
854: PetscCall(TSCreate(comm, &ts));
855: PetscCall(TSSetDM(ts, da));
856: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx));
857: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
858: PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
859: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
860: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx));
861: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx));
862: PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx));
864: PetscCall(TSSetType(ts, TSSSP));
865: /*PetscCall(TSSetType(ts,TSMPRK));*/
866: PetscCall(TSSetMaxTime(ts, 10));
867: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
869: /* Compute initial conditions and starting time step */
870: PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0));
871: PetscCall(FVRHSFunction_2WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
872: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
873: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
874: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
875: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
876: {
877: PetscInt steps;
878: PetscScalar mass_initial, mass_final, mass_difference, mass_differenceg;
879: const PetscScalar *ptr_X, *ptr_X0;
880: const PetscReal hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow;
881: const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;
883: PetscCall(TSSolve(ts, X));
884: PetscCall(TSGetSolveTime(ts, &ptime));
885: PetscCall(TSGetStepNumber(ts, &steps));
886: /* calculate the total mass at initial time and final time */
887: mass_initial = 0.0;
888: mass_final = 0.0;
889: PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
890: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
891: for (i = xs; i < xs + xm; i++) {
892: if (i < ctx.sf || i > ctx.fs - 1) {
893: for (k = 0; k < dof; k++) {
894: mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
895: mass_final = mass_final + hs * ptr_X[i * dof + k];
896: }
897: } else {
898: for (k = 0; k < dof; k++) {
899: mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
900: mass_final = mass_final + hf * ptr_X[i * dof + k];
901: }
902: }
903: }
904: PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
905: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
906: mass_difference = mass_final - mass_initial;
907: PetscCallMPI(MPIU_Allreduce(&mass_difference, &mass_differenceg, 1, MPIU_SCALAR, MPIU_SUM, comm));
908: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_differenceg));
909: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
910: PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1 / ctx.cfl_idt)));
911: if (ctx.exact) {
912: PetscReal nrm1 = 0;
913: PetscCall(SolutionErrorNorms_2WaySplit(&ctx, da, ptime, X, &nrm1));
914: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
915: }
916: if (ctx.simulation) {
917: PetscReal nrm1 = 0;
918: PetscViewer fd;
919: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
920: Vec XR;
921: PetscBool flg;
922: const PetscScalar *ptr_XR;
923: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
924: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
925: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
926: PetscCall(VecDuplicate(X0, &XR));
927: PetscCall(VecLoad(XR, fd));
928: PetscCall(PetscViewerDestroy(&fd));
929: PetscCall(VecGetArrayRead(X, &ptr_X));
930: PetscCall(VecGetArrayRead(XR, &ptr_XR));
931: for (i = xs; i < xs + xm; i++) {
932: if (i < ctx.sf || i > ctx.fs - 1)
933: for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
934: else
935: for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
936: }
937: PetscCall(VecRestoreArrayRead(X, &ptr_X));
938: PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
939: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
940: PetscCall(VecDestroy(&XR));
941: }
942: }
944: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
945: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
946: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
947: if (draw & 0x4) {
948: Vec Y;
949: PetscCall(VecDuplicate(X, &Y));
950: PetscCall(FVSample_2WaySplit(&ctx, da, ptime, Y));
951: PetscCall(VecAYPX(Y, -1, X));
952: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
953: PetscCall(VecDestroy(&Y));
954: }
956: if (view_final) {
957: PetscViewer viewer;
958: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
959: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
960: PetscCall(VecView(X, viewer));
961: PetscCall(PetscViewerPopFormat(viewer));
962: PetscCall(PetscViewerDestroy(&viewer));
963: }
965: /* Clean up */
966: PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
967: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
968: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
969: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
970: PetscCall(VecDestroy(&X));
971: PetscCall(VecDestroy(&X0));
972: PetscCall(VecDestroy(&R));
973: PetscCall(DMDestroy(&da));
974: PetscCall(TSDestroy(&ts));
975: PetscCall(ISDestroy(&ctx.iss));
976: PetscCall(ISDestroy(&ctx.isf));
977: PetscCall(ISDestroy(&ctx.issb));
978: PetscCall(PetscFree(index_slow));
979: PetscCall(PetscFree(index_fast));
980: PetscCall(PetscFree(index_slowbuffer));
981: PetscCall(PetscFunctionListDestroy(&limiters));
982: PetscCall(PetscFunctionListDestroy(&physics));
983: PetscCall(PetscFinalize());
984: return 0;
985: }
987: /*TEST
989: build:
990: requires: !complex
991: depends: finitevolume1d.c
993: test:
994: suffix: 1
995: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22
997: test:
998: suffix: 2
999: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -limit mc -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
1000: output_file: output/ex6_1.out
1002: TEST*/