Actual source code: finitevolume1d.c
1: #include "finitevolume1d.h"
2: #include <petscdmda.h>
3: #include <petscdraw.h>
4: #include <petsc/private/tsimpl.h>
6: #include <petsc/private/kernels/blockinvert.h>
7: const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "INFLOW", "FVBCType", "FVBC_", 0};
9: static inline PetscReal Sgn(PetscReal a)
10: {
11: return (a < 0) ? -1 : 1;
12: }
13: static inline PetscReal Abs(PetscReal a)
14: {
15: return (a < 0) ? 0 : a;
16: }
17: static inline PetscReal Sqr(PetscReal a)
18: {
19: return a * a;
20: }
22: PETSC_UNUSED static inline PetscReal MinAbs(PetscReal a, PetscReal b)
23: {
24: return (PetscAbs(a) < PetscAbs(b)) ? a : b;
25: }
26: static inline PetscReal MinMod2(PetscReal a, PetscReal b)
27: {
28: return (a * b < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscAbs(b));
29: }
30: static inline PetscReal MaxMod2(PetscReal a, PetscReal b)
31: {
32: return (a * b < 0) ? 0 : Sgn(a) * PetscMax(PetscAbs(a), PetscAbs(b));
33: }
34: static inline PetscReal MinMod3(PetscReal a, PetscReal b, PetscReal c)
35: {
36: return (a * b < 0 || a * c < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscMin(PetscAbs(b), PetscAbs(c)));
37: }
39: /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
40: void Limit_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
41: {
42: PetscInt i;
43: for (i = 0; i < info->m; i++) lmt[i] = 0;
44: }
45: void Limit_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
46: {
47: PetscInt i;
48: for (i = 0; i < info->m; i++) lmt[i] = jR[i];
49: }
50: void Limit_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
51: {
52: PetscInt i;
53: for (i = 0; i < info->m; i++) lmt[i] = jL[i];
54: }
55: void Limit_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
56: {
57: PetscInt i;
58: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]);
59: }
60: void Limit_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
61: {
62: PetscInt i;
63: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]);
64: }
65: void Limit_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
66: {
67: PetscInt i;
68: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i]));
69: }
70: void Limit_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
71: {
72: PetscInt i;
73: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]);
74: }
75: void Limit_VanLeer(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
76: { /* phi = (t + abs(t)) / (1 + abs(t)) */
77: PetscInt i;
78: for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * Abs(jR[i]) + Abs(jL[i]) * jR[i]) / (Abs(jL[i]) + Abs(jR[i]) + 1e-15);
79: }
80: void Limit_VanAlbada(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
81: { /* phi = (t + t^2) / (1 + t^2) */
82: PetscInt i;
83: for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
84: }
85: void Limit_VanAlbadaTVD(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
86: { /* phi = (t + t^2) / (1 + t^2) */
87: PetscInt i;
88: for (i = 0; i < info->m; i++) lmt[i] = (jL[i] * jR[i] < 0) ? 0 : (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (Sqr(jL[i]) + Sqr(jR[i]) + 1e-15);
89: }
90: void Limit_Koren(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
91: { /* phi = (t + 2*t^2) / (2 - t + 2*t^2) */
92: PetscInt i;
93: for (i = 0; i < info->m; i++) lmt[i] = ((jL[i] * Sqr(jR[i]) + 2 * Sqr(jL[i]) * jR[i]) / (2 * Sqr(jL[i]) - jL[i] * jR[i] + 2 * Sqr(jR[i]) + 1e-15));
94: }
95: void Limit_KorenSym(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt) /* differentiable */
96: { /* Symmetric version of above */
97: PetscInt i;
98: for (i = 0; i < info->m; i++) lmt[i] = (1.5 * (jL[i] * Sqr(jR[i]) + Sqr(jL[i]) * jR[i]) / (2 * Sqr(jL[i]) - jL[i] * jR[i] + 2 * Sqr(jR[i]) + 1e-15));
99: }
100: void Limit_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
101: { /* Eq 11 of Cada-Torrilhon 2009 */
102: PetscInt i;
103: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]);
104: }
105: static PetscReal CadaTorrilhonPhiHatR_Eq13(PetscReal L, PetscReal R)
106: {
107: return PetscMax(0, PetscMin((L + 2 * R) / 3, PetscMax(-0.5 * L, PetscMin(2 * L, PetscMin((L + 2 * R) / 3, 1.6 * R)))));
108: }
109: void Limit_CadaTorrilhon2(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
110: { /* Cada-Torrilhon 2009, Eq 13 */
111: PetscInt i;
112: for (i = 0; i < info->m; i++) lmt[i] = CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]);
113: }
114: void Limit_CadaTorrilhon3R(PetscReal r, LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
115: { /* Cada-Torrilhon 2009, Eq 22 */
116: /* They recommend 0.001 < r < 1, but larger values are more accurate in smooth regions */
117: const PetscReal eps = 1e-7, hx = info->hx;
118: PetscInt i;
119: for (i = 0; i < info->m; i++) {
120: const PetscReal eta = (Sqr(jL[i]) + Sqr(jR[i])) / Sqr(r * hx);
121: lmt[i] = ((eta < 1 - eps) ? (jL[i] + 2 * jR[i]) / 3 : ((eta > 1 + eps) ? CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]) : 0.5 * ((1 - (eta - 1) / eps) * (jL[i] + 2 * jR[i]) / 3 + (1 + (eta + 1) / eps) * CadaTorrilhonPhiHatR_Eq13(jL[i], jR[i]))));
122: }
123: }
124: void Limit_CadaTorrilhon3R0p1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
125: {
126: Limit_CadaTorrilhon3R(0.1, info, jL, jR, lmt);
127: }
128: void Limit_CadaTorrilhon3R1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
129: {
130: Limit_CadaTorrilhon3R(1, info, jL, jR, lmt);
131: }
132: void Limit_CadaTorrilhon3R10(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
133: {
134: Limit_CadaTorrilhon3R(10, info, jL, jR, lmt);
135: }
136: void Limit_CadaTorrilhon3R100(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
137: {
138: Limit_CadaTorrilhon3R(100, info, jL, jR, lmt);
139: }
141: /* ----------------------- Limiters for split systems ------------------------- */
143: void Limit2_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
144: {
145: PetscInt i;
146: for (i = 0; i < info->m; i++) lmt[i] = 0;
147: }
148: void Limit2_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
149: {
150: PetscInt i;
151: if (n < sf - 1) { /* slow components */
152: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
153: } else if (n == sf - 1) { /* slow component which is next to fast components */
154: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxs / 2.0 + info->hxf / 2.0);
155: } else if (n == sf) { /* fast component which is next to slow components */
156: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxf;
157: } else if (n > sf && n < fs - 1) { /* fast components */
158: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxf;
159: } else if (n == fs - 1) { /* fast component next to slow components */
160: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxf / 2.0 + info->hxs / 2.0);
161: } else if (n == fs) { /* slow component next to fast components */
162: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
163: } else { /* slow components */
164: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
165: }
166: }
167: void Limit2_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
168: {
169: PetscInt i;
170: if (n < sf - 1) {
171: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
172: } else if (n == sf - 1) {
173: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
174: } else if (n == sf) {
175: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxs / 2.0 + info->hxf / 2.0);
176: } else if (n > sf && n < fs - 1) {
177: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxf;
178: } else if (n == fs - 1) {
179: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxf;
180: } else if (n == fs) {
181: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxf / 2.0 + info->hxs / 2.0);
182: } else {
183: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
184: }
185: }
186: void Limit2_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
187: {
188: PetscInt i;
189: if (n < sf - 1) {
190: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxs;
191: } else if (n == sf - 1) {
192: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
193: } else if (n == sf) {
194: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf);
195: } else if (n > sf && n < fs - 1) {
196: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxf;
197: } else if (n == fs - 1) {
198: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
199: } else if (n == fs) {
200: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs);
201: } else {
202: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxs;
203: }
204: }
205: void Limit2_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
206: {
207: PetscInt i;
208: if (n < sf - 1) {
209: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxs;
210: } else if (n == sf - 1) {
211: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
212: } else if (n == sf) {
213: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), jR[i] / info->hxf);
214: } else if (n > sf && n < fs - 1) {
215: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxf;
216: } else if (n == fs - 1) {
217: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
218: } else if (n == fs) {
219: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), jR[i] / info->hxs);
220: } else {
221: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxs;
222: }
223: }
224: void Limit2_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
225: {
226: PetscInt i;
227: if (n < sf - 1) {
228: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxs;
229: } else if (n == sf - 1) {
230: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / info->hxs, 2 * jR[i] / (info->hxs / 2.0 + info->hxf / 2.0)), MinMod2(2 * jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxf / 2.0)));
231: } else if (n == sf) {
232: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), 2 * jR[i] / info->hxf), MinMod2(2 * jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), jR[i] / info->hxf));
233: } else if (n > sf && n < fs - 1) {
234: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxf;
235: } else if (n == fs - 1) {
236: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / info->hxf, 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)), MinMod2(2 * jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)));
237: } else if (n == fs) {
238: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), 2 * jR[i] / info->hxs), MinMod2(2 * jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), jR[i] / info->hxs));
239: } else {
240: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxs;
241: }
242: }
243: void Limit2_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
244: {
245: PetscInt i;
246: if (n < sf - 1) {
247: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxs;
248: } else if (n == sf - 1) {
249: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)), 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
250: } else if (n == sf) {
251: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf), 2 * jR[i] / info->hxf);
252: } else if (n > sf && n < fs - 1) {
253: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxf;
254: } else if (n == fs - 1) {
255: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)), 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
256: } else if (n == fs) {
257: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs), 2 * jR[i] / info->hxs);
258: } else {
259: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxs;
260: }
261: }
262: void Limit2_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sf, const PetscInt fs, PetscInt n, PetscScalar *lmt)
263: { /* Eq 11 of Cada-Torrilhon 2009 */
264: PetscInt i;
265: if (n < sf - 1) {
266: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
267: } else if (n == sf - 1) {
268: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, (jL[i] / info->hxs + 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
269: } else if (n == sf) {
270: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), (jL[i] / (info->hxs / 2.0 + info->hxf / 2.0) + 2 * jR[i] / info->hxf) / 3, 2 * jR[i] / info->hxf);
271: } else if (n > sf && n < fs - 1) {
272: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxf;
273: } else if (n == fs - 1) {
274: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, (jL[i] / info->hxf + 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxf / 2.0 + info->hxs / 2.0));
275: } else if (n == fs) {
276: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxs / 2.0), (jL[i] / (info->hxf / 2.0 + info->hxs / 2.0) + 2 * jR[i] / info->hxs) / 3, 2 * jR[i] / info->hxs);
277: } else {
278: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
279: }
280: }
282: /* ---- Three-way splitting ---- */
283: void Limit3_Upwind(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
284: {
285: PetscInt i;
286: for (i = 0; i < info->m; i++) lmt[i] = 0;
287: }
288: void Limit3_LaxWendroff(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
289: {
290: PetscInt i;
291: if (n < sm - 1 || n > ms) { /* slow components */
292: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxs;
293: } else if (n == sm - 1 || n == ms - 1) { /* slow component which is next to medium components */
294: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxs / 2.0 + info->hxm / 2.0);
295: } else if (n < mf - 1 || n > fm) { /* medium components */
296: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxm;
297: } else if (n == mf - 1 || n == fm) { /* medium component next to fast components */
298: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / (info->hxm / 2.0 + info->hxf / 2.0);
299: } else { /* fast components */
300: for (i = 0; i < info->m; i++) lmt[i] = jR[i] / info->hxf;
301: }
302: }
303: void Limit3_BeamWarming(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
304: {
305: PetscInt i;
306: if (n < sm || n > ms) {
307: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxs;
308: } else if (n == sm || n == ms) {
309: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxs / 2.0 + info->hxf / 2.0);
310: } else if (n < mf || n > fm) {
311: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxm;
312: } else if (n == mf || n == fm) {
313: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / (info->hxm / 2.0 + info->hxf / 2.0);
314: } else {
315: for (i = 0; i < info->m; i++) lmt[i] = jL[i] / info->hxf;
316: }
317: }
318: void Limit3_Fromm(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
319: {
320: PetscInt i;
321: if (n < sm - 1 || n > ms) {
322: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxs;
323: } else if (n == sm - 1) {
324: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
325: } else if (n == sm) {
326: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm);
327: } else if (n == ms - 1) {
328: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
329: } else if (n == ms) {
330: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs);
331: } else if (n < mf - 1 || n > fm) {
332: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxm;
333: } else if (n == mf - 1) {
334: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
335: } else if (n == mf) {
336: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf);
337: } else if (n == fm - 1) {
338: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
339: } else if (n == fm) {
340: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm);
341: } else {
342: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxf;
343: }
344: }
345: void Limit3_Minmod(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
346: {
347: PetscInt i;
348: if (n < sm - 1 || n > ms) {
349: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxs;
350: } else if (n == sm - 1) {
351: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxf / 2.0));
352: } else if (n == sm) {
353: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxs / 2.0 + info->hxf / 2.0), jR[i] / info->hxf);
354: } else if (n == ms - 1) {
355: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
356: } else if (n == ms) {
357: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), jR[i] / info->hxs);
358: } else if (n < mf - 1 || n > fm) {
359: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i], jR[i]) / info->hxm;
360: } else if (n == mf - 1) {
361: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
362: } else if (n == mf) {
363: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), jR[i] / info->hxf);
364: } else if (n == fm - 1) {
365: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
366: } else if (n == fm) {
367: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), jR[i] / info->hxm);
368: } else {
369: for (i = 0; i < info->m; i++) lmt[i] = 0.5 * (jL[i] + jR[i]) / info->hxf;
370: }
371: }
372: void Limit3_Superbee(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
373: {
374: PetscInt i;
375: if (n < sm - 1 || n > ms) {
376: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxs;
377: } else if (n == sm - 1) {
378: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / info->hxs, 2 * jR[i] / (info->hxs / 2.0 + info->hxm / 2.0)), MinMod2(2 * jL[i] / info->hxs, jR[i] / (info->hxs / 2.0 + info->hxm / 2.0)));
379: } else if (n == sm) {
380: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), 2 * jR[i] / info->hxm), MinMod2(2 * jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), jR[i] / info->hxm));
381: } else if (n == ms - 1) {
382: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / info->hxm, 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)), MinMod2(2 * jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)));
383: } else if (n == ms) {
384: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), 2 * jR[i] / info->hxs), MinMod2(2 * jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), jR[i] / info->hxs));
385: } else if (n < mf - 1 || n > fm) {
386: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxm;
387: } else if (n == mf - 1) {
388: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / info->hxm, 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)), MinMod2(2 * jL[i] / info->hxm, jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)));
389: } else if (n == mf) {
390: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), 2 * jR[i] / info->hxf), MinMod2(2 * jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), jR[i] / info->hxf));
391: } else if (n == fm - 1) {
392: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / info->hxf, 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)), MinMod2(2 * jL[i] / info->hxf, jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)));
393: } else if (n == fm) {
394: for (i = 0; i < info->m; i++) lmt[i] = MinMod2(MinMod2(jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), 2 * jR[i] / info->hxm), MinMod2(2 * jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), jR[i] / info->hxm));
395: } else {
396: for (i = 0; i < info->m; i++) lmt[i] = MaxMod2(MinMod2(jL[i], 2 * jR[i]), MinMod2(2 * jL[i], jR[i])) / info->hxf;
397: }
398: }
399: void Limit3_MC(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
400: {
401: PetscInt i;
402: if (n < sm - 1 || n > ms) {
403: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxs;
404: } else if (n == sm - 1) {
405: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, 0.5 * (jL[i] / info->hxs + jR[i] / (info->hxs / 2.0 + info->hxm / 2.0)), 2 * jR[i] / (info->hxs / 2.0 + info->hxm / 2.0));
406: } else if (n == sm) {
407: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), 0.5 * (jL[i] / (info->hxs / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm), 2 * jR[i] / info->hxm);
408: } else if (n == ms - 1) {
409: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)), 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
410: } else if (n == ms) {
411: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxs / 2.0) + jR[i] / info->hxs), 2 * jR[i] / info->hxs);
412: } else if (n < mf - 1 || n > fm) {
413: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxm;
414: } else if (n == mf - 1) {
415: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, 0.5 * (jL[i] / info->hxm + jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)), 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
416: } else if (n == mf) {
417: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), 0.5 * (jL[i] / (info->hxm / 2.0 + info->hxf / 2.0) + jR[i] / info->hxf), 2 * jR[i] / info->hxf);
418: } else if (n == fm - 1) {
419: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, 0.5 * (jL[i] / info->hxf + jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)), 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
420: } else if (n == fm) {
421: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), 0.5 * (jL[i] / (info->hxf / 2.0 + info->hxm / 2.0) + jR[i] / info->hxm), 2 * jR[i] / info->hxm);
422: } else {
423: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], 0.5 * (jL[i] + jR[i]), 2 * jR[i]) / info->hxf;
424: }
425: }
426: void Limit3_Koren3(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, const PetscInt sm, const PetscInt mf, const PetscInt fm, const PetscInt ms, PetscInt n, PetscScalar *lmt)
427: { /* Eq 11 of Cada-Torrilhon 2009 */
428: PetscInt i;
429: if (n < sm - 1 || n > ms) {
430: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
431: } else if (n == sm - 1) {
432: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxs, (jL[i] / info->hxs + 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
433: } else if (n == sm) {
434: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxs / 2.0 + info->hxm / 2.0), (jL[i] / (info->hxs / 2.0 + info->hxm / 2.0) + 2 * jR[i] / info->hxm) / 3, 2 * jR[i] / info->hxm);
435: } else if (n == ms - 1) {
436: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, (jL[i] / info->hxm + 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0)) / 3, 2 * jR[i] / (info->hxm / 2.0 + info->hxs / 2.0));
437: } else if (n == ms) {
438: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxs / 2.0), (jL[i] / (info->hxm / 2.0 + info->hxs / 2.0) + 2 * jR[i] / info->hxs) / 3, 2 * jR[i] / info->hxs);
439: } else if (n < mf - 1 || n > fm) {
440: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxm;
441: } else if (n == mf - 1) {
442: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxm, (jL[i] / info->hxm + 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0)) / 3, 2 * jR[i] / (info->hxm / 2.0 + info->hxf / 2.0));
443: } else if (n == mf) {
444: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxm / 2.0 + info->hxf / 2.0), (jL[i] / (info->hxm / 2.0 + info->hxf / 2.0) + 2 * jR[i] / info->hxf) / 3, 2 * jR[i] / info->hxf);
445: } else if (n == fm - 1) {
446: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / info->hxf, (jL[i] / info->hxf + 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0)) / 3, 2 * jR[i] / (info->hxf / 2.0 + info->hxm / 2.0));
447: } else if (n == fm) {
448: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i] / (info->hxf / 2.0 + info->hxm / 2.0), (jL[i] / (info->hxf / 2.0 + info->hxm / 2.0) + 2 * jR[i] / info->hxm) / 3, 2 * jR[i] / info->hxm);
449: } else {
450: for (i = 0; i < info->m; i++) lmt[i] = MinMod3(2 * jL[i], (jL[i] + 2 * jR[i]) / 3, 2 * jR[i]) / info->hxs;
451: }
452: }
454: PetscErrorCode RiemannListAdd(PetscFunctionList *flist, const char *name, RiemannFunction rsolve)
455: {
456: PetscFunctionBeginUser;
457: PetscCall(PetscFunctionListAdd(flist, name, rsolve));
458: PetscFunctionReturn(PETSC_SUCCESS);
459: }
461: PetscErrorCode RiemannListFind(PetscFunctionList flist, const char *name, RiemannFunction *rsolve)
462: {
463: PetscFunctionBeginUser;
464: PetscCall(PetscFunctionListFind(flist, name, rsolve));
465: PetscCheck(*rsolve, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Riemann solver \"%s\" could not be found", name);
466: PetscFunctionReturn(PETSC_SUCCESS);
467: }
469: PetscErrorCode ReconstructListAdd(PetscFunctionList *flist, const char *name, ReconstructFunction r)
470: {
471: PetscFunctionBeginUser;
472: PetscCall(PetscFunctionListAdd(flist, name, r));
473: PetscFunctionReturn(PETSC_SUCCESS);
474: }
476: PetscErrorCode ReconstructListFind(PetscFunctionList flist, const char *name, ReconstructFunction *r)
477: {
478: PetscFunctionBeginUser;
479: PetscCall(PetscFunctionListFind(flist, name, r));
480: PetscCheck(*r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Reconstruction \"%s\" could not be found", name);
481: PetscFunctionReturn(PETSC_SUCCESS);
482: }
484: PetscErrorCode RiemannListAdd_2WaySplit(PetscFunctionList *flist, const char *name, RiemannFunction_2WaySplit rsolve)
485: {
486: PetscFunctionBeginUser;
487: PetscCall(PetscFunctionListAdd(flist, name, rsolve));
488: PetscFunctionReturn(PETSC_SUCCESS);
489: }
491: PetscErrorCode RiemannListFind_2WaySplit(PetscFunctionList flist, const char *name, RiemannFunction_2WaySplit *rsolve)
492: {
493: PetscFunctionBeginUser;
494: PetscCall(PetscFunctionListFind(flist, name, rsolve));
495: PetscCheck(*rsolve, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Riemann solver \"%s\" could not be found", name);
496: PetscFunctionReturn(PETSC_SUCCESS);
497: }
499: PetscErrorCode ReconstructListAdd_2WaySplit(PetscFunctionList *flist, const char *name, ReconstructFunction_2WaySplit r)
500: {
501: PetscFunctionBeginUser;
502: PetscCall(PetscFunctionListAdd(flist, name, r));
503: PetscFunctionReturn(PETSC_SUCCESS);
504: }
506: PetscErrorCode ReconstructListFind_2WaySplit(PetscFunctionList flist, const char *name, ReconstructFunction_2WaySplit *r)
507: {
508: PetscFunctionBeginUser;
509: PetscCall(PetscFunctionListFind(flist, name, r));
510: PetscCheck(*r, PETSC_COMM_SELF, PETSC_ERR_ARG_UNKNOWN_TYPE, "Reconstruction \"%s\" could not be found", name);
511: PetscFunctionReturn(PETSC_SUCCESS);
512: }
514: /* --------------------------------- Physics ------- */
515: PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
516: {
517: PetscFunctionBeginUser;
518: PetscCall(PetscFree(vctx));
519: PetscFunctionReturn(PETSC_SUCCESS);
520: }
522: /* --------------------------------- Finite Volume Solver --------------- */
523: PetscErrorCode FVRHSFunction(TS ts, PetscReal time, Vec X, Vec F, void *vctx)
524: {
525: FVCtx *ctx = (FVCtx *)vctx;
526: PetscInt i, j, k, Mx, dof, xs, xm;
527: PetscReal hx, cfl_idt = 0;
528: PetscScalar *x, *f, *slope;
529: Vec Xloc;
530: DM da;
532: PetscFunctionBeginUser;
533: PetscCall(TSGetDM(ts, &da));
534: PetscCall(DMGetLocalVector(da, &Xloc)); /* Xloc contains ghost points */
535: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0)); /* Mx is the number of center points */
536: hx = (ctx->xmax - ctx->xmin) / Mx;
537: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc)); /* X is solution vector which does not contain ghost points */
538: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
539: PetscCall(VecZeroEntries(F)); /* F is the right-hand side function corresponds to center points */
540: PetscCall(DMDAVecGetArray(da, Xloc, &x));
541: PetscCall(DMDAVecGetArray(da, F, &f));
542: PetscCall(DMDAGetArray(da, PETSC_TRUE, &slope)); /* contains ghost points */
543: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
545: if (ctx->bctype == FVBC_OUTFLOW) {
546: for (i = xs - 2; i < 0; i++) {
547: for (j = 0; j < dof; j++) x[i * dof + j] = x[j];
548: }
549: for (i = Mx; i < xs + xm + 2; i++) {
550: for (j = 0; j < dof; j++) x[i * dof + j] = x[(xs + xm - 1) * dof + j];
551: }
552: }
554: for (i = xs - 1; i < xs + xm + 1; i++) {
555: struct _LimitInfo info;
556: PetscScalar *cjmpL, *cjmpR;
557: /* Determine the right eigenvectors R, where A = R \Lambda R^{-1} */
558: PetscCall((*ctx->physics.characteristic)(ctx->physics.user, dof, &x[i * dof], ctx->R, ctx->Rinv, ctx->speeds, ctx->xmin + hx * i));
559: /* Evaluate jumps across interfaces (i-1, i) and (i, i+1), put in characteristic basis */
560: PetscCall(PetscArrayzero(ctx->cjmpLR, 2 * dof));
561: cjmpL = &ctx->cjmpLR[0];
562: cjmpR = &ctx->cjmpLR[dof];
563: for (j = 0; j < dof; j++) {
564: PetscScalar jmpL, jmpR;
565: jmpL = x[(i + 0) * dof + j] - x[(i - 1) * dof + j];
566: jmpR = x[(i + 1) * dof + j] - x[(i + 0) * dof + j];
567: for (k = 0; k < dof; k++) {
568: cjmpL[k] += ctx->Rinv[k + j * dof] * jmpL;
569: cjmpR[k] += ctx->Rinv[k + j * dof] * jmpR;
570: }
571: }
572: /* Apply limiter to the left and right characteristic jumps */
573: info.m = dof;
574: info.hx = hx;
575: (*ctx->limit)(&info, cjmpL, cjmpR, ctx->cslope);
576: for (j = 0; j < dof; j++) ctx->cslope[j] /= hx; /* rescale to a slope */
577: for (j = 0; j < dof; j++) {
578: PetscScalar tmp = 0;
579: for (k = 0; k < dof; k++) tmp += ctx->R[j + k * dof] * ctx->cslope[k];
580: slope[i * dof + j] = tmp;
581: }
582: }
584: for (i = xs; i < xs + xm + 1; i++) {
585: PetscReal maxspeed;
586: PetscScalar *uL, *uR;
587: uL = &ctx->uLR[0];
588: uR = &ctx->uLR[dof];
589: for (j = 0; j < dof; j++) {
590: uL[j] = x[(i - 1) * dof + j] + slope[(i - 1) * dof + j] * hx / 2;
591: uR[j] = x[(i - 0) * dof + j] - slope[(i - 0) * dof + j] * hx / 2;
592: }
593: PetscCall((*ctx->physics.riemann)(ctx->physics.user, dof, uL, uR, ctx->flux, &maxspeed, ctx->xmin + hx * i, ctx->xmin, ctx->xmax));
594: cfl_idt = PetscMax(cfl_idt, PetscAbsScalar(maxspeed / hx)); /* Max allowable value of 1/Delta t */
595: if (i > xs) {
596: for (j = 0; j < dof; j++) f[(i - 1) * dof + j] -= ctx->flux[j] / hx;
597: }
598: if (i < xs + xm) {
599: for (j = 0; j < dof; j++) f[i * dof + j] += ctx->flux[j] / hx;
600: }
601: }
602: PetscCall(DMDAVecRestoreArray(da, Xloc, &x));
603: PetscCall(DMDAVecRestoreArray(da, F, &f));
604: PetscCall(DMDARestoreArray(da, PETSC_TRUE, &slope));
605: PetscCall(DMRestoreLocalVector(da, &Xloc));
606: PetscCall(MPIU_Allreduce(&cfl_idt, &ctx->cfl_idt, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)da)));
607: if (0) {
608: /* We need to a way to inform the TS of a CFL constraint, this is a debugging fragment */
609: PetscReal dt, tnow;
610: PetscCall(TSGetTimeStep(ts, &dt));
611: PetscCall(TSGetTime(ts, &tnow));
612: if (dt > 0.5 / ctx->cfl_idt) {
613: if (1) {
614: PetscCall(PetscPrintf(ctx->comm, "Stability constraint exceeded at t=%g, dt %g > %g\n", (double)tnow, (double)dt, (double)(0.5 / ctx->cfl_idt)));
615: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Stability constraint exceeded, %g > %g", (double)dt, (double)(ctx->cfl / ctx->cfl_idt));
616: }
617: }
618: PetscFunctionReturn(PETSC_SUCCESS);
619: }
621: PetscErrorCode FVSample(FVCtx *ctx, DM da, PetscReal time, Vec U)
622: {
623: PetscScalar *u, *uj;
624: PetscInt i, j, k, dof, xs, xm, Mx;
626: PetscFunctionBeginUser;
627: PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
628: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
629: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
630: PetscCall(DMDAVecGetArray(da, U, &u));
631: PetscCall(PetscMalloc1(dof, &uj));
632: for (i = xs; i < xs + xm; i++) {
633: const PetscReal h = (ctx->xmax - ctx->xmin) / Mx, xi = ctx->xmin + h / 2 + i * h;
634: const PetscInt N = 200;
635: /* Integrate over cell i using trapezoid rule with N points. */
636: for (k = 0; k < dof; k++) u[i * dof + k] = 0;
637: for (j = 0; j < N + 1; j++) {
638: PetscScalar xj = xi + h * (j - N / 2) / (PetscReal)N;
639: PetscCall((*ctx->physics.sample)(ctx->physics.user, ctx->initial, ctx->bctype, ctx->xmin, ctx->xmax, time, xj, uj));
640: for (k = 0; k < dof; k++) u[i * dof + k] += ((j == 0 || j == N) ? 0.5 : 1.0) * uj[k] / N;
641: }
642: }
643: PetscCall(DMDAVecRestoreArray(da, U, &u));
644: PetscCall(PetscFree(uj));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
649: {
650: PetscReal xmin, xmax;
651: PetscScalar sum, tvsum, tvgsum;
652: const PetscScalar *x;
653: PetscInt imin, imax, Mx, i, j, xs, xm, dof;
654: Vec Xloc;
655: PetscBool iascii;
657: PetscFunctionBeginUser;
658: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
659: if (iascii) {
660: /* PETSc lacks a function to compute total variation norm (difficult in multiple dimensions), we do it here */
661: PetscCall(DMGetLocalVector(da, &Xloc));
662: PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, Xloc));
663: PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, Xloc));
664: PetscCall(DMDAVecGetArrayRead(da, Xloc, (void *)&x));
665: PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
666: PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
667: tvsum = 0;
668: for (i = xs; i < xs + xm; i++) {
669: for (j = 0; j < dof; j++) tvsum += PetscAbsScalar(x[i * dof + j] - x[(i - 1) * dof + j]);
670: }
671: PetscCall(MPIU_Allreduce(&tvsum, &tvgsum, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)da)));
672: PetscCall(DMDAVecRestoreArrayRead(da, Xloc, (void *)&x));
673: PetscCall(DMRestoreLocalVector(da, &Xloc));
674: PetscCall(VecMin(X, &imin, &xmin));
675: PetscCall(VecMax(X, &imax, &xmax));
676: PetscCall(VecSum(X, &sum));
677: PetscCall(PetscViewerASCIIPrintf(viewer, "Solution range [%8.5f,%8.5f] with minimum at %" PetscInt_FMT ", mean %8.5f, ||x||_TV %8.5f\n", (double)xmin, (double)xmax, imin, (double)(sum / Mx), (double)(tvgsum / Mx)));
678: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "Viewer type not supported");
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }