Source code for pyro.compressible_sr.problems.rt2

"""A RT problem with two distinct modes: short wave length on the
left and long wavelength on the right.  This allows one to see
how the growth rate depends on wavenumber.
"""


import sys

import numpy as np

from pyro.compressible_sr import eos
from pyro.mesh import patch
from pyro.util import msg


[docs] def init_data(my_data, rp): """ initialize the rt problem """ msg.bold("initializing the rt problem...") # make sure that we are passed a valid patch object if not isinstance(my_data, patch.CellCenterData2d): print("ERROR: patch invalid in rt2.py") print(my_data.__class__) sys.exit() # 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") gamma = rp.get_param("eos.gamma") grav = rp.get_param("compressible.grav") dens1 = rp.get_param("rt2.dens1") dens2 = rp.get_param("rt2.dens2") p0 = rp.get_param("rt2.p0") amp = rp.get_param("rt2.amp") sigma = rp.get_param("rt2.sigma") # initialize the components, remember, that ener here is # rho*eint + 0.5*rho*v**2, where eint is the specific # internal energy (erg/g) xmom[:, :] = 0.0 ymom[:, :] = 0.0 dens[:, :] = 0.0 f_l = 18 f_r = 3 # set the density to be stratified in the y-direction myg = my_data.grid ycenter = 0.5*(myg.ymin + myg.ymax) p = myg.scratch_array() j = myg.jlo while j <= myg.jhi: if myg.y[j] < ycenter: dens[:, j] = dens1 p[:, j] = p0 + dens1*grav*myg.y[j] else: dens[:, j] = dens2 p[:, j] = p0 + dens1*grav*ycenter + dens2*grav*(myg.y[j] - ycenter) j += 1 idx_l = myg.x2d < (myg.xmax - myg.xmin)/3.0 idx_r = myg.x2d >= (myg.xmax - myg.xmin)/3.0 ymom[idx_l] = amp*np.sin(4.0*np.pi*f_l*myg.x2d[idx_l] / (myg.xmax-myg.xmin))*np.exp(-(myg.y2d[idx_l]-ycenter)**2/sigma**2) ymom[idx_r] = amp*np.sin(4.0*np.pi*f_r*myg.x2d[idx_r] / (myg.xmax-myg.xmin))*np.exp(-(myg.y2d[idx_r]-ycenter)**2/sigma**2) ymom *= dens rhoh = eos.rhoh_from_rho_p(gamma, dens, p) W = 1./np.sqrt(1-xmom**2-ymom**2) dens[:, :] *= W xmom[:, :] *= rhoh*W**2 ymom[:, :] *= rhoh*W**2 ener[:, :] = rhoh*W**2 - p - dens
[docs] def finalize(): """ print out any information to the user at the end of the run """