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