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