From 1d042dd72ec7643391c8b5e734cc99d31303edb7 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 19 Sep 2024 17:59:41 -0400 Subject: [PATCH 1/2] add an implementation of an HLLC solver with low Mach corrections this follows Minioshima & Miyoshi (2021) to get the correct asymptotic pressure behavior at low Mach numbers --- pyro/compressible/riemann.py | 195 ++++++++++++++++++++++++++++++++++- 1 file changed, 192 insertions(+), 3 deletions(-) diff --git a/pyro/compressible/riemann.py b/pyro/compressible/riemann.py index 834166ffe..789bfe024 100644 --- a/pyro/compressible/riemann.py +++ b/pyro/compressible/riemann.py @@ -863,6 +863,193 @@ 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)) + + # Estimate the star quantities -- we'll just do the PVRS + # estimate here. + + 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 * (un_l - un_r) * factor + + # 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)) + + # 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 +1094,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 +1112,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 +1124,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 From 78f7d98fbbabf2289e9171985faa0df8b6435e60 Mon Sep 17 00:00:00 2001 From: Michael Zingale Date: Thu, 19 Sep 2024 19:17:35 -0400 Subject: [PATCH 2/2] separate out the wave speeds --- pyro/compressible/riemann.py | 208 +++++++++++++++-------------------- 1 file changed, 89 insertions(+), 119 deletions(-) diff --git a/pyro/compressible/riemann.py b/pyro/compressible/riemann.py index 789bfe024..d567c6c32 100644 --- a/pyro/compressible/riemann.py +++ b/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 @@ -956,34 +953,8 @@ def riemann_hllc_lowspeed(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 -- we'll just do the PVRS - # estimate here. - - 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 * (un_l - un_r) * factor - - # 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 @@ -1022,7 +993,6 @@ def riemann_hllc_lowspeed(idir, ng, 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))