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