Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

reorganize compressible riemann solver interface to prepare for SphericalPolar geometry #210

Merged
merged 7 commits into from
Aug 29, 2024
Merged

reorganize compressible riemann solver interface to prepare for SphericalPolar geometry #210

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
90850db
reorganize riemann to prepare for sphericalpolar geometry
zhichen3 Aug 29, 2024
2cb12f1
fix flake8
zhichen3 Aug 29, 2024
3d70863
fix flake8 and pylint
zhichen3 Aug 29, 2024
12750e4
change riemann interface
zhichen3 Aug 29, 2024
5f6e593
update riemann_flux
zhichen3 Aug 29, 2024
1942efa
Merge branch 'main' into compressible_riemann
zhichen3 Aug 29, 2024
00c58e9
fix
zhichen3 Aug 29, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
149 changes: 118 additions & 31 deletions pyro/compressible/riemann.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
import numpy as np
from numba import njit

import pyro.mesh.array_indexer as ai
from pyro.util import msg


@njit(cache=True)
def riemann_cgf(idir, ng,
Expand Down Expand Up @@ -56,12 +59,14 @@ def riemann_cgf(idir, ng,
Returns
-------
out : ndarray
Conserved flux
Conserved states.
"""

# pylint: disable=unused-variable

qx, qy, nvar = U_l.shape

F = np.zeros((qx, qy, nvar))
U_out = np.zeros((qx, qy, nvar))

smallc = 1.e-10
smallrho = 1.e-10
Expand Down Expand Up @@ -286,24 +291,23 @@ def riemann_cgf(idir, ng,
if j == jhi + 1 and upper_solid == 1:
un_state = 0.0

# compute the fluxes
F[i, j, idens] = rho_state * un_state
# update conserved state
U_out[i, j, idens] = rho_state

if idir == 1:
F[i, j, ixmom] = rho_state * un_state**2 + p_state
F[i, j, iymom] = rho_state * ut_state * un_state
U_out[i, j, ixmom] = rho_state * un_state
U_out[i, j, iymom] = rho_state * ut_state
else:
F[i, j, ixmom] = rho_state * ut_state * un_state
F[i, j, iymom] = rho_state * un_state**2 + p_state
U_out[i, j, ixmom] = rho_state * ut_state
U_out[i, j, iymom] = rho_state * un_state

F[i, j, iener] = rhoe_state * un_state + \
0.5 * rho_state * (un_state**2 + ut_state**2) * un_state + \
p_state * un_state
U_out[i, j, iener] = rhoe_state + \
0.5 * rho_state * (un_state**2 + ut_state**2)

if nspec > 0:
F[i, j, irhoX:irhoX + nspec] = xn * rho_state * un_state
U_out[i, j, irhoX:irhoX + nspec] = xn * rho_state

return F
return U_out


@njit(cache=True)
Expand Down Expand Up @@ -623,6 +627,9 @@ def riemann_hllc(idir, ng,
Conserved flux
"""

# Only Cartesian2d is supported in HLLC
coord_type = 0

qx, qy, nvar = U_l.shape

F = np.zeros((qx, qy, nvar))
Expand Down Expand Up @@ -781,7 +788,8 @@ def riemann_hllc(idir, ng,
# R region
U_state[:] = U_r[i, j, :]

F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec,
F[i, j, :] = consFlux(idir, coord_type, gamma,
idens, ixmom, iymom, iener, irhoX, nspec,
U_state)

elif S_c <= 0.0 < S_r:
Expand All @@ -806,7 +814,8 @@ def riemann_hllc(idir, ng,
U_r[i, j, irhoX:irhoX + nspec] / rho_r

# find the flux on the right interface
F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec,
F[i, j, :] = consFlux(idir, coord_type, gamma,
idens, ixmom, iymom, iener, irhoX, nspec,
U_r[i, j, :])

# correct the flux
Expand Down Expand Up @@ -834,7 +843,8 @@ def riemann_hllc(idir, ng,
U_l[i, j, irhoX:irhoX + nspec] / rho_l

# find the flux on the left interface
F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec,
F[i, j, :] = consFlux(idir, coord_type, gamma,
idens, ixmom, iymom, iener, irhoX, nspec,
U_l[i, j, :])

# correct the flux
Expand All @@ -844,16 +854,83 @@ def riemann_hllc(idir, ng,
# L region
U_state[:] = U_l[i, j, :]

F[i, j, :] = consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec,
F[i, j, :] = consFlux(idir, coord_type, gamma,
idens, ixmom, iymom, iener, irhoX, nspec,
U_state)

# 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):
"""
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_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.
rp : RuntimeParameters object
The runtime parameters for the simulation
ivars : Variables object
The Variables object that tells us which indices refer to which
variables
lower_solid, upper_solid : int
Are we at lower or upper solid boundaries?
tc : TimerCollection object
The timers we are using to profile

Returns
-------
out : ndarray, ndarray
Fluxes in x and y direction
"""

tm_riem = tc.timer("riemann")
tm_riem.begin()

myg = my_data.grid

riemann_method = rp.get_param("compressible.riemann")
gamma = rp.get_param("eos.gamma")

riemann_solvers = {"HLLC": riemann_hllc, "CGF": riemann_cgf}

if riemann_method not in riemann_solvers:
msg.fail("ERROR: Riemann solver undefined")

riemannFunc = riemann_solvers[riemann_method]

_f = riemannFunc(idir, myg.ng,
ivars.idens, ivars.ixmom, ivars.iymom,
ivars.iener, ivars.irhox, ivars.naux,
lower_solid, upper_solid,
gamma, U_l, U_r)

# 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


@njit(cache=True)
def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_state):
def consFlux(idir, coord_type, gamma,
idens, ixmom, iymom, iener, irhoX, nspec,
U_state):
r"""
Calculate the conservative flux.

Expand All @@ -879,27 +956,37 @@ def consFlux(idir, gamma, idens, ixmom, iymom, iener, irhoX, nspec, U_state):

F = np.zeros_like(U_state)

u = U_state[ixmom] / U_state[idens]
v = U_state[iymom] / U_state[idens]
u = U_state[..., ixmom] / U_state[..., idens]
v = U_state[..., iymom] / U_state[..., idens]

p = (U_state[iener] - 0.5 * U_state[idens] * (u * u + v * v)) * (gamma - 1.0)
p = (U_state[..., iener] - 0.5 * U_state[..., idens] * (u * u + v * v)) * (gamma - 1.0)

if idir == 1:
F[idens] = U_state[idens] * u
F[ixmom] = U_state[ixmom] * u + p
F[iymom] = U_state[iymom] * u
F[iener] = (U_state[iener] + p) * u
F[..., idens] = U_state[..., idens] * u
F[..., ixmom] = U_state[..., ixmom] * u

# if Cartesian2d, then add pressure to xmom flux
if coord_type == 0:
F[..., ixmom] += p

F[..., iymom] = U_state[..., iymom] * u
F[..., iener] = (U_state[..., iener] + p) * u

if nspec > 0:
F[irhoX:irhoX + nspec] = U_state[irhoX:irhoX + nspec] * u
F[..., irhoX:irhoX + nspec] = U_state[..., irhoX:irhoX + nspec] * u

else:
F[idens] = U_state[idens] * v
F[ixmom] = U_state[ixmom] * v
F[iymom] = U_state[iymom] * v + p
F[iener] = (U_state[iener] + p) * v
F[..., idens] = U_state[..., idens] * v
F[..., ixmom] = U_state[..., ixmom] * v
F[..., iymom] = U_state[..., iymom] * v

# if Cartesian2d, then add pressure to ymom flux
if coord_type == 0:
F[..., iymom] += p

F[..., iener] = (U_state[..., iener] + p) * v

if nspec > 0:
F[irhoX:irhoX + nspec] = U_state[irhoX:irhoX + nspec] * v
F[..., irhoX:irhoX + nspec] = U_state[..., irhoX:irhoX + nspec] * v

return F
55 changes: 13 additions & 42 deletions pyro/compressible/unsplit_fluxes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,6 @@
import pyro.mesh.array_indexer as ai
from pyro.compressible import riemann
from pyro.mesh import reconstruction
from pyro.util import msg


def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt):
Expand Down Expand Up @@ -283,33 +282,15 @@ def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt):
# =========================================================================
# compute transverse fluxes
# =========================================================================
tm_riem = tc.timer("riemann")
tm_riem.begin()

riemann_method = rp.get_param("compressible.riemann")
# Get the flux through x, y interface via riemann solver
F_x = riemann.riemann_flux(1, U_xl, U_xr,
my_data, rp, ivars,
solid.xl, solid.xr, tc)

riemannFunc = None
if riemann_method == "HLLC":
riemannFunc = riemann.riemann_hllc
elif riemann_method == "CGF":
riemannFunc = riemann.riemann_cgf
else:
msg.fail("ERROR: Riemann solver undefined")

_fx = riemannFunc(1, myg.ng,
ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux,
solid.xl, solid.xr,
gamma, U_xl, U_xr)

_fy = riemannFunc(2, myg.ng,
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)

tm_riem.end()
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
Expand Down Expand Up @@ -392,23 +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 = riemann.riemann_flux(1, U_xl, U_xr,
my_data, rp, ivars,
solid.xl, solid.xr, tc)

tm_riem.begin()

_fx = riemannFunc(1, myg.ng,
ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux,
solid.xl, solid.xr,
gamma, U_xl, U_xr)

_fy = riemannFunc(2, myg.ng,
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)

tm_riem.end()
F_y = riemann.riemann_flux(2, U_yl, U_yr,
my_data, rp, ivars,
solid.yl, solid.yr, tc)

# =========================================================================
# apply artificial viscosity
Expand Down
33 changes: 6 additions & 27 deletions pyro/compressible_rk/fluxes.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,7 +23,6 @@
import pyro.mesh.array_indexer as ai
from pyro.compressible import interface, riemann
from pyro.mesh import reconstruction
from pyro.util import msg


def fluxes(my_data, rp, ivars, solid, tc):
Expand Down Expand Up @@ -161,33 +160,13 @@ def fluxes(my_data, rp, ivars, solid, tc):
# =========================================================================
# construct the fluxes normal to the interfaces
# =========================================================================
tm_riem = tc.timer("Riemann")
tm_riem.begin()
F_x = riemann.riemann_flux(1, U_xl, U_xr,
my_data, rp, ivars,
solid.xl, solid.xr, tc)

riemann_method = rp.get_param("compressible.riemann")

riemannFunc = None
if riemann_method == "HLLC":
riemannFunc = riemann.riemann_hllc
elif riemann_method == "CGF":
riemannFunc = riemann.riemann_cgf
else:
msg.fail("ERROR: Riemann solver undefined")

_fx = riemannFunc(1, myg.ng,
ivars.idens, ivars.ixmom, ivars.iymom, ivars.iener, ivars.irhox, ivars.naux,
solid.xl, solid.xr,
gamma, U_xl, U_xr)

_fy = riemannFunc(2, myg.ng,
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)

tm_riem.end()
F_y = riemann.riemann_flux(2, U_yl, U_yr,
my_data, rp, ivars,
solid.yl, solid.yr, tc)

# =========================================================================
# apply artificial viscosity
Expand Down