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