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;
210: PetscScalar *x, *f, *slope;
211: Vec Xloc;
212: DM da;
214: PetscFunctionBeginUser;
215: ctx->cfl_idt = 0;
216: PetscCall(TSGetDM(ts, &da));
217: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
218: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */
219: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
220: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
221: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
222: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
224: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
226: PetscCall(DMDAVecGetArray(da, Xloc, &x));
227: PetscCall(DMDAVecGetArray(da, F, &f));
228: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
230: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
232: if (ctx->bctype == FVBC_OUTFLOW) {
233: for (i = xs - 2; i < 0; i++) {
234: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
235: }
236: for (i = Mx; i < xs + xm + 2; i++) {
237: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
238: }
239: }
240: for (i = xs - 1; i < xs + xm + 1; i++) {
241: struct _LimitInfo info;
242: PetscScalar *cjmpL, *cjmpR;
243: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
244: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
245: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
246: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
247: cjmpL = &ctx->cjmpLR[0];
248: cjmpR = &ctx->cjmpLR[dof];
249: for (j = 0; j < dof; j++) {
250: PetscScalar jmpL, jmpR;
251: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
252: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
253: for (k = 0; k < dof; k++) {
254: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
255: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
256: }
257: }
258: /* Apply limiter to the left and right characteristic jumps */
259: info.m = dof;
260: info.hxs = hxs;
261: info.hxf = hxf;
262: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
263: for (j = 0; j < dof; j++) {
264: PetscScalar tmp = 0;
265: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
266: slope[i * dof + j] = tmp;
267: }
268: }
270: for (i = xs; i < xs + xm + 1; i++) {
271: PetscReal maxspeed;
272: PetscScalar *uL, *uR;
273: uL = &ctx->uLR[0];
274: uR = &ctx->uLR[dof];
275: if (i < sf) { /* slow region */
276: for (j = 0; j < dof; j++) {
277: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
278: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
279: }
280: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
281: if (i > xs) {
282: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
283: }
284: if (i < xs + xm) {
285: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
286: }
287: } else if (i == sf) { /* interface between the slow region and the fast region */
288: for (j = 0; j < dof; j++) {
289: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
290: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
291: }
292: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
293: if (i > xs) {
294: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
295: }
296: if (i < xs + xm) {
297: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
298: }
299: } else if (i > sf && i < fs) { /* fast region */
300: for (j = 0; j < dof; j++) {
301: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
302: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
303: }
304: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
305: if (i > xs) {
306: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
307: }
308: if (i < xs + xm) {
309: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxf;
310: }
311: } else if (i == fs) { /* interface between the fast region and the slow region */
312: for (j = 0; j < dof; j++) {
313: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
314: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
315: }
316: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
317: if (i > xs) {
318: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxf;
319: }
320: if (i < xs + xm) {
321: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
322: }
323: } else { /* slow region */
324: for (j = 0; j < dof; j++) {
325: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
326: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
327: }
328: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
329: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
330: if (i > xs) {
331: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hxs;
332: }
333: if (i < xs + xm) {
334: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hxs;
335: }
336: }
337: }
338: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
339: PetscCall(DMDAVecRestoreArray(da, F, &f));
340: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
341: PetscCall(DMRestoreLocalVector(da, &Xloc));
342: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
343: if (0) {
344: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
345: PetscReal dt, tnow;
346: PetscCall(TSGetTimeStep(ts, &dt));
347: PetscCall(TSGetTime(ts, &tnow));
348: if (dt > 0.5 / ctx->cfl_idt) {
349: 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))));
350: else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
351: }
352: }
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
356: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
357: PetscErrorCode FVRHSFunctionslow_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
358: {
359: FVCtx *ctx = (FVCtx *)vctx;
360: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
361: PetscReal hxs, hxf;
362: PetscScalar *x, *f, *slope;
363: Vec Xloc;
364: DM da;
366: PetscFunctionBeginUser;
367: ctx->cfl_idt = 0;
368: PetscCall(TSGetDM(ts, &da));
369: PetscCall(DMGetLocalVector(da, &Xloc));
370: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
371: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
372: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
373: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
374: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
375: PetscCall(VecZeroEntries(F));
376: PetscCall(DMDAVecGetArray(da, Xloc, &x));
377: PetscCall(VecGetArray(F, &f));
378: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
379: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
381: if (ctx->bctype == FVBC_OUTFLOW) {
382: for (i = xs - 2; i < 0; i++) {
383: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
384: }
385: for (i = Mx; i < xs + xm + 2; i++) {
386: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
387: }
388: }
389: for (i = xs - 1; i < xs + xm + 1; i++) {
390: struct _LimitInfo info;
391: PetscScalar *cjmpL, *cjmpR;
392: if (i < sf - lsbwidth + 1 || i > fs + rsbwidth - 2) { /* slow components and the first and last fast components */
393: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
394: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
395: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
396: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
397: cjmpL = &ctx->cjmpLR[0];
398: cjmpR = &ctx->cjmpLR[dof];
399: for (j = 0; j < dof; j++) {
400: PetscScalar jmpL, jmpR;
401: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
402: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
403: for (k = 0; k < dof; k++) {
404: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
405: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
406: }
407: }
408: /* Apply limiter to the left and right characteristic jumps */
409: info.m = dof;
410: info.hxs = hxs;
411: info.hxf = hxf;
412: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
413: for (j = 0; j < dof; j++) {
414: PetscScalar tmp = 0;
415: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
416: slope[i * dof + j] = tmp;
417: }
418: }
419: }
421: for (i = xs; i < xs + xm + 1; i++) {
422: PetscReal maxspeed;
423: PetscScalar *uL, *uR;
424: uL = &ctx->uLR[0];
425: uR = &ctx->uLR[dof];
426: if (i < sf - lsbwidth) { /* slow region */
427: for (j = 0; j < dof; j++) {
428: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
429: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
430: }
431: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
432: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hxs)); /* Max allowable value of 1/Delta t */
433: if (i > xs) {
434: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
435: }
436: if (i < xs + xm) {
437: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
438: islow++;
439: }
440: }
441: if (i == sf - lsbwidth) { /* interface between the slow region and the fast region */
442: for (j = 0; j < dof; j++) {
443: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
444: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
445: }
446: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
447: if (i > xs) {
448: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
449: }
450: }
451: if (i == fs + rsbwidth) { /* slow region */
452: for (j = 0; j < dof; j++) {
453: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
454: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
455: }
456: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
457: if (i < xs + xm) {
458: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
459: islow++;
460: }
461: }
462: if (i > fs + rsbwidth) { /* slow region */
463: for (j = 0; j < dof; j++) {
464: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
465: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
466: }
467: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
468: if (i > xs) {
469: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
470: }
471: if (i < xs + xm) {
472: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
473: islow++;
474: }
475: }
476: }
477: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
478: PetscCall(VecRestoreArray(F, &f));
479: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
480: PetscCall(DMRestoreLocalVector(da, &Xloc));
481: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
482: PetscFunctionReturn(PETSC_SUCCESS);
483: }
485: PetscErrorCode FVRHSFunctionslowbuffer_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
486: {
487: FVCtx *ctx = (FVCtx *)vctx;
488: PetscInt i, j, k, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs, lsbwidth = ctx->lsbwidth, rsbwidth = ctx->rsbwidth;
489: PetscReal hxs, hxf;
490: PetscScalar *x, *f, *slope;
491: Vec Xloc;
492: DM da;
494: PetscFunctionBeginUser;
495: PetscCall(TSGetDM(ts, &da));
496: PetscCall(DMGetLocalVector(da, &Xloc));
497: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
498: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
499: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
500: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
501: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
502: PetscCall(VecZeroEntries(F));
503: PetscCall(DMDAVecGetArray(da, Xloc, &x));
504: PetscCall(VecGetArray(F, &f));
505: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
506: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
508: if (ctx->bctype == FVBC_OUTFLOW) {
509: for (i = xs - 2; i < 0; i++) {
510: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
511: }
512: for (i = Mx; i < xs + xm + 2; i++) {
513: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
514: }
515: }
516: for (i = xs - 1; i < xs + xm + 1; i++) {
517: struct _LimitInfo info;
518: PetscScalar *cjmpL, *cjmpR;
519: if ((i > sf - lsbwidth - 2 && i < sf + 1) || (i > fs - 2 && i < fs + rsbwidth + 1)) {
520: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
521: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
522: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
523: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
524: cjmpL = &ctx->cjmpLR[0];
525: cjmpR = &ctx->cjmpLR[dof];
526: for (j = 0; j < dof; j++) {
527: PetscScalar jmpL, jmpR;
528: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
529: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
530: for (k = 0; k < dof; k++) {
531: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
532: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
533: }
534: }
535: /* Apply limiter to the left and right characteristic jumps */
536: info.m = dof;
537: info.hxs = hxs;
538: info.hxf = hxf;
539: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
540: for (j = 0; j < dof; j++) {
541: PetscScalar tmp = 0;
542: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
543: slope[i * dof + j] = tmp;
544: }
545: }
546: }
548: for (i = xs; i < xs + xm + 1; i++) {
549: PetscReal maxspeed;
550: PetscScalar *uL, *uR;
551: uL = &ctx->uLR[0];
552: uR = &ctx->uLR[dof];
553: if (i == sf - lsbwidth) {
554: for (j = 0; j < dof; j++) {
555: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
556: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
557: }
558: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
559: if (i < xs + xm) {
560: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
561: islow++;
562: }
563: }
564: if (i > sf - lsbwidth && i < sf) {
565: for (j = 0; j < dof; j++) {
566: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
567: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
568: }
569: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
570: if (i > xs) {
571: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
572: }
573: if (i < xs + xm) {
574: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
575: islow++;
576: }
577: }
578: if (i == sf) { /* interface between the slow region and the fast region */
579: for (j = 0; j < dof; j++) {
580: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
581: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
582: }
583: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
584: if (i > xs) {
585: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
586: }
587: }
588: if (i == fs) { /* interface between the fast region and the slow region */
589: for (j = 0; j < dof; j++) {
590: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
591: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
592: }
593: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
594: if (i < xs + xm) {
595: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
596: islow++;
597: }
598: }
599: if (i > fs && i < fs + rsbwidth) {
600: for (j = 0; j < dof; j++) {
601: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
602: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
603: }
604: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
605: if (i > xs) {
606: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
607: }
608: if (i < xs + xm) {
609: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hxs;
610: islow++;
611: }
612: }
613: if (i == fs + rsbwidth) {
614: for (j = 0; j < dof; j++) {
615: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
616: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
617: }
618: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
619: if (i > xs) {
620: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hxs;
621: }
622: }
623: }
624: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
625: PetscCall(VecRestoreArray(F, &f));
626: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
627: PetscCall(DMRestoreLocalVector(da, &Xloc));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
632: PetscErrorCode FVRHSFunctionfast_2WaySplit(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
633: {
634: FVCtx *ctx = (FVCtx *)vctx;
635: PetscInt i, j, k, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
636: PetscReal hxs, hxf;
637: PetscScalar *x, *f, *slope;
638: Vec Xloc;
639: DM da;
641: PetscFunctionBeginUser;
642: PetscCall(TSGetDM(ts, &da));
643: PetscCall(DMGetLocalVector(da, &Xloc));
644: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
645: hxs = (ctx->xmax - ctx->xmin) * 3.0 / 8.0 / ctx->sf;
646: hxf = (ctx->xmax - ctx->xmin) / 4.0 / (ctx->fs - ctx->sf);
647: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
648: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
649: PetscCall(VecZeroEntries(F));
650: PetscCall(DMDAVecGetArray(da, Xloc, &x));
651: PetscCall(VecGetArray(F, &f));
652: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
653: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
655: if (ctx->bctype == FVBC_OUTFLOW) {
656: for (i = xs - 2; i < 0; i++) {
657: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
658: }
659: for (i = Mx; i < xs + xm + 2; i++) {
660: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
661: }
662: }
663: for (i = xs - 1; i < xs + xm + 1; i++) {
664: struct _LimitInfo info;
665: PetscScalar *cjmpL, *cjmpR;
666: if (i > sf - 2 && i < fs + 1) {
667: PetscCall((*ctx->physics2.characteristic2)(ctx->physics2.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds));
668: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
669: cjmpL = &ctx->cjmpLR[0];
670: cjmpR = &ctx->cjmpLR[dof];
671: for (j = 0; j < dof; j++) {
672: PetscScalar jmpL, jmpR;
673: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
674: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
675: for (k = 0; k < dof; k++) {
676: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
677: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
678: }
679: }
680: /* Apply limiter to the left and right characteristic jumps */
681: info.m = dof;
682: info.hxs = hxs;
683: info.hxf = hxf;
684: (*ctx->limit2)(&info, cjmpL, cjmpR, ctx->sf, ctx->fs, i, ctx->cslope);
685: for (j = 0; j < dof; j++) {
686: PetscScalar tmp = 0;
687: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
688: slope[i * dof + j] = tmp;
689: }
690: }
691: }
693: for (i = xs; i < xs + xm + 1; i++) {
694: PetscReal maxspeed;
695: PetscScalar *uL, *uR;
696: uL = &ctx->uLR[0];
697: uR = &ctx->uLR[dof];
698: if (i == sf) { /* interface between the slow region and the fast region */
699: for (j = 0; j < dof; j++) {
700: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxs / 2;
701: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
702: }
703: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
704: if (i < xs + xm) {
705: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
706: ifast++;
707: }
708: }
709: if (i > sf && i < fs) { /* fast region */
710: for (j = 0; j < dof; j++) {
711: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
712: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxf / 2;
713: }
714: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
715: if (i > xs) {
716: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
717: }
718: if (i < xs + xm) {
719: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hxf;
720: ifast++;
721: }
722: }
723: if (i == fs) { /* interface between the fast region and the slow region */
724: for (j = 0; j < dof; j++) {
725: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hxf / 2;
726: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hxs / 2;
727: }
728: PetscCall((*ctx->physics2.riemann2)(ctx->physics2.user, dof, uL, uR, ctx->flux, &maxspeed));
729: if (i > xs) {
730: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hxf;
731: }
732: }
733: }
734: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
735: PetscCall(VecRestoreArray(F, &f));
736: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
737: PetscCall(DMRestoreLocalVector(da, &Xloc));
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: int main(int argc, char *argv[])
742: {
743: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
744: PetscFunctionList limiters = 0, physics = 0;
745: MPI_Comm comm;
746: TS ts;
747: DM da;
748: Vec X, X0, R;
749: FVCtx ctx;
750: 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;
751: PetscBool view_final = PETSC_FALSE;
752: PetscReal ptime;
754: PetscFunctionBeginUser;
755: PetscCall(PetscInitialize(&argc, &argv, 0, help));
756: comm = PETSC_COMM_WORLD;
757: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
759: /* Register limiters to be available on the command line */
760: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit2_Upwind));
761: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit2_LaxWendroff));
762: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit2_BeamWarming));
763: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit2_Fromm));
764: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit2_Minmod));
765: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit2_Superbee));
766: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit2_MC));
767: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit2_Koren3));
769: /* Register physical models to be available on the command line */
770: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
772: ctx.comm = comm;
773: ctx.cfl = 0.9;
774: ctx.bctype = FVBC_PERIODIC;
775: ctx.xmin = -1.0;
776: ctx.xmax = 1.0;
777: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
778: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
779: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
780: PetscCall(PetscOptionsFList("-limit", "Name of flux imiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
781: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
782: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
783: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
784: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
785: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
786: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
787: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
788: PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
789: PetscOptionsEnd();
791: /* Choose the limiter from the list of registered limiters */
792: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit2));
793: PetscCheck(ctx.limit2, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
795: /* Choose the physics from the list of registered models */
796: {
797: PetscErrorCode (*r)(FVCtx *);
798: PetscCall(PetscFunctionListFind(physics, physname, &r));
799: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
800: /* Create the physics, will set the number of fields and their names */
801: PetscCall((*r)(&ctx));
802: }
804: /* Create a DMDA to manage the parallel grid */
805: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics2.dof, 2, NULL, &da));
806: PetscCall(DMSetFromOptions(da));
807: PetscCall(DMSetUp(da));
808: /* Inform the DMDA of the field names provided by the physics. */
809: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
810: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics2.fieldname[i]));
811: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
812: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
814: /* Set coordinates of cell centers */
815: 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));
817: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
818: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
819: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
821: /* Create a vector to store the solution and to save the initial state */
822: PetscCall(DMCreateGlobalVector(da, &X));
823: PetscCall(VecDuplicate(X, &X0));
824: PetscCall(VecDuplicate(X, &R));
826: /* create index for slow parts and fast parts,
827: count_slow + count_fast = Mx, counts_slow*hs = 0.5, counts_fast*hf = 0.5 */
828: count_slow = Mx / (1.0 + ctx.hratio / 3.0);
829: 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");
830: count_fast = Mx - count_slow;
831: ctx.sf = count_slow / 2;
832: ctx.fs = ctx.sf + count_fast;
833: PetscCall(PetscMalloc1(xm * dof, &index_slow));
834: PetscCall(PetscMalloc1(xm * dof, &index_fast));
835: PetscCall(PetscMalloc1(6 * dof, &index_slowbuffer));
836: if (((AdvectCtx *)ctx.physics2.user)->a > 0) {
837: ctx.lsbwidth = 2;
838: ctx.rsbwidth = 4;
839: } else {
840: ctx.lsbwidth = 4;
841: ctx.rsbwidth = 2;
842: }
843: for (i = xs; i < xs + xm; i++) {
844: if (i < ctx.sf - ctx.lsbwidth || i > ctx.fs + ctx.rsbwidth - 1)
845: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
846: else if ((i >= ctx.sf - ctx.lsbwidth && i < ctx.sf) || (i > ctx.fs - 1 && i <= ctx.fs + ctx.rsbwidth - 1))
847: for (k = 0; k < dof; k++) index_slowbuffer[islowbuffer++] = i * dof + k;
848: else
849: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
850: }
851: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
852: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
853: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islowbuffer, index_slowbuffer, PETSC_COPY_VALUES, &ctx.issb));
855: /* Create a time-stepping object */
856: PetscCall(TSCreate(comm, &ts));
857: PetscCall(TSSetDM(ts, da));
858: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction_2WaySplit, &ctx));
859: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
860: PetscCall(TSRHSSplitSetIS(ts, "slowbuffer", ctx.issb));
861: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
862: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow_2WaySplit, &ctx));
863: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast_2WaySplit, &ctx));
864: PetscCall(TSRHSSplitSetRHSFunction(ts, "slowbuffer", NULL, FVRHSFunctionslowbuffer_2WaySplit, &ctx));
866: PetscCall(TSSetType(ts, TSSSP));
867: /*PetscCall(TSSetType(ts,TSMPRK));*/
868: PetscCall(TSSetMaxTime(ts, 10));
869: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
871: /* Compute initial conditions and starting time step */
872: PetscCall(FVSample_2WaySplit(&ctx, da, 0, X0));
873: PetscCall(FVRHSFunction_2WaySplit(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
874: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
875: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
876: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
877: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
878: {
879: PetscInt steps;
880: PetscScalar mass_initial, mass_final, mass_difference;
881: const PetscScalar *ptr_X, *ptr_X0;
882: const PetscReal hs = (ctx.xmax - ctx.xmin) * 3.0 / 4.0 / count_slow;
883: const PetscReal hf = (ctx.xmax - ctx.xmin) / 4.0 / count_fast;
885: PetscCall(TSSolve(ts, X));
886: PetscCall(TSGetSolveTime(ts, &ptime));
887: PetscCall(TSGetStepNumber(ts, &steps));
888: /* calculate the total mass at initial time and final time */
889: mass_initial = 0.0;
890: mass_final = 0.0;
891: PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
892: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
893: for (i = xs; i < xs + xm; i++) {
894: if (i < ctx.sf || i > ctx.fs - 1) {
895: for (k = 0; k < dof; k++) {
896: mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
897: mass_final = mass_final + hs * ptr_X[i * dof + k];
898: }
899: } else {
900: for (k = 0; k < dof; k++) {
901: mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
902: mass_final = mass_final + hf * ptr_X[i * dof + k];
903: }
904: }
905: }
906: PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
907: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
908: mass_difference = mass_final - mass_initial;
909: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mass_difference, 1, MPIU_SCALAR, MPIU_SUM, comm));
910: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
911: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
912: PetscCall(PetscPrintf(comm, "Maximum allowable stepsize according to CFL %g\n", (double)(1 / ctx.cfl_idt)));
913: if (ctx.exact) {
914: PetscReal nrm1 = 0;
915: PetscCall(SolutionErrorNorms_2WaySplit(&ctx, da, ptime, X, &nrm1));
916: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
917: }
918: if (ctx.simulation) {
919: PetscReal nrm1 = 0;
920: PetscViewer fd;
921: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
922: Vec XR;
923: PetscBool flg;
924: const PetscScalar *ptr_XR;
925: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
926: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
927: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
928: PetscCall(VecDuplicate(X0, &XR));
929: PetscCall(VecLoad(XR, fd));
930: PetscCall(PetscViewerDestroy(&fd));
931: PetscCall(VecGetArrayRead(X, &ptr_X));
932: PetscCall(VecGetArrayRead(XR, &ptr_XR));
933: for (i = xs; i < xs + xm; i++) {
934: if (i < ctx.sf || i > ctx.fs - 1)
935: for (k = 0; k < dof; k++) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
936: else
937: for (k = 0; k < dof; k++) nrm1 = nrm1 + hf * PetscAbs(ptr_X[i * dof + k] - ptr_XR[i * dof + k]);
938: }
939: PetscCall(VecRestoreArrayRead(X, &ptr_X));
940: PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
941: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
942: PetscCall(VecDestroy(&XR));
943: }
944: }
946: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
947: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
948: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
949: if (draw & 0x4) {
950: Vec Y;
951: PetscCall(VecDuplicate(X, &Y));
952: PetscCall(FVSample_2WaySplit(&ctx, da, ptime, Y));
953: PetscCall(VecAYPX(Y, -1, X));
954: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
955: PetscCall(VecDestroy(&Y));
956: }
958: if (view_final) {
959: PetscViewer viewer;
960: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
961: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
962: PetscCall(VecView(X, viewer));
963: PetscCall(PetscViewerPopFormat(viewer));
964: PetscCall(PetscViewerDestroy(&viewer));
965: }
967: /* Clean up */
968: PetscCall((*ctx.physics2.destroy)(ctx.physics2.user));
969: for (i = 0; i < ctx.physics2.dof; i++) PetscCall(PetscFree(ctx.physics2.fieldname[i]));
970: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
971: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
972: PetscCall(VecDestroy(&X));
973: PetscCall(VecDestroy(&X0));
974: PetscCall(VecDestroy(&R));
975: PetscCall(DMDestroy(&da));
976: PetscCall(TSDestroy(&ts));
977: PetscCall(ISDestroy(&ctx.iss));
978: PetscCall(ISDestroy(&ctx.isf));
979: PetscCall(ISDestroy(&ctx.issb));
980: PetscCall(PetscFree(index_slow));
981: PetscCall(PetscFree(index_fast));
982: PetscCall(PetscFree(index_slowbuffer));
983: PetscCall(PetscFunctionListDestroy(&limiters));
984: PetscCall(PetscFunctionListDestroy(&physics));
985: PetscCall(PetscFinalize());
986: return 0;
987: }
989: /*TEST
991: build:
992: requires: !complex
993: depends: finitevolume1d.c
995: test:
996: suffix: 1
997: args: -da_grid_x 60 -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 2a22
999: test:
1000: suffix: 2
1001: args: -da_grid_x 60 -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 2a22 -ts_use_splitrhsfunction 0
1002: output_file: output/ex6_1.out
1004: TEST*/