xref: /petsc/src/ts/tutorials/multirate/finitevolume1d.c (revision 834855d6effb0d027771461c8e947ee1ce5a1e17) 
1 #include "finitevolume1d.h"
2 #include <petscdmda.h>
3 #include <petscdraw.h>
4 #include <petsc/private/tsimpl.h>
5 
6 #include <petsc/private/kernels/blockinvert.h> /* For the Kernel_*_gets_* stuff for BAIJ */
7 const char *FVBCTypes[] = {"PERIODIC", "OUTFLOW", "INFLOW", "FVBCType", "FVBC_", 0};
8 
Sgn(PetscReal a)9 static inline PetscReal Sgn(PetscReal a)
10 {
11   return (a < 0) ? -1 : 1;
12 }
Abs(PetscReal a)13 static inline PetscReal Abs(PetscReal a)
14 {
15   return (a < 0) ? 0 : a;
16 }
Sqr(PetscReal a)17 static inline PetscReal Sqr(PetscReal a)
18 {
19   return a * a;
20 }
21 
MinAbs(PetscReal a,PetscReal b)22 PETSC_UNUSED static inline PetscReal MinAbs(PetscReal a, PetscReal b)
23 {
24   return (PetscAbs(a) < PetscAbs(b)) ? a : b;
25 }
MinMod2(PetscReal a,PetscReal b)26 static inline PetscReal MinMod2(PetscReal a, PetscReal b)
27 {
28   return (a * b < 0) ? 0 : Sgn(a) * PetscMin(PetscAbs(a), PetscAbs(b));
29 }
MaxMod2(PetscReal a,PetscReal b)30 static inline PetscReal MaxMod2(PetscReal a, PetscReal b)
31 {
32   return (a * b < 0) ? 0 : Sgn(a) * PetscMax(PetscAbs(a), PetscAbs(b));
33 }
MinMod3(PetscReal a,PetscReal b,PetscReal c)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 }
38 
39 /* ----------------------- Lots of limiters, these could go in a separate library ------------------------- */
Limit_Upwind(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_LaxWendroff(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_BeamWarming(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_Fromm(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_Minmod(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_Superbee(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_MC(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_VanLeer(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_VanAlbada(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_VanAlbadaTVD(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_Koren(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_KorenSym(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_Koren3(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
CadaTorrilhonPhiHatR_Eq13(PetscReal L,PetscReal R)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 }
Limit_CadaTorrilhon2(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_CadaTorrilhon3R(PetscReal r,LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_CadaTorrilhon3R0p1(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)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 }
Limit_CadaTorrilhon3R1(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)128 void Limit_CadaTorrilhon3R1(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
129 {
130   Limit_CadaTorrilhon3R(1, info, jL, jR, lmt);
131 }
Limit_CadaTorrilhon3R10(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)132 void Limit_CadaTorrilhon3R10(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
133 {
134   Limit_CadaTorrilhon3R(10, info, jL, jR, lmt);
135 }
Limit_CadaTorrilhon3R100(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,PetscScalar * lmt)136 void Limit_CadaTorrilhon3R100(LimitInfo info, const PetscScalar *jL, const PetscScalar *jR, PetscScalar *lmt)
137 {
138   Limit_CadaTorrilhon3R(100, info, jL, jR, lmt);
139 }
140 
141 /* ----------------------- Limiters for split systems ------------------------- */
142 
Limit2_Upwind(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_LaxWendroff(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_BeamWarming(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_Fromm(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_Minmod(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_Superbee(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_MC(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
Limit2_Koren3(LimitInfo info,const PetscScalar * jL,const PetscScalar * jR,const PetscInt sf,const PetscInt fs,PetscInt n,PetscScalar * lmt)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 }
281 
282 /* ---- Three-way splitting ---- */
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 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 }
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 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 }
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 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 }
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 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 }
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 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 }
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 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 }
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 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 }
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 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 }
453 
RiemannListAdd(PetscFunctionList * flist,const char * name,RiemannFunction rsolve)454 PetscErrorCode RiemannListAdd(PetscFunctionList *flist, const char *name, RiemannFunction rsolve)
455 {
456   PetscFunctionBeginUser;
457   PetscCall(PetscFunctionListAdd(flist, name, rsolve));
458   PetscFunctionReturn(PETSC_SUCCESS);
459 }
460 
RiemannListFind(PetscFunctionList flist,const char * name,RiemannFunction * rsolve)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 }
468 
ReconstructListAdd(PetscFunctionList * flist,const char * name,ReconstructFunction r)469 PetscErrorCode ReconstructListAdd(PetscFunctionList *flist, const char *name, ReconstructFunction r)
470 {
471   PetscFunctionBeginUser;
472   PetscCall(PetscFunctionListAdd(flist, name, r));
473   PetscFunctionReturn(PETSC_SUCCESS);
474 }
475 
ReconstructListFind(PetscFunctionList flist,const char * name,ReconstructFunction * r)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 }
483 
RiemannListAdd_2WaySplit(PetscFunctionList * flist,const char * name,RiemannFunction_2WaySplit rsolve)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 }
490 
RiemannListFind_2WaySplit(PetscFunctionList flist,const char * name,RiemannFunction_2WaySplit * rsolve)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 }
498 
ReconstructListAdd_2WaySplit(PetscFunctionList * flist,const char * name,ReconstructFunction_2WaySplit r)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 }
505 
ReconstructListFind_2WaySplit(PetscFunctionList flist,const char * name,ReconstructFunction_2WaySplit * r)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 }
513 
514 /* --------------------------------- Physics ------- */
PhysicsDestroy_SimpleFree(void * vctx)515 PetscErrorCode PhysicsDestroy_SimpleFree(void *vctx)
516 {
517   PetscFunctionBeginUser;
518   PetscCall(PetscFree(vctx));
519   PetscFunctionReturn(PETSC_SUCCESS);
520 }
521 
522 /* --------------------------------- Finite Volume Solver --------------- */
FVRHSFunction(TS ts,PetscReal time,Vec X,Vec F,void * vctx)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;
531 
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));
544 
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   }
553 
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   }
583 
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   PetscCallMPI(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) 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 }
619 
FVSample(FVCtx * ctx,DM da,PetscReal time,Vec U)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;
624 
625   PetscFunctionBeginUser;
626   PetscCheck(ctx->physics.sample, PETSC_COMM_SELF, PETSC_ERR_SUP, "Physics has not provided a sampling function");
627   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
628   PetscCall(DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0));
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 }
646 
SolutionStatsView(DM da,Vec X,PetscViewer viewer)647 PetscErrorCode SolutionStatsView(DM da, Vec X, PetscViewer viewer)
648 {
649   PetscReal          xmin, xmax;
650   PetscScalar        sum, tvsum, tvgsum;
651   const PetscScalar *x;
652   PetscInt           imin, imax, Mx, i, j, xs, xm, dof;
653   Vec                Xloc;
654   PetscBool          isascii;
655 
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, 0, 0, &xm, 0, 0));
665   PetscCall(DMDAGetInfo(da, 0, &Mx, 0, 0, 0, 0, 0, &dof, 0, 0, 0, 0, 0));
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(&tvsum, &tvgsum, 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)(tvgsum / Mx)));
677   PetscFunctionReturn(PETSC_SUCCESS);
678 }
679