include/pops/numerics/fv/numerical_flux.hpp Source File

adc_cpp 0.3.0 Model-free C++23 core for coupled hyperbolic-elliptic systems on adaptive (AMR) meshes, with MPI and GPU (Kokkos) backends

• include
• pops
• numerics
• fv

numerical_flux.hpp

Go to the documentation of this file.

1

38

39#pragma once

40

41#include <pops/core/state/state.hpp>

42#include <pops/core/foundation/types.hpp>

43

44#include <cmath> // std::sqrt (RoeFlux: Roe average); libstdc++ does not pull it transitively

45#include <concepts> // Riemann capabilities (HasHLLCStructure / HasRoeDissipation)

46

47namespace pops {

48

55struct RusanovFlux {

56 template <class Model>

57 POPS_HD typename Model::State operator()(const Model& m, const typename Model::State& UL,

58 const Aux& AL, const typename Model::State& UR,

59 const Aux& AR, int dir) const {

60 const auto FL = m.flux(UL, AL, dir);

61 const auto FR = m.flux(UR, AR, dir);

62 const Real sL = m.max_wave_speed(UL, AL, dir);

63 const Real sR = m.max_wave_speed(UR, AR, dir);

64 const Real alpha = sL > sR ? sL : sR; // device-safe max (no std::max)

65 typename Model::State F;

66 for (int c = 0; c < Model::n_vars; ++c)

67 F[c] = Real(0.5) * (FL[c] + FR[c]) - Real(0.5) * alpha * (UR[c] - UL[c]);

68 return F;

69 }

70};

71

77template <class Model>

78POPS_HD inline void hll_speeds(const Model& m, const typename Model::State& UL, const Aux& AL,

79 const typename Model::State& UR, const Aux& AR, int dir, Real& sL,

80 Real& sR) {

81 Real lL, hL, lR, hR;

82 m.wave_speeds(UL, AL, dir, lL, hL);

83 m.wave_speeds(UR, AR, dir, lR, hR);

84 sL = lL < lR ? lL : lR;

85 sR = hL > hR ? hL : hR;

86}

87

95template <class Model>

96POPS_HD inline typename Model::State hll_flux_with_speeds(const Model& m,

97 const typename Model::State& UL,

98 const Aux& AL,

99 const typename Model::State& UR,

100 const Aux& AR, int dir, Real sL, Real sR) {

101 const auto FL = m.flux(UL, AL, dir);

102 const auto FR = m.flux(UR, AR, dir);

103 if (sL >= 0)

104 return FL;

105 if (sR <= 0)

106 return FR;

107 typename Model::State F;

108 const Real inv = Real(1) / (sR - sL);

109 for (int c = 0; c < Model::n_vars; ++c)

110 F[c] = (sR * FL[c] - sL * FR[c] + sL * sR * (UR[c] - UL[c])) * inv;

111 return F;

112}

113

119struct HLLFlux {

120 template <class Model>

121 POPS_HD typename Model::State operator()(const Model& m, const typename Model::State& UL,

122 const Aux& AL, const typename Model::State& UR,

123 const Aux& AR, int dir) const {

124 Real sL, sR;

125 hll_speeds(m, UL, AL, UR, AR, dir, sL, sR);

126 return hll_flux_with_speeds(m, UL, AL, UR, AR, dir, sL, sR);

127 }

128};

129

130// ---------------------------------------------------------------------------------------------

131// RIEMANN CAPABILITIES: OPTIONAL traits through which a NON-Euler model supplies the physical

132// structure required by a contact-resolving (HLLC) or Roe-like solver. cf. the file header.

133// ---------------------------------------------------------------------------------------------

134

141template <class M>

142concept HasHLLCStructure =

143 requires(const M m, const typename M::State u, const typename M::State v, const Aux a, Real p,

144 Real q, Real sl, Real sr, Real ss, int dir) {

145 { m.pressure(u) } -> std::convertible_to<Real>;

146 m.wave_speeds(u, a, dir, sl, sr); // signed speeds (hll_speeds, outer wave bounds)

147 { m.contact_speed(u, v, p, q, sl, sr, dir) } -> std::convertible_to<Real>;

148 { m.hllc_star_state(u, p, sl, ss, dir) } -> std::same_as<typename M::State>;

149 };

150

155template <class M>

156concept HasRoeDissipation = requires(const M m, const typename M::State ul, const Aux al,

157 const typename M::State ur, const Aux ar, int dir) {

158 { m.roe_dissipation(ul, al, ur, ar, dir) } -> std::same_as<typename M::State>;

159};

160

169struct HLLCFlux {

170 template <class Model>

171 POPS_HD typename Model::State operator()(const Model& m, const typename Model::State& UL,

172 const Aux& AL, const typename Model::State& UR,

173 const Aux& AR, int dir) const {

174 // CAPABILITY PATH (HasHLLCStructure): GENERIC HLLC algorithm -- the contact speed and the star

175 // states come from the MODEL, the core assumes neither layout nor EOS. A canonical Euler model

176 // WITHOUT hooks takes the historical branch below, bit-identical.

177 if constexpr (HasHLLCStructure<Model>) {

178 Real sL, sR;

179 hll_speeds(m, UL, AL, UR, AR, dir, sL, sR);

180 const auto FL = m.flux(UL, AL, dir);

181 const auto FR = m.flux(UR, AR, dir);

182 if (sL >= 0)

183 return FL;

184 if (sR <= 0)

185 return FR;

186 const Real pL = m.pressure(UL), pR = m.pressure(UR);

187 const Real sStar = m.contact_speed(UL, UR, pL, pR, sL, sR, dir);

188 typename Model::State F;

189 if (sStar >= 0) {

190 const typename Model::State Us = m.hllc_star_state(UL, pL, sL, sStar, dir);

191 for (int c = 0; c < Model::n_vars; ++c)

192 F[c] = FL[c] + sL * (Us[c] - UL[c]);

193 } else {

194 const typename Model::State Us = m.hllc_star_state(UR, pR, sR, sStar, dir);

195 for (int c = 0; c < Model::n_vars; ++c)

196 F[c] = FR[c] + sR * (Us[c] - UR[c]);

197 }

198 return F;

199 } else {

200 const int in = (dir == 0) ? 1 : 2; // normal momentum component

201 const int it = (dir == 0) ? 2 : 1; // tangential

202 const Real rL = UL[0], rR = UR[0];

203 const Real unL = UL[in] / rL, unR = UR[in] / rR;

204 const Real pL = m.pressure(UL), pR = m.pressure(UR);

205 Real sL, sR;

206 hll_speeds(m, UL, AL, UR, AR, dir, sL, sR);

207 const auto FL = m.flux(UL, AL, dir);

208 const auto FR = m.flux(UR, AR, dir);

209 if (sL >= 0)

210 return FL;

211 if (sR <= 0)

212 return FR;

213

214 // contact wave speed (Toro 10.37)

215 const Real sStar = (pR - pL + rL * unL * (sL - unL) - rR * unR * (sR - unR)) /

216 (rL * (sL - unL) - rR * (sR - unR));

217 typename Model::State F;

218 if (sStar >= 0) { // left star state

219 const Real fac = rL * (sL - unL) / (sL - sStar);

220 typename Model::State Us;

221 Us[0] = fac;

222 Us[in] = fac * sStar;

223 Us[it] = fac * (UL[it] / rL);

224 Us[3] = fac * (UL[3] / rL + (sStar - unL) * (sStar + pL / (rL * (sL - unL))));

225 for (int c = 0; c < 4; ++c)

226 F[c] = FL[c] + sL * (Us[c] - UL[c]);

227 } else { // right star state

228 const Real fac = rR * (sR - unR) / (sR - sStar);

229 typename Model::State Us;

230 Us[0] = fac;

231 Us[in] = fac * sStar;

232 Us[it] = fac * (UR[it] / rR);

233 Us[3] = fac * (UR[3] / rR + (sStar - unR) * (sStar + pR / (rR * (sR - unR))));

234 for (int c = 0; c < 4; ++c)

235 F[c] = FR[c] + sR * (Us[c] - UR[c]);

236 }

237 return F;

238 } // end of canonical Euler 2D path (else of the if constexpr HasHLLCStructure)

239 }

240};

241

245inline constexpr Real kRoeEntropyFixFraction = Real(0.1);

246

257struct RoeFlux {

258 template <class Model>

259 POPS_HD typename Model::State operator()(const Model& m, const typename Model::State& UL,

260 const Aux& AL, const typename Model::State& UR,

261 const Aux& AR, int dir) const {

262 // CAPABILITY PATH (HasRoeDissipation): GENERIC Roe-like solver -- the dissipation

263 // d = |A_roe| (UR - UL) (linearization + eigenstructure + entropy fix) comes from the MODEL.

264 // A canonical Euler model WITHOUT a hook takes the historical branch below, bit-identical.

265 if constexpr (HasRoeDissipation<Model>) {

266 const auto FL = m.flux(UL, AL, dir);

267 const auto FR = m.flux(UR, AR, dir);

268 const typename Model::State d = m.roe_dissipation(UL, AL, UR, AR, dir);

269 typename Model::State F;

270 for (int c = 0; c < Model::n_vars; ++c)

271 F[c] = Real(0.5) * (FL[c] + FR[c]) - Real(0.5) * d[c];

272 return F;

273 } else {

274 const int in = (dir == 0) ? 1 : 2; // normal momentum

275 const int it = (dir == 0) ? 2 : 1; // tangential

276 const Real rL = UL[0], rR = UR[0];

277 const Real unL = UL[in] / rL, unR = UR[in] / rR;

278 const Real utL = UL[it] / rL, utR = UR[it] / rR;

279 const Real pL = m.pressure(UL), pR = m.pressure(UR);

280 const Real HL = (UL[3] + pL) / rL, HR = (UR[3] + pR) / rR;

281

282 // Roe average (weighted by sqrt(rho))

283 const Real sqL = std::sqrt(rL), sqR = std::sqrt(rR), den = sqL + sqR;

284 const Real un = (sqL * unL + sqR * unR) / den;

285 const Real ut = (sqL * utL + sqR * utR) / den;

286 const Real H = (sqL * HL + sqR * HR) / den;

287 const Real rho = sqL * sqR;

288 const Real q2 = un * un + ut * ut;

289 // gamma-1 derived from the ideal gas: p = (gamma-1) (E - 1/2 rho |v|^2)

290 const Real gm1 = pL / (UL[3] - Real(0.5) * rL * (unL * unL + utL * utL));

291 const Real c2 = gm1 * (H - Real(0.5) * q2);

292 const Real c = std::sqrt(c2);

293

294 // wave jumps and amplitudes

295 const Real dr = rR - rL, dp = pR - pL, dun = unR - unL, dut = utR - utL;

296 const Real a1 = (dp - rho * c * dun) / (Real(2) * c2); // un - c wave

297 const Real a2 = dr - dp / c2; // entropy, un

298 const Real a3 = rho * dut; // shear, un

299 const Real a5 = (dp + rho * c * dun) / (Real(2) * c2); // un + c wave

300

301 // |eigenvalue| with Harten entropy fix on the acoustic waves (1, 5).

302 // kRoeEntropyFixFraction = 0.1 is an EULER/ROE-SPECIFIC entropy policy (width of the parabolic

303 // smoothing as a fraction of the Roe sound speed, the usual value from the literature):

304 // it has no meaning for another hyperbolic system and must not be presented as a generic core

305 // parameter.

306 const Real eps = kRoeEntropyFixFraction * c;

307 auto absfix = [eps](Real l) {

308 const Real al = l < 0 ? -l : l;

309 return al < eps ? Real(0.5) * (l * l / eps + eps) : al;

310 };

311 const Real al1 = absfix(un - c), al2 = (un < 0 ? -un : un), al5 = absfix(un + c);

312

313 // dissipation Sum |lambda_k| a_k r_k, basis (rho, mom_n, mom_t, E)

314 const Real d_rho = al1 * a1 + al2 * a2 + al5 * a5;

315 const Real d_mn = al1 * a1 * (un - c) + al2 * a2 * un + al5 * a5 * (un + c);

316 const Real d_mt = al1 * a1 * ut + al2 * (a2 * ut + a3) + al5 * a5 * ut;

317 const Real d_E =

318 al1 * a1 * (H - un * c) + al2 * (a2 * Real(0.5) * q2 + a3 * ut) + al5 * a5 * (H + un * c);

319

320 const auto FL = m.flux(UL, AL, dir);

321 const auto FR = m.flux(UR, AR, dir);

322 typename Model::State F;

323 F[0] = Real(0.5) * (FL[0] + FR[0]) - Real(0.5) * d_rho;

324 F[in] = Real(0.5) * (FL[in] + FR[in]) - Real(0.5) * d_mn;

325 F[it] = Real(0.5) * (FL[it] + FR[it]) - Real(0.5) * d_mt;

326 F[3] = Real(0.5) * (FL[3] + FR[3]) - Real(0.5) * d_E;

327 return F;

328 } // end of canonical ideal-gas Euler 2D path (else of the if constexpr HasRoeDissipation)

329 }

330};

331

337using EulerHLLCFlux2D = HLLCFlux;

338using EulerRoeFlux2D = RoeFlux;

339

340} // namespace pops

pops::HasHLLCStructure

HLLC capability: the model provides the CONTACT wave speed and the STAR STATE on side k.

Definition numerical_flux.hpp:142

pops::HasRoeDissipation

Roe capability: the model provides its FULL Roe dissipation d = |A_roe(UL, UR)| (UR - UL) – Roe avera...

Definition numerical_flux.hpp:156

pops

Definition amr_hierarchy.hpp:29

pops::Real

double Real

Definition types.hpp:30

pops::hll_speeds

POPS_HD void hll_speeds(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir, Real &sL, Real &sR)

hll_speeds: Davis estimates for the signal speeds of the HLL/HLLC solvers.

Definition numerical_flux.hpp:78

pops::hll_flux_with_speeds

POPS_HD Model::State hll_flux_with_speeds(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir, Real sL, Real sR)

hll_flux_with_speeds: HLL flux from ALREADY estimated signal speeds (sL, sR).

Definition numerical_flux.hpp:96

pops::kRoeEntropyFixFraction

constexpr Real kRoeEntropyFixFraction

Width of the RoeFlux Harten entropy-fix smoothing, as a fraction of the Roe sound speed (eps = kRoeEn...

Definition numerical_flux.hpp:245

state.hpp

Pointwise types of the physics layer: StateVec<N> (conserved state) and Aux (auxiliary fields from th...

pops::Aux

POINTWISE auxiliary fields shared with the physics: single coupling channel.

Definition state.hpp:123

pops::HLLCFlux

HLLCFlux (HLL + Contact wave, Toro): 3 waves, resolves the contact discontinuity.

Definition numerical_flux.hpp:169

pops::HLLCFlux::operator()

POPS_HD Model::State operator()(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir) const

Definition numerical_flux.hpp:171

pops::HLLFlux

HLLFlux (Harten-Lax-van Leer): 2 signal speeds, less diffusive than Rusanov.

Definition numerical_flux.hpp:119

pops::HLLFlux::operator()

POPS_HD Model::State operator()(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir) const

Definition numerical_flux.hpp:121

pops::RoeFlux

RoeFlux: Roe linearization + Harten entropy fix (acoustic waves).

Definition numerical_flux.hpp:257

pops::RoeFlux::operator()

POPS_HD Model::State operator()(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir) const

Definition numerical_flux.hpp:259

pops::RusanovFlux

RusanovFlux (local Lax-Friedrichs): robust flux, compatible with any minimal PhysicalModel.

Definition numerical_flux.hpp:55

pops::RusanovFlux::operator()

POPS_HD Model::State operator()(const Model &m, const typename Model::State &UL, const Aux &AL, const typename Model::State &UR, const Aux &AR, int dir) const

Definition numerical_flux.hpp:57

types.hpp

Base scalar types and the POPS_HD macro (host+device portability).

POPS_HD

#define POPS_HD

Definition types.hpp:25

---

Generated by  1.9.8