Actual source code: ex7.c
1: static const char help[] = "1D periodic Finite Volume solver by a particular slope limiter 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, say\n"
5: " hxs = (xmax - xmin)/2.0*(hratio+1.0)/Mx, \n"
6: " hxf = (xmax - xmin)/2.0*(1.0+1.0/hratio)/Mx, \n"
7: " with x belongs to (xmin,xmax), the number of total mesh points is Mx and the ratio between the meshsize of coarse\n\n"
8: " grids and fine grids is hratio.\n"
9: " exact - Exact Riemann solver which usually needs to perform a Newton iteration to connect\n"
10: " the states across shocks and rarefactions\n"
11: " simulation - use reference solution which is generated by smaller time step size to be true solution,\n"
12: " also the reference solution should be generated by user and stored in a binary file.\n"
13: " characteristic - Limit the characteristic variables, this is usually preferred (default)\n"
14: "Several initial conditions can be chosen with -initial N\n\n"
15: "The problem size should be set with -da_grid_x M\n\n"
16: "This script choose the slope limiter by biased second-order upwind procedure which is proposed by Van Leer in 1994\n"
17: " u(x_(k+1/2),t) = u(x_k,t) + phi(x_(k+1/2),t)*(u(x_k,t)-u(x_(k-1),t)) \n"
18: " limiter phi(x_(k+1/2),t) = max(0,min(r(k+1/2),min(2,gamma(k+1/2)*r(k+1/2)+alpha(k+1/2)))) \n"
19: " r(k+1/2) = (u(x_(k+1))-u(x_k))/(u(x_k)-u(x_(k-1))) \n"
20: " alpha(k+1/2) = (h_k*h_(k+1))/(h_(k-1)+h_k)/(h_(k-1)+h_k+h_(k+1)) \n"
21: " gamma(k+1/2) = h_k*(h_(k-1)+h_k)/(h_k+h_(k+1))/(h_(k-1)+h_k+h_(k+1)) \n";
23: #include <petscts.h>
24: #include <petscdm.h>
25: #include <petscdmda.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: /* --------------------------------- Finite Volume data structures ----------------------------------- */
35: typedef enum {
36: FVBC_PERIODIC,
37: FVBC_OUTFLOW
38: } FVBCType;
39: static const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "FVBCType", "FVBC_", NULL};
41: typedef struct {
42: PetscErrorCode (*sample)(void *, PetscInt, FVBCType, PetscReal, PetscReal, PetscReal, PetscReal, PetscReal *);
43: PetscErrorCode (*flux)(void *, const PetscScalar *, PetscScalar *, PetscReal *);
44: PetscErrorCode (*destroy)(void *);
45: void *user;
46: PetscInt dof;
47: char *fieldname[16];
48: } PhysicsCtx;
50: typedef struct {
51: PhysicsCtx physics;
52: MPI_Comm comm;
53: char prefix[256];
55: /* Local work arrays */
56: PetscScalar *flux; /* Flux across interface */
57: PetscReal *speeds; /* Speeds of each wave */
58: PetscScalar *u; /* value at face */
60: PetscReal cfl_idt; /* Max allowable value of 1/Delta t */
61: PetscReal cfl;
62: PetscReal xmin, xmax;
63: PetscInt initial;
64: PetscBool exact;
65: PetscBool simulation;
66: FVBCType bctype;
67: PetscInt hratio; /* hratio = hslow/hfast */
68: IS isf, iss;
69: PetscInt sf, fs; /* slow-fast and fast-slow interfaces */
70: } FVCtx;
72: /* --------------------------------- Physics ----------------------------------- */
73: static PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
74: {
75: PetscFunctionBeginUser;
76: PetscCall(PetscFree(vctx));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: /* --------------------------------- Advection ----------------------------------- */
81: typedef struct {
82: PetscReal a; /* advective velocity */
83: } AdvectCtx;
85: static PetscErrorCode PhysicsFlux_Advect(void *vctx, const PetscScalar *u, PetscScalar *flux, PetscReal *maxspeed)
86: {
87: AdvectCtx *ctx = (AdvectCtx *)vctx;
88: PetscReal speed;
90: PetscFunctionBeginUser;
91: speed = ctx->a;
92: flux[0] = speed * u[0];
93: *maxspeed = speed;
94: PetscFunctionReturn(PETSC_SUCCESS);
95: }
97: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
98: {
99: AdvectCtx *ctx = (AdvectCtx *)vctx;
100: PetscReal a = ctx->a, x0;
102: PetscFunctionBeginUser;
103: switch (bctype) {
104: case FVBC_OUTFLOW:
105: x0 = x - a * t;
106: break;
107: case FVBC_PERIODIC:
108: x0 = RangeMod(x - a * t, xmin, xmax);
109: break;
110: default:
111: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
112: }
113: switch (initial) {
114: case 0:
115: u[0] = (x0 < 0) ? 1 : -1;
116: break;
117: case 1:
118: u[0] = (x0 < 0) ? -1 : 1;
119: break;
120: case 2:
121: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
122: break;
123: case 3:
124: u[0] = PetscSinReal(2 * PETSC_PI * x0);
125: break;
126: case 4:
127: u[0] = PetscAbs(x0);
128: break;
129: case 5:
130: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
131: break;
132: case 6:
133: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
134: break;
135: case 7:
136: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
137: break;
138: default:
139: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
140: }
141: PetscFunctionReturn(PETSC_SUCCESS);
142: }
144: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
145: {
146: AdvectCtx *user;
148: PetscFunctionBeginUser;
149: PetscCall(PetscNew(&user));
150: ctx->physics.sample = PhysicsSample_Advect;
151: ctx->physics.flux = PhysicsFlux_Advect;
152: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
153: ctx->physics.user = user;
154: ctx->physics.dof = 1;
155: PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
156: user->a = 1;
157: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
158: {
159: PetscCall(PetscOptionsReal("-physics_advect_a", "Speed", "", user->a, &user->a, NULL));
160: }
161: PetscOptionsEnd();
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: /* --------------------------------- Finite Volume Solver ----------------------------------- */
167: static PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
168: {
169: FVCtx *ctx = (FVCtx *)vctx;
170: PetscInt i, j, Mx, dof, xs, xm, sf = ctx->sf, fs = ctx->fs;
171: PetscReal hf, hs;
172: PetscScalar *x, *f, *r, *min, *alpha, *gamma;
173: Vec Xloc;
174: DM da;
176: PetscFunctionBeginUser;
177: ctx->cfl_idt = 0;
178: PetscCall(TSGetDM(ts, &da));
179: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
180: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); /* Mx is the number of center points */
181: hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
182: hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
183: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
184: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
185: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
186: PetscCall(DMDAVecGetArray(da, Xloc, &x));
187: PetscCall(DMDAVecGetArray(da, F, &f));
188: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
189: PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));
191: if (ctx->bctype == FVBC_OUTFLOW) {
192: for (i = xs - 2; i < 0; i++) {
193: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
194: }
195: for (i = Mx; i < xs + xm + 2; i++) {
196: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
197: }
198: }
200: for (i = xs; i < xs + xm + 1; i++) {
201: PetscReal maxspeed;
202: PetscScalar *u;
203: if (i < sf || i > fs + 1) {
204: u = &ctx->u[0];
205: alpha[0] = 1.0 / 6.0;
206: gamma[0] = 1.0 / 3.0;
207: for (j = 0; j < dof; j++) {
208: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
209: min[j] = PetscMin(r[j], 2.0);
210: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
211: }
212: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
213: ctx->cfl_idt = PetscMax(ctx->cfl_idt, PetscAbsScalar(maxspeed / hs));
214: if (i > xs) {
215: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
216: }
217: if (i < xs + xm) {
218: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
219: }
220: } else if (i == sf) {
221: u = &ctx->u[0];
222: alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
223: gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
224: for (j = 0; j < dof; j++) {
225: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
226: min[j] = PetscMin(r[j], 2.0);
227: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
228: }
229: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
230: if (i > xs) {
231: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
232: }
233: if (i < xs + xm) {
234: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
235: }
236: } else if (i == sf + 1) {
237: u = &ctx->u[0];
238: alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf);
239: gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf);
240: for (j = 0; j < dof; j++) {
241: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
242: min[j] = PetscMin(r[j], 2.0);
243: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
244: }
245: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
246: if (i > xs) {
247: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
248: }
249: if (i < xs + xm) {
250: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
251: }
252: } else if (i > sf + 1 && i < fs) {
253: u = &ctx->u[0];
254: alpha[0] = 1.0 / 6.0;
255: gamma[0] = 1.0 / 3.0;
256: for (j = 0; j < dof; j++) {
257: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
258: min[j] = PetscMin(r[j], 2.0);
259: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
260: }
261: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
262: if (i > xs) {
263: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
264: }
265: if (i < xs + xm) {
266: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hf;
267: }
268: } else if (i == fs) {
269: u = &ctx->u[0];
270: alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
271: gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
272: for (j = 0; j < dof; j++) {
273: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
274: min[j] = PetscMin(r[j], 2.0);
275: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
276: }
277: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
278: if (i > xs) {
279: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hf;
280: }
281: if (i < xs + xm) {
282: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
283: }
284: } else if (i == fs + 1) {
285: u = &ctx->u[0];
286: alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs);
287: gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs);
288: for (j = 0; j < dof; j++) {
289: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
290: min[j] = PetscMin(r[j], 2.0);
291: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
292: }
293: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
294: if (i > xs) {
295: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hs;
296: }
297: if (i < xs + xm) {
298: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hs;
299: }
300: }
301: }
302: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
303: PetscCall(DMDAVecRestoreArray(da, F, &f));
304: PetscCall(DMRestoreLocalVector(da, &Xloc));
305: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &ctx->cfl_idt, 1, MPIU_SCALAR, MPIU_MAX, PetscObjectComm((PetscObject)da)));
306: if (0) {
307: /* We need a way to inform the TS of a CFL constraint, this is a debugging fragment */
308: PetscReal dt, tnow;
309: PetscCall(TSGetTimeStep(ts, &dt));
310: PetscCall(TSGetTime(ts, &tnow));
311: 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))));
312: }
313: PetscCall(PetscFree4(r, min, alpha, gamma));
314: PetscFunctionReturn(PETSC_SUCCESS);
315: }
317: static PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
318: {
319: FVCtx *ctx = (FVCtx *)vctx;
320: PetscInt i, j, Mx, dof, xs, xm, islow = 0, sf = ctx->sf, fs = ctx->fs;
321: PetscReal hf, hs;
322: PetscScalar *x, *f, *r, *min, *alpha, *gamma;
323: Vec Xloc;
324: DM da;
326: PetscFunctionBeginUser;
327: PetscCall(TSGetDM(ts, &da));
328: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
329: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); /* Mx is the number of center points */
330: hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
331: hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
332: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
333: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
334: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
335: PetscCall(DMDAVecGetArray(da, Xloc, &x));
336: PetscCall(VecGetArray(F, &f));
337: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
338: PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));
340: if (ctx->bctype == FVBC_OUTFLOW) {
341: for (i = xs - 2; i < 0; i++) {
342: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
343: }
344: for (i = Mx; i < xs + xm + 2; i++) {
345: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
346: }
347: }
349: for (i = xs; i < xs + xm + 1; i++) {
350: PetscReal maxspeed;
351: PetscScalar *u;
352: if (i < sf) {
353: u = &ctx->u[0];
354: alpha[0] = 1.0 / 6.0;
355: gamma[0] = 1.0 / 3.0;
356: for (j = 0; j < dof; j++) {
357: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
358: min[j] = PetscMin(r[j], 2.0);
359: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
360: }
361: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
362: if (i > xs) {
363: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
364: }
365: if (i < xs + xm) {
366: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
367: islow++;
368: }
369: } else if (i == sf) {
370: u = &ctx->u[0];
371: alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
372: gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
373: for (j = 0; j < dof; j++) {
374: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
375: min[j] = PetscMin(r[j], 2.0);
376: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
377: }
378: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
379: if (i > xs) {
380: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
381: }
382: } else if (i == fs) {
383: u = &ctx->u[0];
384: alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
385: gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
386: for (j = 0; j < dof; j++) {
387: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
388: min[j] = PetscMin(r[j], 2.0);
389: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
390: }
391: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
392: if (i < xs + xm) {
393: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
394: islow++;
395: }
396: } else if (i == fs + 1) {
397: u = &ctx->u[0];
398: alpha[0] = hs * hs / (hf + hs) / (hf + hs + hs);
399: gamma[0] = hs * (hf + hs) / (hs + hs) / (hf + hs + hs);
400: for (j = 0; j < dof; j++) {
401: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
402: min[j] = PetscMin(r[j], 2.0);
403: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
404: }
405: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
406: if (i > xs) {
407: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
408: }
409: if (i < xs + xm) {
410: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
411: islow++;
412: }
413: } else if (i > fs + 1) {
414: u = &ctx->u[0];
415: alpha[0] = 1.0 / 6.0;
416: gamma[0] = 1.0 / 3.0;
417: for (j = 0; j < dof; j++) {
418: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
419: min[j] = PetscMin(r[j], 2.0);
420: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
421: }
422: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
423: if (i > xs) {
424: for (j = 0; j < dof; j++) f[(islow - 1) * dof + j] -= ctx->flux[j] / hs;
425: }
426: if (i < xs + xm) {
427: for (j = 0; j < dof; j++) f[islow * dof + j] += ctx->flux[j] / hs;
428: islow++;
429: }
430: }
431: }
432: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
433: PetscCall(VecRestoreArray(F, &f));
434: PetscCall(DMRestoreLocalVector(da, &Xloc));
435: PetscCall(PetscFree4(r, min, alpha, gamma));
436: PetscFunctionReturn(PETSC_SUCCESS);
437: }
439: static PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
440: {
441: FVCtx *ctx = (FVCtx *)vctx;
442: PetscInt i, j, Mx, dof, xs, xm, ifast = 0, sf = ctx->sf, fs = ctx->fs;
443: PetscReal hf, hs;
444: PetscScalar *x, *f, *r, *min, *alpha, *gamma;
445: Vec Xloc;
446: DM da;
448: PetscFunctionBeginUser;
449: PetscCall(TSGetDM(ts, &da));
450: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
451: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL)); /* Mx is the number of center points */
452: hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
453: hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
454: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
455: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
456: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
457: PetscCall(DMDAVecGetArray(da, Xloc, &x));
458: PetscCall(VecGetArray(F, &f));
459: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
460: PetscCall(PetscMalloc4(dof, &r, dof, &min, dof, &alpha, dof, &gamma));
462: if (ctx->bctype == FVBC_OUTFLOW) {
463: for (i = xs - 2; i < 0; i++) {
464: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
465: }
466: for (i = Mx; i < xs + xm + 2; i++) {
467: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
468: }
469: }
471: for (i = xs; i < xs + xm + 1; i++) {
472: PetscReal maxspeed;
473: PetscScalar *u;
474: if (i == sf) {
475: u = &ctx->u[0];
476: alpha[0] = hs * hf / (hs + hs) / (hs + hs + hf);
477: gamma[0] = hs * (hs + hs) / (hs + hf) / (hs + hs + hf);
478: for (j = 0; j < dof; j++) {
479: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
480: min[j] = PetscMin(r[j], 2.0);
481: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
482: }
483: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
484: if (i < xs + xm) {
485: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
486: ifast++;
487: }
488: } else if (i == sf + 1) {
489: u = &ctx->u[0];
490: alpha[0] = hf * hf / (hs + hf) / (hs + hf + hf);
491: gamma[0] = hf * (hs + hf) / (hf + hf) / (hs + hf + hf);
492: for (j = 0; j < dof; j++) {
493: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
494: min[j] = PetscMin(r[j], 2.0);
495: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
496: }
497: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
498: if (i > xs) {
499: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
500: }
501: if (i < xs + xm) {
502: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
503: ifast++;
504: }
505: } else if (i > sf + 1 && i < fs) {
506: u = &ctx->u[0];
507: alpha[0] = 1.0 / 6.0;
508: gamma[0] = 1.0 / 3.0;
509: for (j = 0; j < dof; j++) {
510: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
511: min[j] = PetscMin(r[j], 2.0);
512: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
513: }
514: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
515: if (i > xs) {
516: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
517: }
518: if (i < xs + xm) {
519: for (j = 0; j < dof; j++) f[ifast * dof + j] += ctx->flux[j] / hf;
520: ifast++;
521: }
522: } else if (i == fs) {
523: u = &ctx->u[0];
524: alpha[0] = hf * hs / (hf + hf) / (hf + hf + hs);
525: gamma[0] = hf * (hf + hf) / (hf + hs) / (hf + hf + hs);
526: for (j = 0; j < dof; j++) {
527: r[j] = (x[i * dof + j] - x[(i - 1) * dof + j]) / (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
528: min[j] = PetscMin(r[j], 2.0);
529: u[j] = x[(i - 1) * dof + j] + PetscMax(0, PetscMin(min[j], alpha[0] + gamma[0] * r[j])) * (x[(i - 1) * dof + j] - x[(i - 2) * dof + j]);
530: }
531: PetscCall((*ctx->physics.flux)(ctx->physics.user, u, ctx->flux, &maxspeed));
532: if (i > xs) {
533: for (j = 0; j < dof; j++) f[(ifast - 1) * dof + j] -= ctx->flux[j] / hf;
534: }
535: }
536: }
537: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
538: PetscCall(VecRestoreArray(F, &f));
539: PetscCall(DMRestoreLocalVector(da, &Xloc));
540: PetscCall(PetscFree4(r, min, alpha, gamma));
541: PetscFunctionReturn(PETSC_SUCCESS);
542: }
544: /* --------------------------------- Finite Volume Solver for slow components ----------------------------------- */
546: PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
547: {
548: PetscScalar *u, *uj, xj, xi;
549: PetscInt i, j, k, dof, xs, xm, Mx, count_slow, count_fast;
550: const PetscInt N = 200;
552: PetscFunctionBeginUser;
553: PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
554: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
555: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
556: PetscCall(DMDAVecGetArray(da, U, &u));
557: PetscCall(PetscMalloc1(dof, &uj));
558: const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
559: const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
560: count_slow = Mx / (1 + ctx->hratio);
561: count_fast = Mx - count_slow;
562: for (i = xs; i < xs + xm; i++) {
563: if (i * hs + 0.5 * hs < (ctx->xmax - ctx->xmin) * 0.25) {
564: xi = ctx->xmin + 0.5 * hs + i * hs;
565: /* Integrate over cell i using trapezoid rule with N points. */
566: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
567: for (j = 0; j < N + 1; j++) {
568: xj = xi + hs * (j - N / 2) / (PetscReal)N;
569: PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
570: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
571: }
572: } else if ((ctx->xmax - ctx->xmin) * 0.25 + (i - count_slow / 2) * hf + 0.5 * hf < (ctx->xmax - ctx->xmin) * 0.75) {
573: xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.25 + 0.5 * hf + (i - count_slow / 2) * hf;
574: /* Integrate over cell i using trapezoid rule with N points. */
575: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
576: for (j = 0; j < N + 1; j++) {
577: xj = xi + hf * (j - N / 2) / (PetscReal)N;
578: PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
579: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
580: }
581: } else {
582: xi = ctx->xmin + (ctx->xmax - ctx->xmin) * 0.75 + 0.5 * hs + (i - count_slow / 2 - count_fast) * hs;
583: /* Integrate over cell i using trapezoid rule with N points. */
584: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
585: for (j = 0; j < N + 1; j++) {
586: xj = xi + hs * (j - N / 2) / (PetscReal)N;
587: PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
588: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
589: }
590: }
591: }
592: PetscCall(DMDAVecRestoreArray(da, U, &u));
593: PetscCall(PetscFree(uj));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: static PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
598: {
599: PetscReal xmin, xmax;
600: PetscScalar sum, tvsum;
601: const PetscScalar *x;
602: PetscInt imin, imax, Mx, i, j, xs, xm, dof;
603: Vec Xloc;
604: PetscBool isascii;
606: PetscFunctionBeginUser;
607: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
608: PetscCheck(isascii, PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
609: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
610: PetscCall(DMGetLocalVector(da, &Xloc));
611: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
612: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
613: PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
614: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
615: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
616: tvsum = 0;
617: for (i = xs; i < xs + xm; i++) {
618: for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
619: }
620: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &tvsum, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)da)));
621: PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
622: PetscCall(DMRestoreLocalVector(da, &Xloc));
624: PetscCall(VecMin(X, &imin, &xmin));
625: PetscCall(VecMax(X, &imax, &xmax));
626: PetscCall(VecSum(X, &sum));
627: PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%g,%g] with minimum at %" PetscInt_FMT ", mean %g, ||x||_TV %g\n", (double)xmin, (double)xmax, imin, (double)(sum / Mx), (double)(tvsum / Mx)));
628: PetscFunctionReturn(PETSC_SUCCESS);
629: }
631: static PetscErrorCode SolutionErrorNorms(FVCtx *ctx, DM da, PetscReal t, Vec X, PetscReal *nrm1)
632: {
633: Vec Y;
634: PetscInt i, Mx, count_slow = 0, count_fast = 0;
635: const PetscScalar *ptr_X, *ptr_Y;
637: PetscFunctionBeginUser;
638: PetscCall(VecGetSize(X, &Mx));
639: PetscCall(VecDuplicate(X, &Y));
640: PetscCall(FVSample(ctx, da, t, Y));
641: const PetscReal hs = (ctx->xmax - ctx->xmin) / 2.0 * (ctx->hratio + 1.0) / Mx;
642: const PetscReal hf = (ctx->xmax - ctx->xmin) / 2.0 * (1.0 + 1.0 / ctx->hratio) / Mx;
643: count_slow = (PetscReal)Mx / (1.0 + ctx->hratio);
644: count_fast = Mx - count_slow;
645: PetscCall(VecGetArrayRead(X, &ptr_X));
646: PetscCall(VecGetArrayRead(Y, &ptr_Y));
647: for (i = 0; i < Mx; i++) {
648: if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) *nrm1 += hs * PetscAbs(ptr_X[i] - ptr_Y[i]);
649: else *nrm1 += hf * PetscAbs(ptr_X[i] - ptr_Y[i]);
650: }
651: PetscCall(VecRestoreArrayRead(X, &ptr_X));
652: PetscCall(VecRestoreArrayRead(Y, &ptr_Y));
653: PetscCall(VecDestroy(&Y));
654: PetscFunctionReturn(PETSC_SUCCESS);
655: }
657: int main(int argc, char *argv[])
658: {
659: char physname[256] = "advect", final_fname[256] = "solution.m";
660: PetscFunctionList physics = NULL;
661: MPI_Comm comm;
662: TS ts;
663: DM da;
664: Vec X, X0, R;
665: FVCtx ctx;
666: PetscInt i, k, dof, xs, xm, Mx, draw = 0, count_slow, count_fast, islow = 0, ifast = 0, *index_slow, *index_fast;
667: PetscBool view_final = PETSC_FALSE;
668: PetscReal ptime;
670: PetscFunctionBeginUser;
671: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
672: comm = PETSC_COMM_WORLD;
673: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
675: /* Register physical models to be available on the command line */
676: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
678: ctx.comm = comm;
679: ctx.cfl = 0.9;
680: ctx.bctype = FVBC_PERIODIC;
681: ctx.xmin = -1.0;
682: ctx.xmax = 1.0;
683: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
684: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
685: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
686: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
687: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
688: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
689: PetscCall(PetscOptionsBool("-exact", "Compare errors with exact solution", "", ctx.exact, &ctx.exact, NULL));
690: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
691: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
692: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
693: PetscCall(PetscOptionsInt("-hratio", "Spacing ratio", "", ctx.hratio, &ctx.hratio, NULL));
694: PetscOptionsEnd();
696: /* Choose the physics from the list of registered models */
697: {
698: PetscErrorCode (*r)(FVCtx *);
699: PetscCall(PetscFunctionListFind(physics, physname, &r));
700: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
701: /* Create the physics, will set the number of fields and their names */
702: PetscCall((*r)(&ctx));
703: }
705: /* Create a DMDA to manage the parallel grid */
706: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
707: PetscCall(DMSetFromOptions(da));
708: PetscCall(DMSetUp(da));
709: /* Inform the DMDA of the field names provided by the physics. */
710: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
711: for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
712: PetscCall(DMDAGetInfo(da, NULL, &Mx, NULL, NULL, NULL, NULL, NULL, &dof, NULL, NULL, NULL, NULL, NULL));
713: PetscCall(DMDAGetCorners(da, &xs, NULL, NULL, &xm, NULL, NULL));
715: /* Set coordinates of cell centers */
716: 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));
718: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
719: PetscCall(PetscMalloc3(dof, &ctx.u, dof, &ctx.flux, dof, &ctx.speeds));
721: /* Create a vector to store the solution and to save the initial state */
722: PetscCall(DMCreateGlobalVector(da, &X));
723: PetscCall(VecDuplicate(X, &X0));
724: PetscCall(VecDuplicate(X, &R));
726: /* create index for slow parts and fast parts*/
727: count_slow = Mx / (1 + ctx.hratio);
728: 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) is even");
729: count_fast = Mx - count_slow;
730: ctx.sf = count_slow / 2;
731: ctx.fs = ctx.sf + count_fast;
732: PetscCall(PetscMalloc1(xm * dof, &index_slow));
733: PetscCall(PetscMalloc1(xm * dof, &index_fast));
734: for (i = xs; i < xs + xm; i++) {
735: if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1)
736: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
737: else
738: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
739: }
740: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
741: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
743: /* Create a time-stepping object */
744: PetscCall(TSCreate(comm, &ts));
745: PetscCall(TSSetDM(ts, da));
746: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
747: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
748: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
749: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx));
750: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx));
752: PetscCall(TSSetType(ts, TSMPRK));
753: PetscCall(TSSetMaxTime(ts, 10));
754: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
756: /* Compute initial conditions and starting time step */
757: PetscCall(FVSample(&ctx, da, 0, X0));
758: PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
759: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
760: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
761: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
762: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
763: {
764: PetscInt steps;
765: PetscScalar mass_initial, mass_final, mass_difference;
766: const PetscScalar *ptr_X, *ptr_X0;
767: const PetscReal hs = (ctx.xmax - ctx.xmin) / 2.0 / count_slow;
768: const PetscReal hf = (ctx.xmax - ctx.xmin) / 2.0 / count_fast;
769: PetscCall(TSSolve(ts, X));
770: PetscCall(TSGetSolveTime(ts, &ptime));
771: PetscCall(TSGetStepNumber(ts, &steps));
772: /* calculate the total mass at initial time and final time */
773: mass_initial = 0.0;
774: mass_final = 0.0;
775: PetscCall(DMDAVecGetArrayRead(da, X0, (void *)&ptr_X0));
776: PetscCall(DMDAVecGetArrayRead(da, X, (void *)&ptr_X));
777: for (i = xs; i < xs + xm; i++) {
778: if (i < ctx.sf || i > ctx.fs - 1) {
779: for (k = 0; k < dof; k++) {
780: mass_initial = mass_initial + hs * ptr_X0[i * dof + k];
781: mass_final = mass_final + hs * ptr_X[i * dof + k];
782: }
783: } else {
784: for (k = 0; k < dof; k++) {
785: mass_initial = mass_initial + hf * ptr_X0[i * dof + k];
786: mass_final = mass_final + hf * ptr_X[i * dof + k];
787: }
788: }
789: }
790: PetscCall(DMDAVecRestoreArrayRead(da, X0, (void *)&ptr_X0));
791: PetscCall(DMDAVecRestoreArrayRead(da, X, (void *)&ptr_X));
792: mass_difference = mass_final - mass_initial;
793: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &mass_difference, 1, MPIU_SCALAR, MPIU_SUM, comm));
794: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
795: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
796: if (ctx.exact) {
797: PetscReal nrm1 = 0;
798: PetscCall(SolutionErrorNorms(&ctx, da, ptime, X, &nrm1));
799: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
800: }
801: if (ctx.simulation) {
802: PetscReal nrm1 = 0;
803: PetscViewer fd;
804: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
805: Vec XR;
806: PetscBool flg;
807: const PetscScalar *ptr_XR;
808: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
809: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
810: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
811: PetscCall(VecDuplicate(X0, &XR));
812: PetscCall(VecLoad(XR, fd));
813: PetscCall(PetscViewerDestroy(&fd));
814: PetscCall(VecGetArrayRead(X, &ptr_X));
815: PetscCall(VecGetArrayRead(XR, &ptr_XR));
816: for (i = 0; i < Mx; i++) {
817: if (i < count_slow / 2 || i > count_slow / 2 + count_fast - 1) nrm1 = nrm1 + hs * PetscAbs(ptr_X[i] - ptr_XR[i]);
818: else nrm1 = nrm1 + hf * PetscAbs(ptr_X[i] - ptr_XR[i]);
819: }
820: PetscCall(VecRestoreArrayRead(X, &ptr_X));
821: PetscCall(VecRestoreArrayRead(XR, &ptr_XR));
822: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g\n", (double)nrm1));
823: PetscCall(VecDestroy(&XR));
824: }
825: }
827: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
828: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
829: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
830: if (draw & 0x4) {
831: Vec Y;
832: PetscCall(VecDuplicate(X, &Y));
833: PetscCall(FVSample(&ctx, da, ptime, Y));
834: PetscCall(VecAYPX(Y, -1, X));
835: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
836: PetscCall(VecDestroy(&Y));
837: }
839: if (view_final) {
840: PetscViewer viewer;
841: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
842: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
843: PetscCall(VecView(X, viewer));
844: PetscCall(PetscViewerPopFormat(viewer));
845: PetscCall(PetscViewerDestroy(&viewer));
846: }
848: /* Clean up */
849: PetscCall((*ctx.physics.destroy)(ctx.physics.user));
850: for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
851: PetscCall(PetscFree3(ctx.u, ctx.flux, ctx.speeds));
852: PetscCall(ISDestroy(&ctx.iss));
853: PetscCall(ISDestroy(&ctx.isf));
854: PetscCall(VecDestroy(&X));
855: PetscCall(VecDestroy(&X0));
856: PetscCall(VecDestroy(&R));
857: PetscCall(DMDestroy(&da));
858: PetscCall(TSDestroy(&ts));
859: PetscCall(PetscFree(index_slow));
860: PetscCall(PetscFree(index_fast));
861: PetscCall(PetscFunctionListDestroy(&physics));
862: PetscCall(PetscFinalize());
863: return 0;
864: }
866: /*TEST
868: build:
869: requires: !complex
871: test:
872: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0
874: test:
875: suffix: 2
876: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
877: output_file: output/ex7_1.out
879: test:
880: suffix: 3
881: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
883: test:
884: suffix: 4
885: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
886: output_file: output/ex7_3.out
888: test:
889: suffix: 5
890: nsize: 2
891: args: -da_grid_x 60 -initial 7 -xmin -1 -xmax 1 -hratio 2 -ts_time_step 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
892: output_file: output/ex7_3.out
893: TEST*/