From 8a1e0078b18831c7a0366357e4adeb3419fa8c9a Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Thu, 4 May 2023 11:38:24 -0400 Subject: [PATCH 1/3] add low-Mach HLLC solver --- src/hydro/hydro.cpp | 8 ++ src/hydro/rsolvers/hydro_lhllc.hpp | 184 +++++++++++++++++++++++++++++ src/hydro/rsolvers/rsolvers.hpp | 1 + src/main.hpp | 2 +- 4 files changed, 194 insertions(+), 1 deletion(-) create mode 100644 src/hydro/rsolvers/hydro_lhllc.hpp diff --git a/src/hydro/hydro.cpp b/src/hydro/hydro.cpp index 48a1b17a..e4f07d5c 100644 --- a/src/hydro/hydro.cpp +++ b/src/hydro/hydro.cpp @@ -266,6 +266,8 @@ std::shared_ptr Initialize(ParameterInput *pin) { riemann = RiemannSolver::hlle; } else if (riemann_str == "hllc") { riemann = RiemannSolver::hllc; + } else if (riemann_str == "lhllc") { + riemann = RiemannSolver::lhllc; } else if (riemann_str == "hlld") { riemann = RiemannSolver::hlld; } else if (riemann_str == "none") { @@ -310,6 +312,12 @@ std::shared_ptr Initialize(ParameterInput *pin) { add_flux_fun(flux_functions); add_flux_fun(flux_functions); add_flux_fun(flux_functions); + add_flux_fun(flux_functions); + add_flux_fun(flux_functions); + add_flux_fun(flux_functions); + add_flux_fun(flux_functions); + add_flux_fun(flux_functions); + add_flux_fun(flux_functions); add_flux_fun(flux_functions); add_flux_fun(flux_functions); add_flux_fun(flux_functions); diff --git a/src/hydro/rsolvers/hydro_lhllc.hpp b/src/hydro/rsolvers/hydro_lhllc.hpp new file mode 100644 index 00000000..84edc913 --- /dev/null +++ b/src/hydro/rsolvers/hydro_lhllc.hpp @@ -0,0 +1,184 @@ +//======================================================================================== +// Athena++ astrophysical MHD code +// Copyright(C) 2014 James M. Stone and other code contributors +// Licensed under the 3-clause BSD License, see LICENSE file for details +//======================================================================================== +//! \file lhllc.cpp +//! \brief Low-mach HLLC Riemann solver for hydrodynamics, an extension of the HLLE fluxes to +//! include the contact wave. Only works for adiabatic hydrodynamics. +//! +//! REFERENCES: +//! - E.F. Toro, "Riemann Solvers and numerical methods for fluid dynamics", 2nd ed., +//! Springer-Verlag, Berlin, (1999) chpt. 10. +//! - P. Batten, N. Clarke, C. Lambert, and D. M. Causon, "On the Choice of Wavespeeds +//! for the HLLC Riemann Solver", SIAM J. Sci. & Stat. Comp. 18, 6, 1553-1570, (1997). +//! - T. Minoshima, T. Miyoshi, "A low-dissipation HLLD approximate Riemann solver for a +//! very wide range of Mach numbers", JCP 446 (2021) 110639. + +#ifndef RSOLVERS_HYDRO_LHLLC_HPP_ +#define RSOLVERS_HYDRO_LHLLC_HPP_ + +// C++ headers +#include // max(), min() +#include // sqrt() + +// AthenaPK headers +#include "../../main.hpp" +#include "rsolvers.hpp" + +//---------------------------------------------------------------------------------------- +//! \fn void Hydro::RiemannSolver +//! \brief The HLLC Riemann solver for adiabatic hydrodynamics (use HLLE for isothermal) + +template <> +struct Riemann { + static KOKKOS_INLINE_FUNCTION void + Solve(parthenon::team_mbr_t const &member, const int k, const int j, const int il, + const int iu, const int ivx, const ScratchPad2D &wl, + const ScratchPad2D &wr, VariableFluxPack &cons, + const AdiabaticHydroEOS &eos, const Real c_h) { + // dvn and dvt are stored in the last component of the scratch spaces + const auto nvar = cons.GetDim(4); + + int ivy = IV1 + ((ivx - IV1) + 1) % 3; + int ivz = IV1 + ((ivx - IV1) + 2) % 3; + Real gamma = eos.GetGamma(); + Real gm1 = gamma - 1.0; + Real igm1 = 1.0 / gm1; + + parthenon::par_for_inner(member, il, iu, [&](const int i) { + Real wli[(NHYDRO)], wri[(NHYDRO)]; + Real fl[(NHYDRO)], fr[(NHYDRO)], flxi[(NHYDRO)]; + //--- Step 1. Load L/R states into local variables + wli[IDN] = wl(IDN, i); + wli[IV1] = wl(ivx, i); + wli[IV2] = wl(ivy, i); + wli[IV3] = wl(ivz, i); + wli[IPR] = wl(IPR, i); + + wri[IDN] = wr(IDN, i); + wri[IV1] = wr(ivx, i); + wri[IV2] = wr(ivy, i); + wri[IV3] = wr(ivz, i); + wri[IPR] = wr(IPR, i); + + //--- Step 2. Compute middle state estimates with PVRS (Toro 10.5.2) + + Real al, ar, el, er; + Real cl = eos.SoundSpeed(wli); + Real cr = eos.SoundSpeed(wri); + const Real vsql = SQR(wli[IV1]) + SQR(wli[IV2]) + SQR(wli[IV3]); + const Real vsqr = SQR(wri[IV1]) + SQR(wri[IV2]) + SQR(wri[IV3]); + el = wli[IPR] * igm1 + 0.5 * wli[IDN] * vsql; + er = wri[IPR] * igm1 + 0.5 * wri[IDN] * vsqr; + Real rhoa = .5 * (wli[IDN] + wri[IDN]); // average density + Real ca = .5 * (cl + cr); // average sound speed + Real pmid = .5 * (wli[IPR] + wri[IPR] + (wli[IV1] - wri[IV1]) * rhoa * ca); + + //--- Step 3. Compute sound speed in L,R + + Real ql, qr; + ql = (pmid <= wli[IPR]) + ? 1.0 + : std::sqrt(1.0 + (gamma + 1) / (2 * gamma) * (pmid / wli[IPR] - 1.0)); + qr = (pmid <= wri[IPR]) + ? 1.0 + : std::sqrt(1.0 + (gamma + 1) / (2 * gamma) * (pmid / wri[IPR] - 1.0)); + + //--- Step 4. Compute the max/min wave speeds based on L/R + + al = wli[IV1] - cl * ql; + ar = wri[IV1] + cr * qr; + + const Real S_L = al; + const Real S_R = ar; + + Real bp = ar > 0.0 ? ar : (TINY_NUMBER); + Real bm = al < 0.0 ? al : -(TINY_NUMBER); + + //--- Step 5. Compute the contact wave speed and pressure + + Real vxl = wli[IV1] - al; + Real vxr = wri[IV1] - ar; + + Real tl = wli[IPR] + vxl * wli[IDN] * wli[IV1]; + Real tr = wri[IPR] + vxr * wri[IDN] * wri[IV1]; + + Real ml = wli[IDN] * vxl; + Real mr = -(wri[IDN] * vxr); + + // shock detector (Minoshima & Miyoshi) + Real cmax = std::max(cl, cr); + + // Determine the contact wave speed... + Real am = (tl - tr) / (ml + mr); + + // ...and the pressure at the contact surface (Minoshima & Miyoshi) + Real chi = std::min(1.0, std::sqrt(std::max(vsql, vsqr)) / cmax); + Real phi = chi * (2.0 - chi); + + const Real rho_L = wli[IDN]; + const Real rho_R = wri[IDN]; + const Real u_L = wli[IV1]; + const Real u_R = wri[IV1]; + const Real P_L = wli[IPR]; + const Real P_R = wri[IPR]; + + Real cp = (((S_R - u_R) * rho_R * P_L) - ((S_L - u_L) * rho_L * P_R) + + (phi * rho_L * rho_R * (S_R - u_R) * (S_L - u_L) * (u_R - u_L))) / + (((S_R - u_R) * rho_R) - ((S_L - u_L) * rho_L)); + + cp = cp > 0.0 ? cp : 0.0; + + //--- Step 6. Compute L/R fluxes along the line bm, bp + + vxl = wli[IV1] - bm; + vxr = wri[IV1] - bp; + + fl[IDN] = wli[IDN] * vxl; + fr[IDN] = wri[IDN] * vxr; + + fl[IV1] = wli[IDN] * wli[IV1] * vxl + wli[IPR]; + fr[IV1] = wri[IDN] * wri[IV1] * vxr + wri[IPR]; + + fl[IV2] = wli[IDN] * wli[IV2] * vxl; + fr[IV2] = wri[IDN] * wri[IV2] * vxr; + + fl[IV3] = wli[IDN] * wli[IV3] * vxl; + fr[IV3] = wri[IDN] * wri[IV3] * vxr; + + fl[IEN] = el * vxl + wli[IPR] * wli[IV1]; + fr[IEN] = er * vxr + wri[IPR] * wri[IV1]; + + //--- Step 8. Compute flux weights or scales + + Real sl, sr, sm; + if (am >= 0.0) { + sl = am / (am - bm); + sr = 0.0; + sm = -bm / (am - bm); + } else { + sl = 0.0; + sr = -am / (bp - am); + sm = bp / (bp - am); + } + + //--- Step 9. Compute the HLLC flux at interface, including weighted contribution + // of the flux along the contact + + flxi[IDN] = sl * fl[IDN] + sr * fr[IDN]; + flxi[IV1] = sl * fl[IV1] + sr * fr[IV1] + sm * cp; + flxi[IV2] = sl * fl[IV2] + sr * fr[IV2]; + flxi[IV3] = sl * fl[IV3] + sr * fr[IV3]; + flxi[IEN] = sl * fl[IEN] + sr * fr[IEN] + sm * cp * am; + + cons.flux(ivx, IDN, k, j, i) = flxi[IDN]; + cons.flux(ivx, ivx, k, j, i) = flxi[IV1]; + cons.flux(ivx, ivy, k, j, i) = flxi[IV2]; + cons.flux(ivx, ivz, k, j, i) = flxi[IV3]; + cons.flux(ivx, IEN, k, j, i) = flxi[IEN]; + }); + } +}; + +#endif // RSOLVERS_HYDRO_LHLLC_HPP_ diff --git a/src/hydro/rsolvers/rsolvers.hpp b/src/hydro/rsolvers/rsolvers.hpp index 91a00943..18b178e6 100644 --- a/src/hydro/rsolvers/rsolvers.hpp +++ b/src/hydro/rsolvers/rsolvers.hpp @@ -30,6 +30,7 @@ struct Riemann; #include "hydro_dc_llf.hpp" #include "hydro_hllc.hpp" #include "hydro_hlle.hpp" +#include "hydro_lhllc.hpp" // "none" solvers for runs/testing without fluid evolution, i.e., just reset fluxes template <> diff --git a/src/main.hpp b/src/main.hpp index b73efbf8..c83e737c 100644 --- a/src/main.hpp +++ b/src/main.hpp @@ -30,7 +30,7 @@ enum { // array indices for 1D primitives: velocity, transverse components of field enum { IV1 = 1, IV2 = 2, IV3 = 3, IPR = 4 }; -enum class RiemannSolver { undefined, none, hlle, llf, hllc, hlld }; +enum class RiemannSolver { undefined, none, hlle, llf, hllc, lhllc, hlld }; enum class Reconstruction { undefined, dc, plm, ppm, wenoz, weno3, limo3 }; enum class Integrator { undefined, rk1, rk2, vl2, rk3 }; enum class Fluid { undefined, euler, glmmhd }; From 4d9bb8ec473f75609276d42acfd6633a90cbe9ca Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Thu, 4 May 2023 11:42:22 -0400 Subject: [PATCH 2/3] remove comment about dvn, dvt --- src/hydro/rsolvers/hydro_lhllc.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/hydro/rsolvers/hydro_lhllc.hpp b/src/hydro/rsolvers/hydro_lhllc.hpp index 84edc913..3170b5f0 100644 --- a/src/hydro/rsolvers/hydro_lhllc.hpp +++ b/src/hydro/rsolvers/hydro_lhllc.hpp @@ -37,7 +37,6 @@ struct Riemann { const int iu, const int ivx, const ScratchPad2D &wl, const ScratchPad2D &wr, VariableFluxPack &cons, const AdiabaticHydroEOS &eos, const Real c_h) { - // dvn and dvt are stored in the last component of the scratch spaces const auto nvar = cons.GetDim(4); int ivy = IV1 + ((ivx - IV1) + 1) % 3; From 8b8946df75bf145ebacf2c896c484349f1dc99e0 Mon Sep 17 00:00:00 2001 From: Ben Wibking Date: Thu, 4 May 2023 11:44:45 -0400 Subject: [PATCH 3/3] fix formatting --- src/hydro/rsolvers/hydro_lhllc.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/hydro/rsolvers/hydro_lhllc.hpp b/src/hydro/rsolvers/hydro_lhllc.hpp index 3170b5f0..89b3a905 100644 --- a/src/hydro/rsolvers/hydro_lhllc.hpp +++ b/src/hydro/rsolvers/hydro_lhllc.hpp @@ -4,8 +4,8 @@ // Licensed under the 3-clause BSD License, see LICENSE file for details //======================================================================================== //! \file lhllc.cpp -//! \brief Low-mach HLLC Riemann solver for hydrodynamics, an extension of the HLLE fluxes to -//! include the contact wave. Only works for adiabatic hydrodynamics. +//! \brief Low-mach HLLC Riemann solver for hydrodynamics, an extension of the HLLE fluxes +//! to include the contact wave. Only works for adiabatic hydrodynamics. //! //! REFERENCES: //! - E.F. Toro, "Riemann Solvers and numerical methods for fluid dynamics", 2nd ed., @@ -124,8 +124,8 @@ struct Riemann { const Real P_R = wri[IPR]; Real cp = (((S_R - u_R) * rho_R * P_L) - ((S_L - u_L) * rho_L * P_R) + - (phi * rho_L * rho_R * (S_R - u_R) * (S_L - u_L) * (u_R - u_L))) / - (((S_R - u_R) * rho_R) - ((S_L - u_L) * rho_L)); + (phi * rho_L * rho_R * (S_R - u_R) * (S_L - u_L) * (u_R - u_L))) / + (((S_R - u_R) * rho_R) - ((S_L - u_L) * rho_L)); cp = cp > 0.0 ? cp : 0.0;