"""Implementation of the Colella 2nd order unsplit Godunov scheme. This
is a 2-dimensional implementation only. We assume that the grid is
uniform, but it is relatively straightforward to relax this
assumption.
There are several different options for this solver (they are all
discussed in the Colella paper).
* limiter: 0 = no limiting; 1 = 2nd order MC limiter; 2 = 4th order MC limiter
* riemann: HLLC or CGF (for Colella, Glaz, and Freguson solver)
* use_flattening: set to 1 to use the multidimensional flattening at shocks
* delta, z0, z1: flattening parameters (we use Colella 1990 defaults)
The grid indices look like::
j+3/2--+---------+---------+---------+
| | | |
j+1 _| | | |
| | | |
| | | |
j+1/2--+---------XXXXXXXXXXX---------+
| X X |
j _| X X |
| X X |
| X X |
j-1/2--+---------XXXXXXXXXXX---------+
| | | |
j-1 _| | | |
| | | |
| | | |
j-3/2--+---------+---------+---------+
| | | | | | |
i-1 i i+1
i-3/2 i-1/2 i+1/2 i+3/2
We wish to solve
.. math::
U_t + F^x_x + F^y_y = H
we want U_{i+1/2}^{n+1/2} -- the interface values that are input to
the Riemann problem through the faces for each zone.
Taylor expanding yields::
n+1/2 dU dU
U = U + 0.5 dx -- + 0.5 dt --
i+1/2,j,L i,j dx dt
dU dF^x dF^y
= U + 0.5 dx -- - 0.5 dt ( ---- + ---- - H )
i,j dx dx dy
dU dF^x dF^y
= U + 0.5 ( dx -- - dt ---- ) - 0.5 dt ---- + 0.5 dt H
i,j dx dx dy
dt dU dF^y
= U + 0.5 dx ( 1 - -- A^x ) -- - 0.5 dt ---- + 0.5 dt H
i,j dx dx dy
dt _ dF^y
= U + 0.5 ( 1 - -- A^x ) DU - 0.5 dt ---- + 0.5 dt H
i,j dx dy
+----------+-----------+ +----+----+ +---+---+
| | |
this is the monotonized this is the source term
central difference term transverse
flux term
There are two components, the central difference in the normal to the
interface, and the transverse flux difference. This is done for the
left and right sides of all 4 interfaces in a zone, which are then
used as input to the Riemann problem, yielding the 1/2 time interface
values::
n+1/2
U
i+1/2,j
Then, the zone average values are updated in the usual finite-volume
way::
n+1 n dt x n+1/2 x n+1/2
U = U + -- { F (U ) - F (U ) }
i,j i,j dx i-1/2,j i+1/2,j
dt y n+1/2 y n+1/2
+ -- { F (U ) - F (U ) }
dy i,j-1/2 i,j+1/2
Updating U_{i,j}:
* We want to find the state to the left and right (or top and bottom)
of each interface, ex. U_{i+1/2,j,[lr]}^{n+1/2}, and use them to
solve a Riemann problem across each of the four interfaces.
* U_{i+1/2,j,[lr]}^{n+1/2} is comprised of two parts, the computation
of the monotonized central differences in the normal direction
(eqs. 2.8, 2.10) and the computation of the transverse derivatives,
which requires the solution of a Riemann problem in the transverse
direction (eqs. 2.9, 2.14).
* the monotonized central difference part is computed using the
primitive variables.
* We compute the central difference part in both directions before
doing the transverse flux differencing, since for the high-order
transverse flux implementation, we use these as the input to the
transverse Riemann problem.
"""
import pyro.compressible_sr.interface as ifc
import pyro.mesh.array_indexer as ai
from pyro.compressible_sr import c2p
from pyro.mesh import reconstruction
from pyro.util import msg
[docs]
def unsplit_fluxes(my_data, my_aux, rp, ivars, solid, tc, dt):
"""
unsplitFluxes returns the fluxes through the x and y interfaces by
doing an unsplit reconstruction of the interface values and then
solving the Riemann problem through all the interfaces at once
currently we assume a gamma-law EOS
The runtime parameter grav is assumed to be the gravitational
acceleration in the y-direction
Parameters
----------
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
vars : Variables object
The Variables object that tells us which indices refer to which
variables
tc : TimerCollection object
The timers we are using to profile
dt : float
The timestep we are advancing through.
Returns
-------
out : ndarray, ndarray
The fluxes on the x- and y-interfaces
"""
tm_flux = tc.timer("unsplitFluxes")
tm_flux.begin()
myg = my_data.grid
gamma = rp.get_param("eos.gamma")
#=========================================================================
# compute the primitive variables
#=========================================================================
# Q = (rho, u, v, p, {X})
dens = my_data.get_var("density")
ymom = my_data.get_var("y-momentum")
q = cons_to_prim_wrapper(my_data.data, gamma, ivars, myg)
#=========================================================================
# compute the flattening coefficients
#=========================================================================
# there is a single flattening coefficient (xi) for all directions
use_flattening = rp.get_param("compressible.use_flattening")
if use_flattening:
msg.fail("ERROR: Flattening is not allowed for the sr compressible problem")
# monotonized central differences
tm_limit = tc.timer("limiting")
tm_limit.begin()
limiter = rp.get_param("compressible.limiter")
ldx = myg.scratch_array(nvar=ivars.nvar)
ldy = myg.scratch_array(nvar=ivars.nvar)
for n in range(ivars.nvar):
ldx[:, :, n] = reconstruction.limit(my_data.data[:, :, n], myg, 1, limiter)
ldy[:, :, n] = reconstruction.limit(my_data.data[:, :, n], myg, 2, limiter)
tm_limit.end()
#=========================================================================
# x-direction
#=========================================================================
# left and right primitive variable states
tm_states = tc.timer("interfaceStates")
tm_states.begin()
_U_l, _U_r = ifc.states(1, myg.ng, myg.dx, dt,
ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix,
ivars.naux, gamma, q, my_data.data, ldx)
U_xl = ai.ArrayIndexer(d=_U_l, grid=myg)
U_xr = ai.ArrayIndexer(d=_U_r, grid=myg)
tm_states.end()
#=========================================================================
# y-direction
#=========================================================================
# left and right primitive variable states
tm_states.begin()
_U_l, _U_r = ifc.states(2, myg.ng, myg.dy, dt,
ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix,
ivars.naux, gamma, q, my_data.data, ldy)
U_yl = ai.ArrayIndexer(d=_U_l, grid=myg)
U_yr = ai.ArrayIndexer(d=_U_r, grid=myg)
tm_states.end()
#=========================================================================
# apply source terms
#=========================================================================
grav = rp.get_param("compressible.grav")
ymom_src = my_aux.get_var("ymom_src")
ymom_src.v()[:, :] = dens.v()*grav
my_aux.fill_BC("ymom_src")
E_src = my_aux.get_var("E_src")
E_src.v()[:, :] = ymom.v()*grav
my_aux.fill_BC("E_src")
U_xl.v(buf=1, n=ivars.iymom)[:, :] += 0.5*dt*ymom_src.ip(-1, buf=1)
U_xl.v(buf=1, n=ivars.iener)[:, :] += 0.5*dt*E_src.ip(-1, buf=1)
U_xr.v(buf=1, n=ivars.iymom)[:, :] += 0.5*dt*ymom_src.v(buf=1)
U_xr.v(buf=1, n=ivars.iener)[:, :] += 0.5*dt*E_src.v(buf=1)
U_yl.v(buf=1, n=ivars.iymom)[:, :] += 0.5*dt*ymom_src.jp(-1, buf=1)
U_yl.v(buf=1, n=ivars.iener)[:, :] += 0.5*dt*E_src.jp(-1, buf=1)
U_yr.v(buf=1, n=ivars.iymom)[:, :] += 0.5*dt*ymom_src.v(buf=1)
U_yr.v(buf=1, n=ivars.iener)[:, :] += 0.5*dt*E_src.v(buf=1)
#=========================================================================
# compute transverse fluxes
#=========================================================================
tm_riem = tc.timer("riemann")
tm_riem.begin()
riemann = rp.get_param("compressible.riemann")
if riemann == "HLLC":
riemannFunc = ifc.riemann_hllc
elif riemann == "CGF":
riemannFunc = ifc.riemann_cgf
else:
msg.fail("ERROR: Riemann solver undefined")
q_xl = cons_to_prim_wrapper(U_xl, gamma, ivars, myg)
q_xr = cons_to_prim_wrapper(U_xr, gamma, ivars, myg)
q_yl = cons_to_prim_wrapper(U_yl, gamma, ivars, myg)
q_yr = cons_to_prim_wrapper(U_yr, gamma, ivars, myg)
# print("before transverse fluxes")
# print(f'rho = {q_xl.jp(0, n=ivars.irho, buf=2)}')
_fx = riemannFunc(1, myg.ng,
ivars.idens, ivars.ixmom,
ivars.iymom, ivars.iener, ivars.irhox,
ivars.irho, ivars.iu, ivars.iv, ivars.ip,
ivars.ix, ivars.naux,
solid.xl, solid.xr,
gamma, U_xl, U_xr, q_xl, q_xr)
_fy = riemannFunc(2, myg.ng,
ivars.idens, ivars.ixmom,
ivars.iymom, ivars.iener, ivars.irhox,
ivars.irho, ivars.iu, ivars.iv, ivars.ip,
ivars.ix, ivars.naux,
solid.yl, solid.yr,
gamma, U_yl, U_yr, q_yl, q_yr)
F_x = ai.ArrayIndexer(d=_fx, grid=myg)
F_y = ai.ArrayIndexer(d=_fy, grid=myg)
tm_riem.end()
#=========================================================================
# construct the interface values of U now
#=========================================================================
"""
finally, we can construct the state perpendicular to the interface
by adding the central difference part to the transverse flux
difference.
The states that we represent by indices i,j are shown below
(1,2,3,4):
j+3/2--+----------+----------+----------+
| | | |
| | | |
j+1 -+ | | |
| | | |
| | | | 1: U_xl[i,j, :] = U
j+1/2--+----------XXXXXXXXXXXX----------+ i-1/2,j,L
| X X |
| X X |
j -+ 1 X 2 X | 2: U_xr[i,j, :] = U
| X X | i-1/2,j,R
| X 4 X |
j-1/2--+----------XXXXXXXXXXXX----------+
| | 3 | | 3: U_yl[i,j, :] = U
| | | | i,j-1/2,L
j-1 -+ | | |
| | | |
| | | | 4: U_yr[i,j, :] = U
j-3/2--+----------+----------+----------+ i,j-1/2,R
| | | | | | |
i-1 i i+1
i-3/2 i-1/2 i+1/2 i+3/2
remember that the fluxes are stored on the left edge, so
F_x[i,j, :] = F_x
i-1/2, j
F_y[i,j, :] = F_y
i, j-1/2
"""
tm_transverse = tc.timer("transverse flux addition")
tm_transverse.begin()
dtdx = dt/myg.dx
dtdy = dt/myg.dy
b = (2, 1)
for n in range(ivars.nvar):
U_xl.v(buf=b, n=n)[:, :] += \
- 0.5*dtdy*(F_y.ip_jp(-1, 1, buf=b, n=n) - F_y.ip(-1, buf=b, n=n))
U_xr.v(buf=b, n=n)[:, :] += \
- 0.5*dtdy*(F_y.jp(1, buf=b, n=n) - F_y.v(buf=b, n=n))
U_yl.v(buf=b, n=n)[:, :] += \
- 0.5*dtdx*(F_x.ip_jp(1, -1, buf=b, n=n) - F_x.jp(-1, buf=b, n=n))
U_yr.v(buf=b, n=n)[:, :] += \
- 0.5*dtdx*(F_x.ip(1, buf=b, n=n) - F_x.v(buf=b, n=n))
tm_transverse.end()
#=========================================================================
# construct the fluxes normal to the interfaces
#=========================================================================
# up until now, F_x and F_y stored the transverse fluxes, now we
# overwrite with the fluxes normal to the interfaces
tm_riem.begin()
q_xl = cons_to_prim_wrapper(U_xl, gamma, ivars, myg)
q_xr = cons_to_prim_wrapper(U_xr, gamma, ivars, myg)
q_yl = cons_to_prim_wrapper(U_yl, gamma, ivars, myg)
q_yr = cons_to_prim_wrapper(U_yr, gamma, ivars, myg)
# print("before normal fluxes")
# print(f'F_rho = {F_x.jp(0, n=ivars.ixmom, buf=2)}')
_fx = riemannFunc(1, myg.ng,
ivars.idens, ivars.ixmom,
ivars.iymom, ivars.iener, ivars.irhox,
ivars.irho, ivars.iu, ivars.iv, ivars.ip,
ivars.ix, ivars.naux,
solid.xl, solid.xr,
gamma, U_xl, U_xr, q_xl, q_xr)
_fy = riemannFunc(2, myg.ng,
ivars.idens, ivars.ixmom,
ivars.iymom, ivars.iener, ivars.irhox,
ivars.irho, ivars.iu, ivars.iv, ivars.ip,
ivars.ix, ivars.naux,
solid.yl, solid.yr,
gamma, U_yl, U_yr, q_yl, q_yr)
F_x = ai.ArrayIndexer(d=_fx, grid=myg)
F_y = ai.ArrayIndexer(d=_fy, grid=myg)
tm_riem.end()
#=========================================================================
# apply artificial viscosity
#=========================================================================
cvisc = rp.get_param("compressible.cvisc")
_ax, _ay = ifc.artificial_viscosity(
myg.ng, myg.dx, myg.dy,
cvisc, q.v(n=ivars.iu, buf=myg.ng), q.v(n=ivars.iv, buf=myg.ng))
avisco_x = ai.ArrayIndexer(d=_ax, grid=myg)
avisco_y = ai.ArrayIndexer(d=_ay, grid=myg)
b = (2, 1)
for n in range(ivars.nvar):
# F_x = F_x + avisco_x * (U(i-1,j) - U(i,j))
var = my_data.get_var_by_index(n)
F_x.v(buf=b, n=n)[:, :] += \
avisco_x.v(buf=b)*(var.ip(-1, buf=b) - var.v(buf=b))
# F_y = F_y + avisco_y * (U(i,j-1) - U(i,j))
F_y.v(buf=b, n=n)[:, :] += \
avisco_y.v(buf=b)*(var.jp(-1, buf=b) - var.v(buf=b))
tm_flux.end()
return F_x, F_y
[docs]
def cons_to_prim_wrapper(U, gamma, ivars, myg):
"""
wrapper for fortran cons to prim routine
"""
q = myg.scratch_array(nvar=ivars.nq)
c2p.cons_to_prim(U, ivars.irho,
ivars.iu, ivars.iv,
ivars.ip, ivars.ix, ivars.irhox,
ivars.idens, ivars.ixmom, ivars.iymom,
ivars.iener, ivars.naux,
gamma, q)
return q