add an implementation of an HLLC solver with low Mach corrections (#282)

this follows Minioshima & Miyoshi (2021) to get the correct asymptotic
pressure behavior at low Mach numbers.  This was inspired by Quokka's
implementation.

Note: at the moment, this does not work with characteristic tracing, so
it should be used with `compressible_rk`.

Author: zingale

pyro/compressible/riemann.py

@@ -593,6 +593,91 @@ def riemann_prim(idir, ng,
     return q_int
 
 
+@njit(cache=True)
+def estimate_wave_speed(rho_l, u_l, p_l, c_l,
+                        rho_r, u_r, p_r, c_r,
+                        gamma):
+
+    # Estimate the star quantities -- use one of three methods to
+    # do this -- the primitive variable Riemann solver, the two
+    # shock approximation, or the two rarefaction approximation.
+    # Pick the method based on the pressure states at the
+    # interface.
+
+    p_max = max(p_l, p_r)
+    p_min = min(p_l, p_r)
+
+    Q = p_max / p_min
+
+    rho_avg = 0.5 * (rho_l + rho_r)
+    c_avg = 0.5 * (c_l + c_r)
+
+    # primitive variable Riemann solver (Toro, 9.3)
+    factor = rho_avg * c_avg
+
+    pstar = 0.5 * (p_l + p_r) + 0.5 * (u_l - u_r) * factor
+    ustar = 0.5 * (u_l + u_r) + 0.5 * (p_l - p_r) / factor
+
+    if Q > 2 and (pstar < p_min or pstar > p_max):
+
+        # use a more accurate Riemann solver for the estimate here
+
+        if pstar < p_min:
+
+            # 2-rarefaction Riemann solver
+            z = (gamma - 1.0) / (2.0 * gamma)
+            p_lr = (p_l / p_r)**z
+
+            ustar = (p_lr * u_l / c_l + u_r / c_r +
+                     2.0 * (p_lr - 1.0) / (gamma - 1.0)) / \
+                     (p_lr / c_l + 1.0 / c_r)
+
+            pstar = 0.5 * (p_l * (1.0 + (gamma - 1.0) * (u_l - ustar) /
+                                  (2.0 * c_l))**(1.0 / z) +
+                           p_r * (1.0 + (gamma - 1.0) * (ustar - u_r) /
+                                  (2.0 * c_r))**(1.0 / z))
+
+        else:
+
+            # 2-shock Riemann solver
+            A_r = 2.0 / ((gamma + 1.0) * rho_r)
+            B_r = p_r * (gamma - 1.0) / (gamma + 1.0)
+
+            A_l = 2.0 / ((gamma + 1.0) * rho_l)
+            B_l = p_l * (gamma - 1.0) / (gamma + 1.0)
+
+            # guess of the pressure
+            p_guess = max(0.0, pstar)
+
+            g_l = np.sqrt(A_l / (p_guess + B_l))
+            g_r = np.sqrt(A_r / (p_guess + B_r))
+
+            pstar = (g_l * p_l + g_r * p_r - (u_r - u_l)) / (g_l + g_r)
+
+            ustar = 0.5 * (u_l + u_r) + \
+                        0.5 * ((pstar - p_r) * g_r - (pstar - p_l) * g_l)
+
+    # estimate the nonlinear wave speeds
+
+    if pstar <= p_l:
+        # rarefaction
+        S_l = u_l - c_l
+    else:
+        # shock
+        S_l = u_l - c_l * np.sqrt(1.0 + ((gamma + 1.0) / (2.0 * gamma)) *
+                                  (pstar / p_l - 1.0))
+
+    if pstar <= p_r:
+        # rarefaction
+        S_r = u_r + c_r
+    else:
+        # shock
+        S_r = u_r + c_r * np.sqrt(1.0 + ((gamma + 1.0) / (2.0 / gamma)) *
+                                  (pstar / p_r - 1.0))
+
+    return S_l, S_r
+
+
 @njit(cache=True)
 def riemann_hllc(idir, ng,
                  idens, ixmom, iymom, iener, irhoX, nspec,
@@ -683,96 +768,8 @@ def riemann_hllc(idir, ng,
             c_l = max(smallc, np.sqrt(gamma * p_l / rho_l))
             c_r = max(smallc, np.sqrt(gamma * p_r / rho_r))
 
-            # Estimate the star quantities -- use one of three methods to
-            # do this -- the primitive variable Riemann solver, the two
-            # shock approximation, or the two rarefaction approximation.
-            # Pick the method based on the pressure states at the
-            # interface.
-
-            p_max = max(p_l, p_r)
-            p_min = min(p_l, p_r)
-
-            Q = p_max / p_min
-
-            rho_avg = 0.5 * (rho_l + rho_r)
-            c_avg = 0.5 * (c_l + c_r)
-
-            # primitive variable Riemann solver (Toro, 9.3)
-            factor = rho_avg * c_avg
-            # factor2 = rho_avg / c_avg
-
-            pstar = 0.5 * (p_l + p_r) + 0.5 * (un_l - un_r) * factor
-            ustar = 0.5 * (un_l + un_r) + 0.5 * (p_l - p_r) / factor
-
-            # rhostar_l = rho_l + (un_l - ustar) * factor2
-            # rhostar_r = rho_r + (ustar - un_r) * factor2
-
-            if Q > 2 and (pstar < p_min or pstar > p_max):
-
-                # use a more accurate Riemann solver for the estimate here
-
-                if pstar < p_min:
-
-                    # 2-rarefaction Riemann solver
-                    z = (gamma - 1.0) / (2.0 * gamma)
-                    p_lr = (p_l / p_r)**z
-
-                    ustar = (p_lr * un_l / c_l + un_r / c_r +
-                             2.0 * (p_lr - 1.0) / (gamma - 1.0)) / \
-                            (p_lr / c_l + 1.0 / c_r)
-
-                    pstar = 0.5 * (p_l * (1.0 + (gamma - 1.0) * (un_l - ustar) /
-                                          (2.0 * c_l))**(1.0 / z) +
-                                   p_r * (1.0 + (gamma - 1.0) * (ustar - un_r) /
-                                          (2.0 * c_r))**(1.0 / z))
-
-                    # rhostar_l = rho_l * (pstar / p_l)**(1.0 / gamma)
-                    # rhostar_r = rho_r * (pstar / p_r)**(1.0 / gamma)
-
-                else:
-
-                    # 2-shock Riemann solver
-                    A_r = 2.0 / ((gamma + 1.0) * rho_r)
-                    B_r = p_r * (gamma - 1.0) / (gamma + 1.0)
-
-                    A_l = 2.0 / ((gamma + 1.0) * rho_l)
-                    B_l = p_l * (gamma - 1.0) / (gamma + 1.0)
-
-                    # guess of the pressure
-                    p_guess = max(0.0, pstar)
-
-                    g_l = np.sqrt(A_l / (p_guess + B_l))
-                    g_r = np.sqrt(A_r / (p_guess + B_r))
-
-                    pstar = (g_l * p_l + g_r * p_r -
-                             (un_r - un_l)) / (g_l + g_r)
-
-                    ustar = 0.5 * (un_l + un_r) + \
-                        0.5 * ((pstar - p_r) * g_r - (pstar - p_l) * g_l)
-
-                    # rhostar_l = rho_l * (pstar / p_l + (gamma - 1.0) / (gamma + 1.0)) / \
-                    #     ((gamma - 1.0) / (gamma + 1.0) * (pstar / p_l) + 1.0)
-                    #
-                    # rhostar_r = rho_r * (pstar / p_r + (gamma - 1.0) / (gamma + 1.0)) / \
-                    #     ((gamma - 1.0) / (gamma + 1.0) * (pstar / p_r) + 1.0)
-
-            # estimate the nonlinear wave speeds
-
-            if pstar <= p_l:
-                # rarefaction
-                S_l = un_l - c_l
-            else:
-                # shock
-                S_l = un_l - c_l * np.sqrt(1.0 + ((gamma + 1.0) / (2.0 * gamma)) *
-                                           (pstar / p_l - 1.0))
-
-            if pstar <= p_r:
-                # rarefaction
-                S_r = un_r + c_r
-            else:
-                # shock
-                S_r = un_r + c_r * np.sqrt(1.0 + ((gamma + 1.0) / (2.0 / gamma)) *
-                                           (pstar / p_r - 1.0))
+            S_l, S_r = estimate_wave_speed(rho_l, un_l, p_l, c_l,
+                                           rho_r, un_r, p_r, c_r, gamma)
 
             #  We could just take S_c = u_star as the estimate for the
             #  contact speed, but we can actually do this more accurately
@@ -863,6 +860,166 @@ def riemann_hllc(idir, ng,
     return F
 
 
+@njit(cache=True)
+def riemann_hllc_lowspeed(idir, ng,
+                          idens, ixmom, iymom, iener, irhoX, nspec,
+                          lower_solid, upper_solid,  # pylint: disable=unused-argument
+                          gamma, U_l, U_r):
+    r"""
+    This is the HLLC Riemann solver based on Toro (2009) alternate formulation
+    (Eqs. 10.43 and 10.44) and the low Mach number asymptotic fix of
+    Minoshima & Miyoshi (2021).  It is also based on the Quokka implementation.
+
+    Parameters
+    ----------
+    idir : int
+        Are we predicting to the edges in the x-direction (1) or y-direction (2)?
+    ng : int
+        The number of ghost cells
+    nspec : int
+        The number of species
+    idens, ixmom, iymom, iener, irhoX : int
+        The indices of the density, x-momentum, y-momentum, internal energy density
+        and species partial densities in the conserved state vector.
+    lower_solid, upper_solid : int
+        Are we at lower or upper solid boundaries?
+    gamma : float
+        Adiabatic index
+    U_l, U_r : ndarray
+        Conserved state on the left and right cell edges.
+
+    Returns
+    -------
+    out : ndarray
+        Conserved flux
+    """
+
+    # Only Cartesian2d is supported in HLLC
+    coord_type = 0
+
+    qx, qy, nvar = U_l.shape
+
+    F = np.zeros((qx, qy, nvar))
+
+    smallc = 1.e-10
+    smallp = 1.e-10
+
+    nx = qx - 2 * ng
+    ny = qy - 2 * ng
+    ilo = ng
+    ihi = ng + nx
+    jlo = ng
+    jhi = ng + ny
+
+    for i in range(ilo - 1, ihi + 1):
+        for j in range(jlo - 1, jhi + 1):
+
+            D_star = np.zeros(nvar)
+
+            # primitive variable states
+            rho_l = U_l[i, j, idens]
+
+            # un = normal velocity; ut = transverse velocity
+            iun = -1
+            if idir == 1:
+                un_l = U_l[i, j, ixmom] / rho_l
+                ut_l = U_l[i, j, iymom] / rho_l
+                iun = ixmom
+            else:
+                un_l = U_l[i, j, iymom] / rho_l
+                ut_l = U_l[i, j, ixmom] / rho_l
+                iun = iymom
+
+            rhoe_l = U_l[i, j, iener] - 0.5 * rho_l * (un_l**2 + ut_l**2)
+
+            p_l = rhoe_l * (gamma - 1.0)
+            p_l = max(p_l, smallp)
+
+            rho_r = U_r[i, j, idens]
+
+            if idir == 1:
+                un_r = U_r[i, j, ixmom] / rho_r
+                ut_r = U_r[i, j, iymom] / rho_r
+            else:
+                un_r = U_r[i, j, iymom] / rho_r
+                ut_r = U_r[i, j, ixmom] / rho_r
+
+            rhoe_r = U_r[i, j, iener] - 0.5 * rho_r * (un_r**2 + ut_r**2)
+
+            p_r = rhoe_r * (gamma - 1.0)
+            p_r = max(p_r, smallp)
+
+            # compute the sound speeds
+            c_l = max(smallc, np.sqrt(gamma * p_l / rho_l))
+            c_r = max(smallc, np.sqrt(gamma * p_r / rho_r))
+
+            S_l, S_r = estimate_wave_speed(rho_l, un_l, p_l, c_l,
+                                           rho_r, un_r, p_r, c_r, gamma)
+
+            # We could just take S_c = u_star as the estimate for the
+            # contact speed, but we can actually do this more accurately
+            # by using the Rankine-Hugonoit jump conditions across each
+            # of the waves (see Toro 2009 Eq, 10.37, Batten et al. SIAM
+            # J. Sci. and Stat. Comp., 18:1553 (1997)
+
+            S_c = (p_r - p_l +
+                   rho_l * un_l * (S_l - un_l) -
+                   rho_r * un_r * (S_r - un_r)) / \
+                (rho_l * (S_l - un_l) - rho_r * (S_r - un_r))
+
+            # D* is used to control the pressure addition to the star flux
+            D_star[iun] = 1.0
+            D_star[iener] = S_c
+
+            # compute the fluxes corresponding to the left and right states
+
+            U_state_l = U_l[i, j, :]
+            F_l = consFlux(idir, coord_type, gamma,
+                           idens, ixmom, iymom, iener, irhoX, nspec,
+                           U_state_l)
+
+            U_state_r = U_r[i, j, :]
+            F_r = consFlux(idir, coord_type, gamma,
+                           idens, ixmom, iymom, iener, irhoX, nspec,
+                           U_state_r)
+
+            # compute the star pressure with the low Mach correction
+            # Minoshima & Miyoshi (2021) Eq. 23.  This is actually averaging
+            # the left and right pressures (see also Toro 2009 Eq. 10.42)
+
+            vmag_l = np.sqrt(un_l**2 + ut_l**2)
+            vmag_r = np.sqrt(un_r**2 + ut_r**2)
+
+            cs_max = max(c_l, c_r)
+            chi = min(1.0, max(vmag_l, vmag_r) / cs_max)
+            phi = chi * (2.0 - chi)
+            pstar_lr = 0.5 * (p_l + p_r) + \
+                0.5 * phi * (rho_l * (S_l - un_l) * (S_c - un_l) +
+                             rho_r * (S_r - un_r) * (S_c - un_r))
+
+            # figure out which region we are in and compute the state and
+            # the interface fluxes using the HLLC Riemann solver
+            if S_r <= 0.0:
+                # R region
+                F[i, j, :] = F_r[:]
+
+            elif S_c <= 0.0 < S_r:
+                # R* region
+                F[i, j, :] = (S_c * (S_r * U_state_r - F_r) + S_r * pstar_lr * D_star) / (S_r - S_c)
+
+            elif S_l < 0.0 < S_c:
+                # L* region
+                F[i, j, :] = (S_c * (S_l * U_state_l - F_l) + S_l * pstar_lr * D_star) / (S_l - S_c)
+
+            else:
+                # L region
+                F[i, j, :] = F_l[:]
+
+            # we should deal with solid boundaries somehow here
+
+    return F
+
+
 def riemann_flux(idir, U_l, U_r, my_data, rp, ivars,
                  lower_solid, upper_solid, tc, return_cons=False):
     """
@@ -907,7 +1064,9 @@ def riemann_flux(idir, U_l, U_r, my_data, rp, ivars,
     riemann_method = rp.get_param("compressible.riemann")
     gamma = rp.get_param("eos.gamma")
 
-    riemann_solvers = {"HLLC": riemann_hllc, "CGF": riemann_cgf}
+    riemann_solvers = {"HLLC": riemann_hllc,
+                       "HLLC_lm": riemann_hllc_lowspeed,
+                       "CGF": riemann_cgf}
 
     if riemann_method not in riemann_solvers:
         msg.fail("ERROR: Riemann solver undefined")
@@ -923,7 +1082,7 @@ def riemann_flux(idir, U_l, U_r, my_data, rp, ivars,
                      gamma, U_l, U_r)
 
     # If riemann_method is not HLLC, then construct flux using conserved states
-    if riemann_method != "HLLC":
+    if riemann_method not in ["HLLC", "HLLC_lm"]:
         _f = consFlux(idir, myg.coord_type, gamma,
                       ivars.idens, ivars.ixmom, ivars.iymom,
                       ivars.iener, ivars.irhox, ivars.naux,
@@ -935,7 +1094,7 @@ def riemann_flux(idir, U_l, U_r, my_data, rp, ivars,
     F = ai.ArrayIndexer(d=_f, grid=myg)
     tm_riem.end()
 
-    if riemann_method != "HLLC" and return_cons:
+    if riemann_method not in ["HLLC", "HLLC_lm"] and return_cons:
         U = ai.ArrayIndexer(d=_u, grid=myg)
         return F, U