Actual source code: ex5.c
1: /*
2: Note:
3: To performance a convergence study, a reference solution can be generated with a small stepsize and stored into a binary file, e.g. -ts_type ssp -ts_dt 1e-5 -ts_max_steps 30000 -ts_view_solution binary:reference.bin
4: Errors can be computed in the following runs with -simulation -f reference.bin
6: Multirate RK options:
7: -ts_rk_dtratio is the ratio between larger time step size and small time step size
8: -ts_rk_multirate_type has three choices: none (for single step RK)
9: combined (for multirate method and user just need to provide the same RHS with the single step RK)
10: partitioned (for multiraet method and user need to provide two RHS, one is for fast components and the orther is for slow components
11: */
13: static const char help[] = "1D periodic Finite Volume solver in slope-limiter form with semidiscrete time stepping.\n"
14: " advection - Variable coefficient scalar advection\n"
15: " u_t + (a*u)_x = 0\n"
16: " a is a piecewise function with two values say a0 and a1 (both a0 and a1 are constant).\n"
17: " in this toy problem we assume a=a0 when x<0 and a=a1 when x>0 with a0<a1 which has the same discontinuous point with initial condition.\n"
18: " we don't provide an exact solution, so you need to generate reference solution to do convergnece staudy,\n"
19: " more precisely, you need firstly generate a reference to a binary file say file.bin, then on command line,\n"
20: " you should type -simulation -f file.bin.\n"
21: " you can choose the number of grids by -da_grid_x.\n"
22: " you can choose the value of a by -physics_advect_a1 and -physics_advect_a2.\n";
24: #include <petscts.h>
25: #include <petscdm.h>
26: #include <petscdmda.h>
27: #include <petscdraw.h>
28: #include <petsc/private/tsimpl.h>
30: #include <petsc/private/kernels/blockinvert.h>
32: #include "finitevolume1d.h"
34: static inline PetscReal RangeMod(PetscReal a, PetscReal xmin, PetscReal xmax)
35: {
36: PetscReal range = xmax - xmin;
37: return xmin + PetscFmodReal(range + PetscFmodReal(a, range), range);
38: }
40: /* --------------------------------- Advection ----------------------------------- */
42: typedef struct {
43: PetscReal a[2]; /* advective velocity */
44: } AdvectCtx;
46: static PetscErrorCode PhysicsRiemann_Advect(void *vctx, PetscInt m, const PetscScalar *uL, const PetscScalar *uR, PetscScalar *flux, PetscReal *maxspeed, PetscReal x, PetscReal xmin, PetscReal xmax)
47: {
48: AdvectCtx *ctx = (AdvectCtx *)vctx;
49: PetscReal *speed;
51: PetscFunctionBeginUser;
52: speed = ctx->a;
53: if (x == 0 || x == xmin || x == xmax)
54: flux[0] = PetscMax(0, (speed[0] + speed[1]) / 2.0) * uL[0] + PetscMin(0, (speed[0] + speed[1]) / 2.0) * uR[0]; /* if condition need to be changed base on different problem, '0' is the discontinuous point of a */
55: else if (x < 0) flux[0] = PetscMax(0, speed[0]) * uL[0] + PetscMin(0, speed[0]) * uR[0]; /* else if condition need to be changed based on different problem, 'x = 0' is discontinuous point of a */
56: else flux[0] = PetscMax(0, speed[1]) * uL[0] + PetscMin(0, speed[1]) * uR[0];
57: *maxspeed = *speed;
58: PetscFunctionReturn(PETSC_SUCCESS);
59: }
61: static PetscErrorCode PhysicsCharacteristic_Advect(void *vctx, PetscInt m, const PetscScalar *u, PetscScalar *X, PetscScalar *Xi, PetscReal *speeds, PetscReal x)
62: {
63: AdvectCtx *ctx = (AdvectCtx *)vctx;
65: PetscFunctionBeginUser;
66: X[0] = 1.;
67: Xi[0] = 1.;
68: if (x < 0) speeds[0] = ctx->a[0]; /* x is discontinuous point of a */
69: else speeds[0] = ctx->a[1];
70: PetscFunctionReturn(PETSC_SUCCESS);
71: }
73: static PetscErrorCode PhysicsSample_Advect(void *vctx, PetscInt initial, FVBCType bctype, PetscReal xmin, PetscReal xmax, PetscReal t, PetscReal x, PetscReal *u)
74: {
75: AdvectCtx *ctx = (AdvectCtx *)vctx;
76: PetscReal *a = ctx->a, x0;
78: PetscFunctionBeginUser;
79: if (x < 0) { /* x is cell center */
80: switch (bctype) {
81: case FVBC_OUTFLOW:
82: x0 = x - a[0] * t;
83: break;
84: case FVBC_PERIODIC:
85: x0 = RangeMod(x - a[0] * t, xmin, xmax);
86: break;
87: default:
88: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
89: }
90: switch (initial) {
91: case 0:
92: u[0] = (x0 < 0) ? 1 : -1;
93: break;
94: case 1:
95: u[0] = (x0 < 0) ? -1 : 1;
96: break;
97: case 2:
98: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
99: break;
100: case 3:
101: u[0] = PetscSinReal(2 * PETSC_PI * x0);
102: break;
103: case 4:
104: u[0] = PetscAbs(x0);
105: break;
106: case 5:
107: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
108: break;
109: case 6:
110: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
111: break;
112: case 7:
113: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
114: break;
115: default:
116: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
117: }
118: } else {
119: switch (bctype) {
120: case FVBC_OUTFLOW:
121: x0 = x - a[1] * t;
122: break;
123: case FVBC_PERIODIC:
124: x0 = RangeMod(x - a[1] * t, xmin, xmax);
125: break;
126: default:
127: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown BCType");
128: }
129: switch (initial) {
130: case 0:
131: u[0] = (x0 < 0) ? 1 : -1;
132: break;
133: case 1:
134: u[0] = (x0 < 0) ? -1 : 1;
135: break;
136: case 2:
137: u[0] = (0 < x0 && x0 < 1) ? 1 : 0;
138: break;
139: case 3:
140: u[0] = PetscSinReal(2 * PETSC_PI * x0);
141: break;
142: case 4:
143: u[0] = PetscAbs(x0);
144: break;
145: case 5:
146: u[0] = (x0 < 0 || x0 > 0.5) ? 0 : PetscSqr(PetscSinReal(2 * PETSC_PI * x0));
147: break;
148: case 6:
149: u[0] = (x0 < 0) ? 0 : ((x0 < 1) ? x0 : ((x0 < 2) ? 2 - x0 : 0));
150: break;
151: case 7:
152: u[0] = PetscPowReal(PetscSinReal(PETSC_PI * x0), 10.0);
153: break;
154: default:
155: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "unknown initial condition");
156: }
157: }
158: PetscFunctionReturn(PETSC_SUCCESS);
159: }
161: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
162: {
163: AdvectCtx *user;
165: PetscFunctionBeginUser;
166: PetscCall(PetscNew(&user));
167: ctx->physics.sample = PhysicsSample_Advect;
168: ctx->physics.riemann = PhysicsRiemann_Advect;
169: ctx->physics.characteristic = PhysicsCharacteristic_Advect;
170: ctx->physics.destroy = PhysicsDestroy_SimpleFree;
171: ctx->physics.user = user;
172: ctx->physics.dof = 1;
173: PetscCall(PetscStrallocpy("u", &ctx->physics.fieldname[0]));
174: user->a[0] = 1;
175: user->a[1] = 1;
176: PetscOptionsBegin(ctx->comm, ctx->prefix, "Options for advection", "");
177: {
178: PetscCall(PetscOptionsReal("-physics_advect_a1", "Speed1", "", user->a[0], &user->a[0], NULL));
179: PetscCall(PetscOptionsReal("-physics_advect_a2", "Speed2", "", user->a[1], &user->a[1], NULL));
180: }
181: PetscOptionsEnd();
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: /* --------------------------------- Finite Volume Solver for slow parts ----------------------------------- */
187: PetscErrorCode FVRHSFunctionslow(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
188: {
189: FVCtx *ctx = (FVCtx *)vctx;
190: PetscInt i, j, k, Mx, dof, xs, xm, len_slow;
191: PetscReal hx, cfl_idt = 0;
192: PetscScalar *x, *f, *slope;
193: Vec Xloc;
194: DM da;
196: PetscFunctionBeginUser;
197: PetscCall(TSGetDM(ts, &da));
198: PetscCall(DMGetLocalVector(da, &Xloc));
199: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
200: hx = (ctx->xmax - ctx->xmin) / Mx;
201: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
202: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
203: PetscCall(VecZeroEntries(F));
204: PetscCall(DMDAVecGetArray(da, Xloc, &x));
205: PetscCall(VecGetArray(F, &f));
206: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
207: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
208: PetscCall(ISGetSize(ctx->iss, &len_slow));
210: if (ctx->bctype == FVBC_OUTFLOW) {
211: for (i = xs - 2; i < 0; i++) {
212: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
213: }
214: for (i = Mx; i < xs + xm + 2; i++) {
215: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
216: }
217: }
218: for (i = xs - 1; i < xs + xm + 1; i++) {
219: struct _LimitInfo info;
220: PetscScalar *cjmpL, *cjmpR;
221: if (i < len_slow + 1) {
222: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
223: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
224: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
225: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
226: cjmpL = &ctx->cjmpLR[0];
227: cjmpR = &ctx->cjmpLR[dof];
228: for (j = 0; j < dof; j++) {
229: PetscScalar jmpL, jmpR;
230: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
231: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
232: for (k = 0; k < dof; k++) {
233: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
234: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
235: }
236: }
237: /* Apply limiter to the left and right characteristic jumps */
238: info.m = dof;
239: info.hx = hx;
240: (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
241: for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
242: for (j = 0; j < dof; j++) {
243: PetscScalar tmp = 0;
244: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
245: slope[i * dof + j] = tmp;
246: }
247: }
248: }
250: for (i = xs; i < xs + xm + 1; i++) {
251: PetscReal maxspeed;
252: PetscScalar *uL, *uR;
253: if (i < len_slow) { /* slow parts can be changed based on a */
254: uL = &ctx->uLR[0];
255: uR = &ctx->uLR[dof];
256: for (j = 0; j < dof; j++) {
257: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
258: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
259: }
260: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
261: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
262: if (i > xs) {
263: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
264: }
265: if (i < xs + xm) {
266: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx;
267: }
268: }
269: if (i == len_slow) { /* slow parts can be changed based on a */
270: uL = &ctx->uLR[0];
271: uR = &ctx->uLR[dof];
272: for (j = 0; j < dof; j++) {
273: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
274: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
275: }
276: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
277: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
278: if (i > xs) {
279: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
280: }
281: }
282: }
283: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
284: PetscCall(VecRestoreArray(F, &f));
285: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
286: PetscCall(DMRestoreLocalVector(da, &Xloc));
287: PetscFunctionReturn(PETSC_SUCCESS);
288: }
290: /* --------------------------------- Finite Volume Solver for fast parts ----------------------------------- */
292: PetscErrorCode FVRHSFunctionfast(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
293: {
294: FVCtx *ctx = (FVCtx *)vctx;
295: PetscInt i, j, k, Mx, dof, xs, xm, len_slow;
296: PetscReal hx, cfl_idt = 0;
297: PetscScalar *x, *f, *slope;
298: Vec Xloc;
299: DM da;
301: PetscFunctionBeginUser;
302: PetscCall(TSGetDM(ts, &da));
303: PetscCall(DMGetLocalVector(da, &Xloc));
304: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
305: hx = (ctx->xmax - ctx->xmin) / Mx;
306: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
307: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
308: PetscCall(VecZeroEntries(F));
309: PetscCall(DMDAVecGetArray(da, Xloc, &x));
310: PetscCall(VecGetArray(F, &f));
311: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
312: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
313: PetscCall(ISGetSize(ctx->iss, &len_slow));
315: if (ctx->bctype == FVBC_OUTFLOW) {
316: for (i = xs - 2; i < 0; i++) {
317: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
318: }
319: for (i = Mx; i < xs + xm + 2; i++) {
320: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
321: }
322: }
323: for (i = xs - 1; i < xs + xm + 1; i++) {
324: struct _LimitInfo info;
325: PetscScalar *cjmpL, *cjmpR;
326: if (i > len_slow - 2) {
327: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
328: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
329: cjmpL = &ctx->cjmpLR[0];
330: cjmpR = &ctx->cjmpLR[dof];
331: for (j = 0; j < dof; j++) {
332: PetscScalar jmpL, jmpR;
333: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
334: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
335: for (k = 0; k < dof; k++) {
336: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
337: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
338: }
339: }
340: /* Apply limiter to the left and right characteristic jumps */
341: info.m = dof;
342: info.hx = hx;
343: (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
344: for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
345: for (j = 0; j < dof; j++) {
346: PetscScalar tmp = 0;
347: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
348: slope[i * dof + j] = tmp;
349: }
350: }
351: }
353: for (i = xs; i < xs + xm + 1; i++) {
354: PetscReal maxspeed;
355: PetscScalar *uL, *uR;
356: if (i > len_slow) {
357: uL = &ctx->uLR[0];
358: uR = &ctx->uLR[dof];
359: for (j = 0; j < dof; j++) {
360: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
361: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
362: }
363: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
364: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
365: if (i > xs) {
366: for (j = 0; j < dof; j++) f[(i - len_slow - 1) * dof + j] -= ctx->flux[j] / hx;
367: }
368: if (i < xs + xm) {
369: for (j = 0; j < dof; j++) f[(i - len_slow) * dof + j] += ctx->flux[j] / hx;
370: }
371: }
372: if (i == len_slow) {
373: uL = &ctx->uLR[0];
374: uR = &ctx->uLR[dof];
375: for (j = 0; j < dof; j++) {
376: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
377: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
378: }
379: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
380: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
381: if (i < xs + xm) {
382: for (j = 0; j < dof; j++) f[(i - len_slow) * dof + j] += ctx->flux[j] / hx;
383: }
384: }
385: }
386: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
387: PetscCall(VecRestoreArray(F, &f));
388: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
389: PetscCall(DMRestoreLocalVector(da, &Xloc));
390: PetscFunctionReturn(PETSC_SUCCESS);
391: }
393: PetscErrorCode FVRHSFunctionslow2(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
394: {
395: FVCtx *ctx = (FVCtx *)vctx;
396: PetscInt i, j, k, Mx, dof, xs, xm, len_slow1, len_slow2;
397: PetscReal hx, cfl_idt = 0;
398: PetscScalar *x, *f, *slope;
399: Vec Xloc;
400: DM da;
402: PetscFunctionBeginUser;
403: PetscCall(TSGetDM(ts, &da));
404: PetscCall(DMGetLocalVector(da, &Xloc));
405: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
406: hx = (ctx->xmax - ctx->xmin) / Mx;
407: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
408: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
409: PetscCall(VecZeroEntries(F));
410: PetscCall(DMDAVecGetArray(da, Xloc, &x));
411: PetscCall(VecGetArray(F, &f));
412: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
413: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
414: PetscCall(ISGetSize(ctx->iss, &len_slow1));
415: PetscCall(ISGetSize(ctx->iss2, &len_slow2));
416: if (ctx->bctype == FVBC_OUTFLOW) {
417: for (i = xs - 2; i < 0; i++) {
418: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
419: }
420: for (i = Mx; i < xs + xm + 2; i++) {
421: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
422: }
423: }
424: for (i = xs - 1; i < xs + xm + 1; i++) {
425: struct _LimitInfo info;
426: PetscScalar *cjmpL, *cjmpR;
427: if (i < len_slow1 + len_slow2 + 1 && i > len_slow1 - 2) {
428: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
429: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
430: cjmpL = &ctx->cjmpLR[0];
431: cjmpR = &ctx->cjmpLR[dof];
432: for (j = 0; j < dof; j++) {
433: PetscScalar jmpL, jmpR;
434: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
435: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
436: for (k = 0; k < dof; k++) {
437: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
438: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
439: }
440: }
441: /* Apply limiter to the left and right characteristic jumps */
442: info.m = dof;
443: info.hx = hx;
444: (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
445: for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
446: for (j = 0; j < dof; j++) {
447: PetscScalar tmp = 0;
448: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
449: slope[i * dof + j] = tmp;
450: }
451: }
452: }
454: for (i = xs; i < xs + xm + 1; i++) {
455: PetscReal maxspeed;
456: PetscScalar *uL, *uR;
457: if (i < len_slow1 + len_slow2 && i > len_slow1) { /* slow parts can be changed based on a */
458: uL = &ctx->uLR[0];
459: uR = &ctx->uLR[dof];
460: for (j = 0; j < dof; j++) {
461: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
462: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
463: }
464: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
465: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
466: if (i > xs) {
467: for (j = 0; j < dof; j++) f[(i - len_slow1 - 1) * dof + j] -= ctx->flux[j] / hx;
468: }
469: if (i < xs + xm) {
470: for (j = 0; j < dof; j++) f[(i - len_slow1) * dof + j] += ctx->flux[j] / hx;
471: }
472: }
473: if (i == len_slow1 + len_slow2) { /* slow parts can be changed based on a */
474: uL = &ctx->uLR[0];
475: uR = &ctx->uLR[dof];
476: for (j = 0; j < dof; j++) {
477: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
478: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
479: }
480: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
481: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
482: if (i > xs) {
483: for (j = 0; j < dof; j++) f[(i - len_slow1 - 1) * dof + j] -= ctx->flux[j] / hx;
484: }
485: }
486: if (i == len_slow1) { /* slow parts can be changed based on a */
487: uL = &ctx->uLR[0];
488: uR = &ctx->uLR[dof];
489: for (j = 0; j < dof; j++) {
490: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
491: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
492: }
493: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
494: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
495: if (i < xs + xm) {
496: for (j = 0; j < dof; j++) f[(i - len_slow1) * dof + j] += ctx->flux[j] / hx;
497: }
498: }
499: }
501: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
502: PetscCall(VecRestoreArray(F, &f));
503: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
504: PetscCall(DMRestoreLocalVector(da, &Xloc));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: PetscErrorCode FVRHSFunctionfast2(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
509: {
510: FVCtx *ctx = (FVCtx *)vctx;
511: PetscInt i, j, k, Mx, dof, xs, xm, len_slow1, len_slow2;
512: PetscReal hx, cfl_idt = 0;
513: PetscScalar *x, *f, *slope;
514: Vec Xloc;
515: DM da;
517: PetscFunctionBeginUser;
518: PetscCall(TSGetDM(ts, &da));
519: PetscCall(DMGetLocalVector(da, &Xloc));
520: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
521: hx = (ctx->xmax - ctx->xmin) / Mx;
522: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
523: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
524: PetscCall(VecZeroEntries(F));
525: PetscCall(DMDAVecGetArray(da, Xloc, &x));
526: PetscCall(VecGetArray(F, &f));
527: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope));
528: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
529: PetscCall(ISGetSize(ctx->iss, &len_slow1));
530: PetscCall(ISGetSize(ctx->iss2, &len_slow2));
532: if (ctx->bctype == FVBC_OUTFLOW) {
533: for (i = xs - 2; i < 0; i++) {
534: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
535: }
536: for (i = Mx; i < xs + xm + 2; i++) {
537: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
538: }
539: }
540: for (i = xs - 1; i < xs + xm + 1; i++) {
541: struct _LimitInfo info;
542: PetscScalar *cjmpL, *cjmpR;
543: if (i > len_slow1 + len_slow2 - 2) {
544: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
545: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
546: cjmpL = &ctx->cjmpLR[0];
547: cjmpR = &ctx->cjmpLR[dof];
548: for (j = 0; j < dof; j++) {
549: PetscScalar jmpL, jmpR;
550: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
551: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
552: for (k = 0; k < dof; k++) {
553: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
554: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
555: }
556: }
557: /* Apply limiter to the left and right characteristic jumps */
558: info.m = dof;
559: info.hx = hx;
560: (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
561: for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
562: for (j = 0; j < dof; j++) {
563: PetscScalar tmp = 0;
564: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
565: slope[i * dof + j] = tmp;
566: }
567: }
568: }
570: for (i = xs; i < xs + xm + 1; i++) {
571: PetscReal maxspeed;
572: PetscScalar *uL, *uR;
573: if (i > len_slow1 + len_slow2) {
574: uL = &ctx->uLR[0];
575: uR = &ctx->uLR[dof];
576: for (j = 0; j < dof; j++) {
577: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
578: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
579: }
580: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
581: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
582: if (i > xs) {
583: for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2 - 1) * dof + j] -= ctx->flux[j] / hx;
584: }
585: if (i < xs + xm) {
586: for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2) * dof + j] += ctx->flux[j] / hx;
587: }
588: }
589: if (i == len_slow1 + len_slow2) {
590: uL = &ctx->uLR[0];
591: uR = &ctx->uLR[dof];
592: for (j = 0; j < dof; j++) {
593: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
594: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
595: }
596: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
597: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
598: if (i < xs + xm) {
599: for (j = 0; j < dof; j++) f[(i - len_slow1 - len_slow2) * dof + j] += ctx->flux[j] / hx;
600: }
601: }
602: }
603: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
604: PetscCall(VecRestoreArray(F, &f));
605: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
606: PetscCall(DMRestoreLocalVector(da, &Xloc));
607: PetscFunctionReturn(PETSC_SUCCESS);
608: }
610: int main(int argc, char *argv[])
611: {
612: char lname[256] = "mc", physname[256] = "advect", final_fname[256] = "solution.m";
613: PetscFunctionList limiters = 0, physics = 0;
614: MPI_Comm comm;
615: TS ts;
616: DM da;
617: Vec X, X0, R;
618: FVCtx ctx;
619: PetscInt i, k, dof, xs, xm, Mx, draw = 0, *index_slow, *index_fast, islow = 0, ifast = 0;
620: PetscBool view_final = PETSC_FALSE;
621: PetscReal ptime;
623: PetscFunctionBeginUser;
624: PetscCall(PetscInitialize(&argc, &argv, 0, help));
625: comm = PETSC_COMM_WORLD;
626: PetscCall(PetscMemzero(&ctx, sizeof(ctx)));
628: /* Register limiters to be available on the command line */
629: PetscCall(PetscFunctionListAdd(&limiters, "upwind", Limit_Upwind));
630: PetscCall(PetscFunctionListAdd(&limiters, "lax-wendroff", Limit_LaxWendroff));
631: PetscCall(PetscFunctionListAdd(&limiters, "beam-warming", Limit_BeamWarming));
632: PetscCall(PetscFunctionListAdd(&limiters, "fromm", Limit_Fromm));
633: PetscCall(PetscFunctionListAdd(&limiters, "minmod", Limit_Minmod));
634: PetscCall(PetscFunctionListAdd(&limiters, "superbee", Limit_Superbee));
635: PetscCall(PetscFunctionListAdd(&limiters, "mc", Limit_MC));
636: PetscCall(PetscFunctionListAdd(&limiters, "vanleer", Limit_VanLeer));
637: PetscCall(PetscFunctionListAdd(&limiters, "vanalbada", Limit_VanAlbada));
638: PetscCall(PetscFunctionListAdd(&limiters, "vanalbadatvd", Limit_VanAlbadaTVD));
639: PetscCall(PetscFunctionListAdd(&limiters, "koren", Limit_Koren));
640: PetscCall(PetscFunctionListAdd(&limiters, "korensym", Limit_KorenSym));
641: PetscCall(PetscFunctionListAdd(&limiters, "koren3", Limit_Koren3));
642: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon2", Limit_CadaTorrilhon2));
643: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r0p1", Limit_CadaTorrilhon3R0p1));
644: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r1", Limit_CadaTorrilhon3R1));
645: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r10", Limit_CadaTorrilhon3R10));
646: PetscCall(PetscFunctionListAdd(&limiters, "cada-torrilhon3-r100", Limit_CadaTorrilhon3R100));
648: /* Register physical models to be available on the command line */
649: PetscCall(PetscFunctionListAdd(&physics, "advect", PhysicsCreate_Advect));
651: ctx.comm = comm;
652: ctx.cfl = 0.9;
653: ctx.bctype = FVBC_PERIODIC;
654: ctx.xmin = -1.0;
655: ctx.xmax = 1.0;
656: PetscOptionsBegin(comm, NULL, "Finite Volume solver options", "");
657: PetscCall(PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL));
658: PetscCall(PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL));
659: PetscCall(PetscOptionsFList("-limit", "Name of flux limiter to use", "", limiters, lname, lname, sizeof(lname), NULL));
660: PetscCall(PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL));
661: PetscCall(PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final));
662: PetscCall(PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL));
663: PetscCall(PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL));
664: PetscCall(PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL));
665: PetscCall(PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL));
666: PetscCall(PetscOptionsBool("-recursive_split", "Split the domain recursively", "", ctx.recursive, &ctx.recursive, NULL));
667: PetscOptionsEnd();
669: /* Choose the limiter from the list of registered limiters */
670: PetscCall(PetscFunctionListFind(limiters, lname, &ctx.limit));
671: PetscCheck(ctx.limit, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Limiter '%s' not found", lname);
673: /* Choose the physics from the list of registered models */
674: {
675: PetscErrorCode (*r)(FVCtx *);
676: PetscCall(PetscFunctionListFind(physics, physname, &r));
677: PetscCheck(r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Physics '%s' not found", physname);
678: /* Create the physics, will set the number of fields and their names */
679: PetscCall((*r)(&ctx));
680: }
682: /* Create a DMDA to manage the parallel grid */
683: PetscCall(DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da));
684: PetscCall(DMSetFromOptions(da));
685: PetscCall(DMSetUp(da));
686: /* Inform the DMDA of the field names provided by the physics. */
687: /* The names will be shown in the title bars when run with -ts_monitor_draw_solution */
688: for (i = 0; i < ctx.physics.dof; i++) PetscCall(DMDASetFieldName(da, i, ctx.physics.fieldname[i]));
689: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
690: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
692: /* Set coordinates of cell centers */
693: 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));
695: /* Allocate work space for the Finite Volume solver (so it doesn't have to be reallocated on each function evaluation) */
696: PetscCall(PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope));
697: PetscCall(PetscMalloc3(2 * dof, &ctx.uLR, dof, &ctx.flux, dof, &ctx.speeds));
699: /* Create a vector to store the solution and to save the initial state */
700: PetscCall(DMCreateGlobalVector(da, &X));
701: PetscCall(VecDuplicate(X, &X0));
702: PetscCall(VecDuplicate(X, &R));
704: /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
705: PetscCall(PetscMalloc1(xm * dof, &index_slow));
706: PetscCall(PetscMalloc1(xm * dof, &index_fast));
707: for (i = xs; i < xs + xm; i++) {
708: if (ctx.xmin + i * (ctx.xmax - ctx.xmin) / (PetscReal)Mx + 0.5 * (ctx.xmax - ctx.xmin) / (PetscReal)Mx < 0)
709: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
710: else
711: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
712: } /* this step need to be changed based on discontinuous point of a */
713: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss));
714: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf));
716: /* Create a time-stepping object */
717: PetscCall(TSCreate(comm, &ts));
718: PetscCall(TSSetDM(ts, da));
719: PetscCall(TSSetRHSFunction(ts, R, FVRHSFunction, &ctx));
720: PetscCall(TSRHSSplitSetIS(ts, "slow", ctx.iss));
721: PetscCall(TSRHSSplitSetIS(ts, "fast", ctx.isf));
722: PetscCall(TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx));
723: PetscCall(TSRHSSplitSetRHSFunction(ts, "fast", NULL, FVRHSFunctionfast, &ctx));
725: if (ctx.recursive) {
726: TS subts;
727: islow = 0;
728: ifast = 0;
729: for (i = xs; i < xs + xm; i++) {
730: PetscReal coord = ctx.xmin + i * (ctx.xmax - ctx.xmin) / (PetscReal)Mx + 0.5 * (ctx.xmax - ctx.xmin) / (PetscReal)Mx;
731: if (coord >= 0 && coord < ctx.xmin + (ctx.xmax - ctx.xmin) * 3 / 4.)
732: for (k = 0; k < dof; k++) index_slow[islow++] = i * dof + k;
733: if (coord >= ctx.xmin + (ctx.xmax - ctx.xmin) * 3 / 4.)
734: for (k = 0; k < dof; k++) index_fast[ifast++] = i * dof + k;
735: } /* this step need to be changed based on discontinuous point of a */
736: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss2));
737: PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf2));
739: PetscCall(TSRHSSplitGetSubTS(ts, "fast", &subts));
740: PetscCall(TSRHSSplitSetIS(subts, "slow", ctx.iss2));
741: PetscCall(TSRHSSplitSetIS(subts, "fast", ctx.isf2));
742: PetscCall(TSRHSSplitSetRHSFunction(subts, "slow", NULL, FVRHSFunctionslow2, &ctx));
743: PetscCall(TSRHSSplitSetRHSFunction(subts, "fast", NULL, FVRHSFunctionfast2, &ctx));
744: }
746: /*PetscCall(TSSetType(ts,TSSSP));*/
747: PetscCall(TSSetType(ts, TSMPRK));
748: PetscCall(TSSetMaxTime(ts, 10));
749: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
751: /* Compute initial conditions and starting time step */
752: PetscCall(FVSample(&ctx, da, 0, X0));
753: PetscCall(FVRHSFunction(ts, 0, X0, X, (void *)&ctx)); /* Initial function evaluation, only used to determine max speed */
754: PetscCall(VecCopy(X0, X)); /* The function value was not used so we set X=X0 again */
755: PetscCall(TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt));
756: PetscCall(TSSetFromOptions(ts)); /* Take runtime options */
757: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
758: {
759: PetscInt steps;
760: PetscScalar mass_initial, mass_final, mass_difference;
762: PetscCall(TSSolve(ts, X));
763: PetscCall(TSGetSolveTime(ts, &ptime));
764: PetscCall(TSGetStepNumber(ts, &steps));
765: PetscCall(PetscPrintf(comm, "Final time %g, steps %" PetscInt_FMT "\n", (double)ptime, steps));
766: /* calculate the total mass at initial time and final time */
767: mass_initial = 0.0;
768: mass_final = 0.0;
769: PetscCall(VecSum(X0, &mass_initial));
770: PetscCall(VecSum(X, &mass_final));
771: mass_difference = (ctx.xmax - ctx.xmin) / (PetscScalar)Mx * (mass_final - mass_initial);
772: PetscCall(PetscPrintf(comm, "Mass difference %g\n", (double)mass_difference));
773: if (ctx.simulation) {
774: PetscViewer fd;
775: char filename[PETSC_MAX_PATH_LEN] = "binaryoutput";
776: Vec XR;
777: PetscReal nrm1, nrmsup;
778: PetscBool flg;
780: PetscCall(PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg));
781: PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_USER, "Must indicate binary file with the -f option");
782: PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd));
783: PetscCall(VecDuplicate(X0, &XR));
784: PetscCall(VecLoad(XR, fd));
785: PetscCall(PetscViewerDestroy(&fd));
786: PetscCall(VecAYPX(XR, -1, X));
787: PetscCall(VecNorm(XR, NORM_1, &nrm1));
788: PetscCall(VecNorm(XR, NORM_INFINITY, &nrmsup));
789: nrm1 /= Mx;
790: PetscCall(PetscPrintf(comm, "Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n", (double)nrm1, (double)nrmsup));
791: PetscCall(VecDestroy(&XR));
792: }
793: }
795: PetscCall(SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD));
796: if (draw & 0x1) PetscCall(VecView(X0, PETSC_VIEWER_DRAW_WORLD));
797: if (draw & 0x2) PetscCall(VecView(X, PETSC_VIEWER_DRAW_WORLD));
798: if (draw & 0x4) {
799: Vec Y;
800: PetscCall(VecDuplicate(X, &Y));
801: PetscCall(FVSample(&ctx, da, ptime, Y));
802: PetscCall(VecAYPX(Y, -1, X));
803: PetscCall(VecView(Y, PETSC_VIEWER_DRAW_WORLD));
804: PetscCall(VecDestroy(&Y));
805: }
807: if (view_final) {
808: PetscViewer viewer;
809: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer));
810: PetscCall(PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB));
811: PetscCall(VecView(X, viewer));
812: PetscCall(PetscViewerPopFormat(viewer));
813: PetscCall(PetscViewerDestroy(&viewer));
814: }
816: /* Clean up */
817: PetscCall((*ctx.physics.destroy)(ctx.physics.user));
818: for (i = 0; i < ctx.physics.dof; i++) PetscCall(PetscFree(ctx.physics.fieldname[i]));
819: PetscCall(PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope));
820: PetscCall(PetscFree3(ctx.uLR, ctx.flux, ctx.speeds));
821: PetscCall(ISDestroy(&ctx.iss));
822: PetscCall(ISDestroy(&ctx.isf));
823: PetscCall(ISDestroy(&ctx.iss2));
824: PetscCall(ISDestroy(&ctx.isf2));
825: PetscCall(VecDestroy(&X));
826: PetscCall(VecDestroy(&X0));
827: PetscCall(VecDestroy(&R));
828: PetscCall(DMDestroy(&da));
829: PetscCall(TSDestroy(&ts));
830: PetscCall(PetscFree(index_slow));
831: PetscCall(PetscFree(index_fast));
832: PetscCall(PetscFunctionListDestroy(&limiters));
833: PetscCall(PetscFunctionListDestroy(&physics));
834: PetscCall(PetscFinalize());
835: return 0;
836: }
838: /*TEST
839: build:
840: requires: !complex
841: depends: finitevolume1d.c
843: test:
844: args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0
846: test:
847: suffix: 2
848: args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1
849: output_file: output/ex5_1.out
851: test:
852: suffix: 3
853: args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 0
855: test:
856: suffix: 4
857: args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type mprk -ts_mprk_type 2a22 -ts_use_splitrhsfunction 1
858: output_file: output/ex5_3.out
860: test:
861: suffix: 5
862: args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 0 -recursive_split
864: test:
865: suffix: 6
866: args: -da_grid_x 60 -initial 1 -xmin -1 -xmax 1 -limit mc -physics_advect_a1 1 -physics_advect_a2 2 -ts_dt 0.025 -ts_max_steps 24 -ts_type rk -ts_rk_type 2a -ts_rk_dtratio 2 -ts_rk_multirate -ts_use_splitrhsfunction 1 -recursive_split
867: TEST*/