"""A shock hitting a ramp at an oblique angle. This is based on
Woodward & Colella 1984."""
import math
import numpy as np
from pyro.util import msg
DEFAULT_INPUTS = "inputs.ramp"
# these defaults comes from the third test problem in
# Woodward and Colella. Journal of Computational Physics, 54, 115-174, 1984
#
# Also, the numbers are consistent with the ones in Castro
PROBLEM_PARAMS = {"ramp.rhol": 8.0, # post-shock initial density
"ramp.ul": 7.1447096, # post-shock initial x-velocity
"ramp.vl": -4.125, # post-shock initial y-velocity
"ramp.pl": 116.5, # post-shock initial pressure
"ramp.rhor": 1.4, # pre-shock initial density
"ramp.ur": 0.0, # pre-shock initial x-velocity
"ramp.vr": 0.0, # pre-shock initial y-velocity
"ramp.pr": 1.0} # pre-shock initial pressure
[docs]
def init_data(my_data, rp):
""" initialize the double Mach reflection problem """
if rp.get_param("driver.verbose"):
msg.bold("initializing the double Mach reflection problem...")
# get the density, momenta, and energy as separate variables
dens = my_data.get_var("density")
xmom = my_data.get_var("x-momentum")
ymom = my_data.get_var("y-momentum")
ener = my_data.get_var("energy")
# initialize the components, remember, that ener here is
# rho*eint + 0.5*rho*v**2, where eint is the specific
# internal energy (erg/g)
r_l = rp.get_param("ramp.rhol")
u_l = rp.get_param("ramp.ul")
v_l = rp.get_param("ramp.vl")
p_l = rp.get_param("ramp.pl")
r_r = rp.get_param("ramp.rhor")
u_r = rp.get_param("ramp.ur")
v_r = rp.get_param("ramp.vr")
p_r = rp.get_param("ramp.pr")
gamma = rp.get_param("eos.gamma")
energy_l = p_l/(gamma - 1.0) + 0.5*r_l*(u_l*u_l + v_l*v_l)
energy_r = p_r/(gamma - 1.0) + 0.5*r_r*(u_r*u_r + v_r*v_r)
# there is probably an easier way to do this, but for now, we
# will just do an explicit loop. Also, we really want to set
# the pressure and get the internal energy from that, and then
# compute the total energy (which is what we store). For now
# we will just fake this
myg = my_data.grid
dens[:, :] = 1.4
for j in range(myg.jlo, myg.jhi+1):
cy_up = myg.y[j] + 0.5*myg.dy*math.sqrt(3)
cy_down = myg.y[j] - 0.5*myg.dy*math.sqrt(3)
cy = np.array([cy_down, cy_up])
for i in range(myg.ilo, myg.ihi+1):
dens[i, j] = 0.0
xmom[i, j] = 0.0
ymom[i, j] = 0.0
ener[i, j] = 0.0
sf_up = math.tan(math.pi/3.0)*(myg.x[i] + 0.5*myg.dx*math.sqrt(3)-1.0/6.0)
sf_down = math.tan(math.pi/3.0)*(myg.x[i] - 0.5*myg.dx*math.sqrt(3)-1.0/6.0)
sf = np.array([sf_down, sf_up]) # initial shock front
for y in cy:
for shockfront in sf:
if y >= shockfront:
dens[i, j] = dens[i, j] + 0.25*r_l
xmom[i, j] = xmom[i, j] + 0.25*r_l*u_l
ymom[i, j] = ymom[i, j] + 0.25*r_l*v_l
ener[i, j] = ener[i, j] + 0.25*energy_l
else:
dens[i, j] = dens[i, j] + 0.25*r_r
xmom[i, j] = xmom[i, j] + 0.25*r_r*u_r
ymom[i, j] = ymom[i, j] + 0.25*r_r*v_r
ener[i, j] = ener[i, j] + 0.25*energy_r
[docs]
def finalize():
""" print out any information to the user at the end of the run """
