From 12750e4120fb4ea2dd0d85c0175d04fc61ab9b36 Mon Sep 17 00:00:00 2001 From: Zhi Date: Thu, 29 Aug 2024 15:39:51 -0400 Subject: [PATCH] change riemann interface --- pyro/compressible/riemann.py | 98 ++++++++--------------------- pyro/compressible/unsplit_fluxes.py | 18 ++++-- pyro/compressible_rk/fluxes.py | 9 ++- 3 files changed, 47 insertions(+), 78 deletions(-) diff --git a/pyro/compressible/riemann.py b/pyro/compressible/riemann.py index c265896fa..47d07beec 100644 --- a/pyro/compressible/riemann.py +++ b/pyro/compressible/riemann.py @@ -6,10 +6,10 @@ @njit(cache=True) -def riemann_cons(idir, ng, - idens, ixmom, iymom, iener, irhoX, nspec, - lower_solid, upper_solid, - gamma, U_l, U_r): +def riemann_cgf(idir, ng, + idens, ixmom, iymom, iener, irhoX, nspec, + lower_solid, upper_solid, + gamma, U_l, U_r): r""" Solve riemann shock tube problem for a general equation of state using the method of Colella, Glaz, and Ferguson. See @@ -310,53 +310,6 @@ def riemann_cons(idir, ng, return U_out -@njit(cache=True) -def riemann_cgf(idir, ng, coord_type, - idens, ixmom, iymom, iener, irhoX, nspec, - lower_solid, upper_solid, - gamma, U_l, U_r): - r""" - This uses riemann_cons to directly output fluxes instead of conserved states. - This is mainly used to be consistent with riemann_hllc. - - 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 - Fluxes - """ - - # get conserved states from U_l and U_r using riemann solver - U_state = riemann_cons(idir, ng, - idens, ixmom, iymom, iener, irhoX, nspec, - lower_solid, upper_solid, - gamma, U_l, U_r) - - # Construct the corresponding flux - F = consFlux(idir, coord_type, gamma, - idens, ixmom, iymom, iener, irhoX, nspec, - U_state) - - return F - - @njit(cache=True) def riemann_prim(idir, ng, irho, iu, iv, ip, iX, nspec, @@ -907,17 +860,17 @@ def riemann_hllc(idir, ng, coord_type, return F -def riemann_flux(U_xl, U_xr, U_yl, U_yr, - my_data, rp, ivars, solid, tc): +def riemann_flux(idir, U_l, U_r, my_data, rp, ivars, + lower_solid, upper_solid, tc): """ - This constructs the unsplit fluxes through the x and y interfaces - using the left and right conserved states by using the riemann solver. + This is the general interface that constructs the unsplit fluxes through + the x and y interfaces using the left and right conserved states by + using the riemann solver. Parameters ---------- - U_xl, U_xr, U_yl, U_yr: ndarray, ndarray, ndarray, ndarray - Conserved states in the left and right x-interface - and left and right y-interface. + U_l, U_r: ndarray, ndarray + Conserved states in the left and right interface my_data : CellCenterData2d object The data object containing the grid and advective scalar that we are advecting. @@ -926,8 +879,8 @@ def riemann_flux(U_xl, U_xr, U_yl, U_yr, ivars : Variables object The Variables object that tells us which indices refer to which variables - solid: A container class - This is used in Riemann solver to indicate which side has solid boundary + lower_solid, upper_solid : int + Are we at lower or upper solid boundaries? tc : TimerCollection object The timers we are using to profile @@ -953,22 +906,23 @@ def riemann_flux(U_xl, U_xr, U_yl, U_yr, else: msg.fail("ERROR: Riemann solver undefined") - _fx = riemannFunc(1, myg.ng, myg.coord_type, - ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, - solid.xl, solid.xr, - gamma, U_xl, U_xr) + _f = riemannFunc(idir, myg.ng, myg.coord_type, + ivars.idens, ivars.ixmom, ivars.iymom, + ivars.iener, ivars.irhox, ivars.naux, + lower_solid, upper_solid, + gamma, U_l, U_r) - _fy = riemannFunc(2, myg.ng, myg.coord_type, - ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux, - solid.yl, solid.yr, - gamma, U_yl, U_yr) - - F_x = ai.ArrayIndexer(d=_fx, grid=myg) - F_y = ai.ArrayIndexer(d=_fy, grid=myg) + # If riemann_method is not HLLC, then it outputs interface conserved states + if riemann_method != "HLLC": + _f = consFlux(idir, myg.coord_type, gamma, + ivars.idens, ivars.ixmom, ivars.iymom, + ivars.iener, ivars.irhox, ivars.naux, + _f) + F = ai.ArrayIndexer(d=_f, grid=myg) tm_riem.end() - return F_x, F_y + return F @njit(cache=True) diff --git a/pyro/compressible/unsplit_fluxes.py b/pyro/compressible/unsplit_fluxes.py index dce095dfe..9f9d8c452 100644 --- a/pyro/compressible/unsplit_fluxes.py +++ b/pyro/compressible/unsplit_fluxes.py @@ -284,8 +284,13 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): # ========================================================================= # Get the flux through x, y interface via riemann solver - F_x, F_y = riemann.riemann_flux(U_xl, U_xr, U_yl, U_yr, - my_data, rp, ivars, solid, tc) + F_x = riemann.riemann_flux(1, U_xl, U_xr, + my_data, rp, ivars, + solid.xl, solid.xr, tc) + + F_y = riemann.riemann_flux(2, U_yl, U_yr, + my_data, rp, ivars, + solid.yl, solid.yr, tc) # ========================================================================= # construct the interface values of U now @@ -368,8 +373,13 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt): # up until now, F_x and F_y stored the transverse fluxes, now we # overwrite with the fluxes normal to the interfaces - F_x, F_y = riemann.riemann_flux(U_xl, U_xr, U_yl, U_yr, - my_data, rp, ivars, solid, tc) + F_x = riemann.riemann_flux(1, U_xl, U_xr, + my_data, rp, ivars, + solid.xl, solid.xr, tc) + + F_y = riemann.riemann_flux(2, U_yl, U_yr, + my_data, rp, ivars, + solid.yl, solid.yr, tc) # ========================================================================= # apply artificial viscosity diff --git a/pyro/compressible_rk/fluxes.py b/pyro/compressible_rk/fluxes.py index 4ab49750b..cd0a665f3 100644 --- a/pyro/compressible_rk/fluxes.py +++ b/pyro/compressible_rk/fluxes.py @@ -160,8 +160,13 @@ def fluxes(my_data, rp, ivars, solid, tc): # ========================================================================= # construct the fluxes normal to the interfaces # ========================================================================= - F_x, F_y = riemann.riemann_flux(U_xl, U_xr, U_yl, U_yr, - my_data, rp, ivars, solid, tc) + F_x = riemann.riemann_flux(1, U_xl, U_xr, + my_data, rp, ivars, + solid.xl, solid.xr, tc) + + F_y = riemann.riemann_flux(2, U_yl, U_yr, + my_data, rp, ivars, + solid.yl, solid.yr, tc) # ========================================================================= # apply artificial viscosity