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 commend 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;
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: return 0;
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;
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: return 0;
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;
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: return 0;
159: }
161: static PetscErrorCode PhysicsCreate_Advect(FVCtx *ctx)
162: {
163: AdvectCtx *user;
166: 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: 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: PetscOptionsReal("-physics_advect_a1", "Speed1", "", user->a[0], &user->a[0], NULL);
179: PetscOptionsReal("-physics_advect_a2", "Speed2", "", user->a[1], &user->a[1], NULL);
180: }
181: PetscOptionsEnd();
182: return 0;
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;
197: TSGetDM(ts, &da);
198: DMGetLocalVector(da, &Xloc);
199: DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
200: hx = (ctx->xmax - ctx->xmin) / Mx;
201: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
202: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
203: VecZeroEntries(F);
204: DMDAVecGetArray(da, Xloc, &x);
205: VecGetArray(F, &f);
206: DMDAGetArray(da, PETSC_TRUE, &slope);
207: DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
208: 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: (*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: 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: (*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: (*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: DMDAVecRestoreArray(da, Xloc, &x);
284: VecRestoreArray(F, &f);
285: DMDARestoreArray(da, PETSC_TRUE, &slope);
286: DMRestoreLocalVector(da, &Xloc);
287: return 0;
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;
302: TSGetDM(ts, &da);
303: DMGetLocalVector(da, &Xloc);
304: DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
305: hx = (ctx->xmax - ctx->xmin) / Mx;
306: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
307: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
308: VecZeroEntries(F);
309: DMDAVecGetArray(da, Xloc, &x);
310: VecGetArray(F, &f);
311: DMDAGetArray(da, PETSC_TRUE, &slope);
312: DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
313: 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: (*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i);
328: 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: (*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: (*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: DMDAVecRestoreArray(da, Xloc, &x);
387: VecRestoreArray(F, &f);
388: DMDARestoreArray(da, PETSC_TRUE, &slope);
389: DMRestoreLocalVector(da, &Xloc);
390: return 0;
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;
403: TSGetDM(ts, &da);
404: DMGetLocalVector(da, &Xloc);
405: DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
406: hx = (ctx->xmax - ctx->xmin) / Mx;
407: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
408: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
409: VecZeroEntries(F);
410: DMDAVecGetArray(da, Xloc, &x);
411: VecGetArray(F, &f);
412: DMDAGetArray(da, PETSC_TRUE, &slope);
413: DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
414: ISGetSize(ctx->iss, &len_slow1);
415: 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: (*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i);
429: 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: (*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: (*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: (*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: DMDAVecRestoreArray(da, Xloc, &x);
502: VecRestoreArray(F, &f);
503: DMDARestoreArray(da, PETSC_TRUE, &slope);
504: DMRestoreLocalVector(da, &Xloc);
505: return 0;
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;
518: TSGetDM(ts, &da);
519: DMGetLocalVector(da, &Xloc);
520: DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
521: hx = (ctx->xmax - ctx->xmin) / Mx;
522: DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc);
523: DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc);
524: VecZeroEntries(F);
525: DMDAVecGetArray(da, Xloc, &x);
526: VecGetArray(F, &f);
527: DMDAGetArray(da, PETSC_TRUE, &slope);
528: DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
529: ISGetSize(ctx->iss, &len_slow1);
530: 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: (*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i);
545: 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: (*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: (*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: DMDAVecRestoreArray(da, Xloc, &x);
604: VecRestoreArray(F, &f);
605: DMDARestoreArray(da, PETSC_TRUE, &slope);
606: DMRestoreLocalVector(da, &Xloc);
607: return 0;
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;
624: PetscInitialize(&argc, &argv, 0, help);
625: comm = PETSC_COMM_WORLD;
626: PetscMemzero(&ctx, sizeof(ctx));
628: /* Register limiters to be available on the command line */
629: PetscFunctionListAdd(&limiters, "upwind", Limit_Upwind);
630: PetscFunctionListAdd(&limiters, "lax-wendroff", Limit_LaxWendroff);
631: PetscFunctionListAdd(&limiters, "beam-warming", Limit_BeamWarming);
632: PetscFunctionListAdd(&limiters, "fromm", Limit_Fromm);
633: PetscFunctionListAdd(&limiters, "minmod", Limit_Minmod);
634: PetscFunctionListAdd(&limiters, "superbee", Limit_Superbee);
635: PetscFunctionListAdd(&limiters, "mc", Limit_MC);
636: PetscFunctionListAdd(&limiters, "vanleer", Limit_VanLeer);
637: PetscFunctionListAdd(&limiters, "vanalbada", Limit_VanAlbada);
638: PetscFunctionListAdd(&limiters, "vanalbadatvd", Limit_VanAlbadaTVD);
639: PetscFunctionListAdd(&limiters, "koren", Limit_Koren);
640: PetscFunctionListAdd(&limiters, "korensym", Limit_KorenSym);
641: PetscFunctionListAdd(&limiters, "koren3", Limit_Koren3);
642: PetscFunctionListAdd(&limiters, "cada-torrilhon2", Limit_CadaTorrilhon2);
643: PetscFunctionListAdd(&limiters, "cada-torrilhon3-r0p1", Limit_CadaTorrilhon3R0p1);
644: PetscFunctionListAdd(&limiters, "cada-torrilhon3-r1", Limit_CadaTorrilhon3R1);
645: PetscFunctionListAdd(&limiters, "cada-torrilhon3-r10", Limit_CadaTorrilhon3R10);
646: PetscFunctionListAdd(&limiters, "cada-torrilhon3-r100", Limit_CadaTorrilhon3R100);
648: /* Register physical models to be available on the command line */
649: 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: PetscOptionsReal("-xmin", "X min", "", ctx.xmin, &ctx.xmin, NULL);
658: PetscOptionsReal("-xmax", "X max", "", ctx.xmax, &ctx.xmax, NULL);
659: PetscOptionsFList("-limit", "Name of flux limiter to use", "", limiters, lname, lname, sizeof(lname), NULL);
660: PetscOptionsInt("-draw", "Draw solution vector, bitwise OR of (1=initial,2=final,4=final error)", "", draw, &draw, NULL);
661: PetscOptionsString("-view_final", "Write final solution in ASCII MATLAB format to given file name", "", final_fname, final_fname, sizeof(final_fname), &view_final);
662: PetscOptionsInt("-initial", "Initial condition (depends on the physics)", "", ctx.initial, &ctx.initial, NULL);
663: PetscOptionsBool("-simulation", "Compare errors with reference solution", "", ctx.simulation, &ctx.simulation, NULL);
664: PetscOptionsReal("-cfl", "CFL number to time step at", "", ctx.cfl, &ctx.cfl, NULL);
665: PetscOptionsEnum("-bc_type", "Boundary condition", "", FVBCTypes, (PetscEnum)ctx.bctype, (PetscEnum *)&ctx.bctype, NULL);
666: 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: PetscFunctionListFind(limiters, lname, &ctx.limit);
673: /* Choose the physics from the list of registered models */
674: {
675: PetscErrorCode (*r)(FVCtx *);
676: PetscFunctionListFind(physics, physname, &r);
678: /* Create the physics, will set the number of fields and their names */
679: (*r)(&ctx);
680: }
682: /* Create a DMDA to manage the parallel grid */
683: DMDACreate1d(comm, DM_BOUNDARY_PERIODIC, 50, ctx.physics.dof, 2, NULL, &da);
684: DMSetFromOptions(da);
685: 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++) DMDASetFieldName(da, i, ctx.physics.fieldname[i]);
689: DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0);
690: DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0);
692: /* Set coordinates of cell centers */
693: 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: PetscMalloc4(dof * dof, &ctx.R, dof * dof, &ctx.Rinv, 2 * dof, &ctx.cjmpLR, 1 * dof, &ctx.cslope);
697: 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: DMCreateGlobalVector(da, &X);
701: VecDuplicate(X, &X0);
702: VecDuplicate(X, &R);
704: /*-------------------------------- create index for slow parts and fast parts ----------------------------------------*/
705: PetscMalloc1(xm * dof, &index_slow);
706: 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: ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss);
714: ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf);
716: /* Create a time-stepping object */
717: TSCreate(comm, &ts);
718: TSSetDM(ts, da);
719: TSSetRHSFunction(ts, R, FVRHSFunction, &ctx);
720: TSRHSSplitSetIS(ts, "slow", ctx.iss);
721: TSRHSSplitSetIS(ts, "fast", ctx.isf);
722: TSRHSSplitSetRHSFunction(ts, "slow", NULL, FVRHSFunctionslow, &ctx);
723: 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: ISCreateGeneral(PETSC_COMM_WORLD, islow, index_slow, PETSC_COPY_VALUES, &ctx.iss2);
737: ISCreateGeneral(PETSC_COMM_WORLD, ifast, index_fast, PETSC_COPY_VALUES, &ctx.isf2);
739: TSRHSSplitGetSubTS(ts, "fast", &subts);
740: TSRHSSplitSetIS(subts, "slow", ctx.iss2);
741: TSRHSSplitSetIS(subts, "fast", ctx.isf2);
742: TSRHSSplitSetRHSFunction(subts, "slow", NULL, FVRHSFunctionslow2, &ctx);
743: TSRHSSplitSetRHSFunction(subts, "fast", NULL, FVRHSFunctionfast2, &ctx);
744: }
746: /*TSSetType(ts,TSSSP);*/
747: TSSetType(ts, TSMPRK);
748: TSSetMaxTime(ts, 10);
749: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
751: /* Compute initial conditions and starting time step */
752: FVSample(&ctx, da, 0, X0);
753: FVRHSFunction(ts, 0, X0, X, (void *)&ctx); /* Initial function evaluation, only used to determine max speed */
754: VecCopy(X0, X); /* The function value was not used so we set X=X0 again */
755: TSSetTimeStep(ts, ctx.cfl / ctx.cfl_idt);
756: TSSetFromOptions(ts); /* Take runtime options */
757: SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD);
758: {
759: PetscInt steps;
760: PetscScalar mass_initial, mass_final, mass_difference;
762: TSSolve(ts, X);
763: TSGetSolveTime(ts, &ptime);
764: TSGetStepNumber(ts, &steps);
765: 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: VecSum(X0, &mass_initial);
770: VecSum(X, &mass_final);
771: mass_difference = (ctx.xmax - ctx.xmin) / (PetscScalar)Mx * (mass_final - mass_initial);
772: 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: PetscOptionsGetString(NULL, NULL, "-f", filename, sizeof(filename), &flg);
782: PetscViewerBinaryOpen(PETSC_COMM_WORLD, filename, FILE_MODE_READ, &fd);
783: VecDuplicate(X0, &XR);
784: VecLoad(XR, fd);
785: PetscViewerDestroy(&fd);
786: VecAYPX(XR, -1, X);
787: VecNorm(XR, NORM_1, &nrm1);
788: VecNorm(XR, NORM_INFINITY, &nrmsup);
789: nrm1 /= Mx;
790: PetscPrintf(comm, "Error ||x-x_e||_1 %g ||x-x_e||_sup %g\n", (double)nrm1, (double)nrmsup);
791: VecDestroy(&XR);
792: }
793: }
795: SolutionStatsView(da, X, PETSC_VIEWER_STDOUT_WORLD);
796: if (draw & 0x1) VecView(X0, PETSC_VIEWER_DRAW_WORLD);
797: if (draw & 0x2) VecView(X, PETSC_VIEWER_DRAW_WORLD);
798: if (draw & 0x4) {
799: Vec Y;
800: VecDuplicate(X, &Y);
801: FVSample(&ctx, da, ptime, Y);
802: VecAYPX(Y, -1, X);
803: VecView(Y, PETSC_VIEWER_DRAW_WORLD);
804: VecDestroy(&Y);
805: }
807: if (view_final) {
808: PetscViewer viewer;
809: PetscViewerASCIIOpen(PETSC_COMM_WORLD, final_fname, &viewer);
810: PetscViewerPushFormat(viewer, PETSC_VIEWER_ASCII_MATLAB);
811: VecView(X, viewer);
812: PetscViewerPopFormat(viewer);
813: PetscViewerDestroy(&viewer);
814: }
816: /* Clean up */
817: (*ctx.physics.destroy)(ctx.physics.user);
818: for (i = 0; i < ctx.physics.dof; i++) PetscFree(ctx.physics.fieldname[i]);
819: PetscFree4(ctx.R, ctx.Rinv, ctx.cjmpLR, ctx.cslope);
820: PetscFree3(ctx.uLR, ctx.flux, ctx.speeds);
821: ISDestroy(&ctx.iss);
822: ISDestroy(&ctx.isf);
823: ISDestroy(&ctx.iss2);
824: ISDestroy(&ctx.isf2);
825: VecDestroy(&X);
826: VecDestroy(&X0);
827: VecDestroy(&R);
828: DMDestroy(&da);
829: TSDestroy(&ts);
830: PetscFree(index_slow);
831: PetscFree(index_fast);
832: PetscFunctionListDestroy(&limiters);
833: PetscFunctionListDestroy(&physics);
834: 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*/