Skip to content

Commit

Permalink
remove old code
Browse files Browse the repository at this point in the history
  • Loading branch information
@zhichen3
zhichen3 committed Aug 1, 2024
1 parent 32cf030 commit 95538df
Show file tree
Hide file tree
Showing 2 changed files with 0 additions and 340 deletions.
3 changes: 0 additions & 3 deletions pyro/compressible/simulation.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -286,9 +286,6 @@ def evolve(self):
self.cc_data, self.rp,
self.ivars)

# Flux_x, Flux_y = flx.unsplit_fluxes(self.cc_data, self.aux_data, self.rp,
# self.ivars, self.solid, self.tc, self.dt)

old_dens = dens.copy()
old_xmom = xmom.copy()
old_ymom = ymom.copy()
Expand Down
337 changes: 0 additions & 337 deletions pyro/compressible/unsplit_fluxes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -572,340 +572,3 @@ def apply_artificial_viscosity(F_x, F_y, q,
avisco_y.v(buf=b)*(var.jp(-1, buf=b) - var.v(buf=b))

return F_x, F_y


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})

q = comp.cons_to_prim(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:
xi_x = reconstruction.flatten(myg, q, 1, ivars, rp)
xi_y = reconstruction.flatten(myg, q, 2, ivars, rp)

xi = reconstruction.flatten_multid(myg, q, xi_x, xi_y, ivars)
else:
xi = 1.0

# 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] = xi*reconstruction.limit(q[:, :, n], myg, 1, limiter)
ldy[:, :, n] = xi*reconstruction.limit(q[:, :, n], myg, 2, limiter)

tm_limit.end()

# =========================================================================
# x-direction
# =========================================================================

# left and right primitive variable states
tm_states = tc.timer("interfaceStates")
tm_states.begin()

V_l, V_r = ifc.states(1, myg.ng, myg.Lx, dt,
ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix,
ivars.naux,
gamma,
q, ldx)

tm_states.end()

# transform interface states back into conserved variables
U_xl = comp.prim_to_cons(V_l, gamma, ivars, myg)
U_xr = comp.prim_to_cons(V_r, gamma, ivars, myg)

# =========================================================================
# y-direction
# =========================================================================

# left and right primitive variable states
tm_states.begin()

_V_l, _V_r = ifc.states(2, myg.ng, myg.Ly, dt,
ivars.irho, ivars.iu, ivars.iv, ivars.ip, ivars.ix,
ivars.naux,
gamma,
q, ldy)
V_l = ai.ArrayIndexer(d=_V_l, grid=myg)
V_r = ai.ArrayIndexer(d=_V_r, grid=myg)

tm_states.end()

# transform interface states back into conserved variables
U_yl = comp.prim_to_cons(V_l, gamma, ivars, myg)
U_yr = comp.prim_to_cons(V_r, gamma, ivars, myg)

# =========================================================================
# apply source terms
# =========================================================================

dens = my_data.get_var("density")
xmom = my_data.get_var("x-momentum")
ymom = my_data.get_var("y-momentum")
pres = my_data.get_var("pressure")

dens_src = my_aux.get_var("dens_src")
xmom_src = my_aux.get_var("xmom_src")
ymom_src = my_aux.get_var("ymom_src")
E_src = my_aux.get_var("E_src")

grav = rp.get_param("compressible.grav")

# Calculate external source (gravity), geometric, and pressure terms
if myg.coord_type == 1:
# assume gravity points in r-direction in spherical.
dens_src.v()[:, :] = 0.0
xmom_src.v()[:, :] = dens.v()*grav + \
ymom.v()**2 / (dens.v()*myg.x2d.v()) - \
(pres.ip(1) - pres.v()) / myg.Lx.v()
ymom_src.v()[:, :] = -(pres.jp(1) - pres.v()) / myg.Ly.v() - \
dens.v()*xmom.v()*ymom.v() / (dens.v()*myg.x2d.v())
E_src.v()[:, :] = xmom.v()*grav

else:
# assume gravity points in y-direction in cartesian
dens_src.v()[:, :] = 0.0
xmom_src.v()[:, :] = 0.0
ymom_src.v()[:, :] = dens.v()*grav
E_src.v()[:, :] = ymom.v()*grav

my_aux.fill_BC("dens_src")
my_aux.fill_BC("xmom_src")
my_aux.fill_BC("ymom_src")
my_aux.fill_BC("E_src")

# ymom_xl[i,j] += 0.5*dt*dens[i-1,j]*grav
U_xl.v(buf=1, n=ivars.ixmom)[:, :] += 0.5*dt*xmom_src.ip(-1, buf=1)
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)

# ymom_xr[i,j] += 0.5*dt*dens[i,j]*grav
U_xr.v(buf=1, n=ivars.ixmom)[:, :] += 0.5*dt*xmom_src.v(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)

# ymom_yl[i,j] += 0.5*dt*dens[i,j-1]*grav
U_yl.v(buf=1, n=ivars.ixmom)[:, :] += 0.5*dt*xmom_src.jp(-1, 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)

# ymom_yr[i,j] += 0.5*dt*dens[i,j]*grav
U_yr.v(buf=1, n=ivars.ixmom)[:, :] += 0.5*dt*xmom_src.v(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")

riemannFunc = None
if riemann == "HLLC":
riemannFunc = ifc.riemann_hllc
elif riemann == "CGF":
riemannFunc = ifc.riemann_cgf
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)

_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)

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()

b = (2, 1)
hdtV = 0.5*dt / myg.V

for n in range(ivars.nvar):

# U_xl[i,j,:] = U_xl[i,j,:] - 0.5*dt/dy * (F_y[i-1,j+1,:] - F_y[i-1,j,:])
U_xl.v(buf=b, n=n)[:, :] += \
- hdtV.v(buf=b)*(F_y.ip_jp(-1, 1, buf=b, n=n)*myg.Ay_r.ip(-1, buf=b) -
F_y.ip(-1, buf=b, n=n)*myg.Ay_l.ip(-1, buf=b))

# U_xr[i,j,:] = U_xr[i,j,:] - 0.5*dt/dy * (F_y[i,j+1,:] - F_y[i,j,:])
U_xr.v(buf=b, n=n)[:, :] += \
- hdtV.v(buf=b)*(F_y.jp(1, buf=b, n=n)*myg.Ay_r.v(buf=b) -
F_y.v(buf=b, n=n)*myg.Ay_l.v(buf=b))

# U_yl[i,j,:] = U_yl[i,j,:] - 0.5*dt/dx * (F_x[i+1,j-1,:] - F_x[i,j-1,:])
U_yl.v(buf=b, n=n)[:, :] += \
- hdtV.v(buf=b)*(F_x.ip_jp(1, -1, buf=b, n=n)*myg.Ax_r.jp(-1, buf=b) -
F_x.jp(-1, buf=b, n=n)*myg.Ax_l.jp(-1, buf=b))

# U_yr[i,j,:] = U_yr[i,j,:] - 0.5*dt/dx * (F_x[i+1,j,:] - F_x[i,j,:])
U_yr.v(buf=b, n=n)[:, :] += \
- hdtV.v(buf=b)*(F_x.ip(1, buf=b, n=n)*myg.Ax_r.v(buf=b) -
F_x.v(buf=b, n=n)*myg.Ax_l.v(buf=b))

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()

_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)

_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)

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

0 comments on commit 95538df

Please sign in to comment.