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