def init_data(myd, rp):
    """ initialize the tophat advection problem """

    msg.bold("initializing the tophat advection problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(myd, patch.CellCenterData2d):
        print("ERROR: patch invalid in tophat.py")
        print(myd.__class__)
        sys.exit()

    dens = myd.get_var("density")

    xmin = myd.grid.xmin
    xmax = myd.grid.xmax

    ymin = myd.grid.ymin
    ymax = myd.grid.ymax

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)

    dens[:, :] = 0.0

    R = 0.1

    inside = (myd.grid.x2d - xctr)**2 + (myd.grid.y2d - yctr)**2 < R**2

    dens[inside] = 1.0

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 rt.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("rt.dens1")
    dens2 = rp.get_param("rt.dens2")
    p0 = rp.get_param("rt.p0")
    amp = rp.get_param("rt.amp")
    sigma = rp.get_param("rt.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

    # 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

    ymom[:, :] = amp*np.cos(2.0*np.pi*myg.x2d/(myg.xmax-myg.xmin))*np.exp(-(myg.y2d-ycenter)**2/sigma**2)

    ymom *= dens

    # set the energy (P = cs2*dens)
    ener[:, :] = p[:, :]/(gamma - 1.0) + \
        0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]

def init_data(my_data, rp):
    """ initialize the tophat advection problem """

    msg.bold("initializing the tophat advection problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print("ERROR: patch invalid in tophat.py")
        print(my_data.__class__)
        sys.exit()

    dens = my_data.get_var("density")

    xmin = my_data.grid.xmin
    xmax = my_data.grid.xmax

    ymin = my_data.grid.ymin
    ymax = my_data.grid.ymax

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)

    dens[:,:] = 0.0

    for i in range(my_data.grid.ilo, my_data.grid.ihi+1):
        for j in range(my_data.grid.jlo, my_data.grid.jhi+1):

            if (numpy.sqrt((my_data.grid.x[i]-xctr)**2 +
                           (my_data.grid.y[j]-yctr)**2) < 0.1):
                dens[i,j] = 1.0

    def __init__(self, solver_name):
        """
        Constructor

        Parameters
        ----------
        solver_name : str
            Name of solver to use
        """

        msg.bold('pyro ...')

        if solver_name not in valid_solvers:
            msg.fail("ERROR: %s is not a valid solver" % solver_name)

        self.pyro_home = os.path.dirname(os.path.realpath(__file__)) + '/'

        # import desired solver under "solver" namespace
        self.solver = importlib.import_module(solver_name)
        self.solver_name = solver_name

        # -------------------------------------------------------------------------
        # runtime parameters
        # -------------------------------------------------------------------------

        # parameter defaults
        self.rp = runparams.RuntimeParameters()
        self.rp.load_params(self.pyro_home + "_defaults")
        self.rp.load_params(self.pyro_home + solver_name + "/_defaults")

        self.tc = profile.TimerCollection()

        self.is_initialized = False

def init_data(my_data, rp):
    """ initialize the smooth advection problem """

    msg.bold("initializing the smooth advection problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print("ERROR: patch invalid in smooth.py")
        print(my_data.__class__)
        sys.exit()

    dens = my_data.get_var("density")

    xmin = my_data.grid.xmin
    xmax = my_data.grid.xmax

    ymin = my_data.grid.ymin
    ymax = my_data.grid.ymax

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)
    
    i = my_data.grid.ilo
    while i <= my_data.grid.ihi:

        j = my_data.grid.jlo
        while j <= my_data.grid.jhi:

            dens[i,j] = 1.0 + numpy.exp(-60.0*((my_data.grid.x[i]-xctr)**2 + \
                                               (my_data.grid.y[j]-yctr)**2))
                    
            j += 1
        i += 1

def initData(my_data):
    """ initialize the incompressible shear problem """

    msg.bold("initializing the incompressible shear problem...")

    rp = my_data.rp

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print my_data.__class__
        msg.fail("ERROR: patch invalid in shear.py")


    # get the necessary runtime parameters
    rho_s = rp.get_param("shear.rho_s")
    delta_s = rp.get_param("shear.delta_s")

    
    # get the velocities
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    myg = my_data.grid

    if (myg.xmin != 0 or myg.xmax != 1 or
        myg.ymin != 0 or myg.ymax != 1):
        msg.fail("ERROR: domain should be a unit square")
        
    y_half = 0.5*(myg.ymin + myg.ymax)

    print 'y_half = ', y_half
    print 'delta_s = ', delta_s
    print 'rho_s = ', rho_s
    
    # there is probably an easier way to do this without loops, but
    # for now, we will just do an explicit loop.
    i = myg.ilo
    while i <= myg.ihi:

        j = myg.jlo
        while j <= myg.jhi:

            if (myg.y[j] <= y_half):
                u[i,j] = numpy.tanh(rho_s*(myg.y[j] - 0.25))
            else:
                u[i,j] = numpy.tanh(rho_s*(0.75 - myg.y[j]))
            
            v[i,j] = delta_s*numpy.sin(2.0*math.pi*myg.x[i])
            
            j += 1
        i += 1
        
    
    print "extrema: ", numpy.min(u.flat), numpy.max(u.flat)

def init_data(my_data, rp):
    """ initialize the HSE problem """

    msg.bold("initializing the HSE problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print("ERROR: patch invalid in hse.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")

    dens0 = rp.get_param("hse.dens0")
    print("dens0 = ", dens0)
    H = rp.get_param("hse.h")

    # isothermal sound speed (squared)
    cs2 = H*abs(grav)

    # 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

    # set the density to be stratified in the y-direction
    myg = my_data.grid

    p = myg.scratch_array()

    for j in range(myg.jlo, myg.jhi+1):
        dens[:, j] = dens0*np.exp(-myg.y[j]/H)
        if j == myg.jlo:
            p[:, j] = dens[:, j]*cs2
        else:
            p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav

    # set the energy
    ener[:, :] = p[:, :]/(gamma - 1.0) + \
        0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]

def init_data(my_data, rp):
    """ initialize a smooth advection problem for testing convergence """

    msg.bold("initializing the advect problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print("ERROR: patch invalid in advect.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")

    # initialize the components, remember, that ener here is rho*eint
    # + 0.5*rho*v**2, where eint is the specific internal energy
    # (erg/g)
    dens.d[:,:] = 1.0
    xmom.d[:,:] = 0.0
    ymom.d[:,:] = 0.0


    gamma = rp.get_param("eos.gamma")

    xmin = rp.get_param("mesh.xmin")
    xmax = rp.get_param("mesh.xmax")

    ymin = rp.get_param("mesh.ymin")
    ymax = rp.get_param("mesh.ymax")

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)

    # this is identical to the advection/smooth problem
    dens.d[:,:] = 1.0 + np.exp(-60.0*((my_data.grid.x2d-xctr)**2 + 
                                      (my_data.grid.y2d-yctr)**2))


    # velocity is diagonal
    u = 1.0
    v = 1.0
    xmom.d[:,:] = dens.d[:,:]*u
    ymom.d[:,:] = dens.d[:,:]*v

    # pressure is constant
    p = 1.0
    ener.d[:,:] = p/(gamma - 1.0) + 0.5*(xmom.d[:,:]**2 + ymom.d[:,:]**2)/dens.d[:,:]

def init_data(my_data, rp):
    """ initialize the incompressible shear problem """

    msg.bold("initializing the incompressible shear problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print(my_data.__class__)
        msg.fail("ERROR: patch invalid in shear.py")

    # get the necessary runtime parameters
    eps = rp.get_param("vortex.eps")

    print('eps = ', eps)

    # get the velocities
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    myg = my_data.grid

    u.d[:,:] = -np.sin(math.pi*myg.y2d)
    v.d[:,:] = np.sin(math.pi*myg.x2d)
    #u.d[:,:] = -np.sin(2.0*math.pi*myg.x2d)*np.cos(2.0*math.pi*myg.y2d)*ran
    #v.d[:,:] = np.cos(2.0*math.pi*myg.x2d)*np.sin(2.0*math.pi*myg.y2d)*ran

    if eps != 0.0:
    #perturbed velocity1 at (0,0)
      r2 = myg.x2d**2+myg.y2d**2
      dvx1l = -eps**3*myg.y2d/r2*(1-np.exp(-r2/eps**2))
      dvy1l = eps**3*myg.x2d/r2*(1-np.exp(-r2/eps**2))

    #perturbed velocity1 at (2pi,0)
      r2 = (myg.x2d - 2.0)**2+myg.y2d**2
      dvx1r = -eps**3*myg.y2d/r2*(1-np.exp(-r2/eps**2))
      dvy1r = eps**3*(myg.x2d-2.0)/r2*(1-np.exp(-r2/eps**2))


    #perturbed velocity2 at (pi,0)
      r2 = (myg.x2d - 1.0)**2+myg.y2d**2
      dvx2 = eps**3*myg.y2d/r2*(1-np.exp(-r2/eps**2))
      dvy2 = -eps**3*(myg.x2d-1.0)/r2*(1-np.exp(-r2/eps**2))

      u.d[:,:] = u.d[:,:] + dvx1l + dvx1r + dvx2
      v.d[:,:] = v.d[:,:] + dvy1l + dvy1r + dvy2

    print("extrema: ", u.min(), u.max())

def init_data(myd, rp):
    """initialize the acoustic_pulse problem.  This comes from
    McCourquodale & Coella 2011"""

    msg.bold("initializing the acoustic pulse problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(myd, fv.FV2d):
        print("ERROR: patch invalid in acoustic_pulse.py")
        print(myd.__class__)
        sys.exit()

    # get the density, momenta, and energy as separate variables
    dens = myd.get_var("density")
    xmom = myd.get_var("x-momentum")
    ymom = myd.get_var("y-momentum")
    ener = myd.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)
    xmom[:, :] = 0.0
    ymom[:, :] = 0.0

    gamma = rp.get_param("eos.gamma")

    rho0 = rp.get_param("acoustic_pulse.rho0")
    drho0 = rp.get_param("acoustic_pulse.drho0")

    xmin = rp.get_param("mesh.xmin")
    xmax = rp.get_param("mesh.xmax")

    ymin = rp.get_param("mesh.ymin")
    ymax = rp.get_param("mesh.ymax")

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)

    dist = np.sqrt((myd.grid.x2d - xctr)**2 +
                   (myd.grid.y2d - yctr)**2)

    dens[:, :] = rho0
    idx = dist <= 0.5
    dens[idx] = rho0 + drho0*np.exp(-16*dist[idx]**2) * np.cos(np.pi*dist[idx])**6

    p = (dens/rho0)**gamma
    ener[:, :] = p/(gamma - 1)

def init_data(my_data, rp):
    """ initialize the slotted advection problem """
    msg.bold("initializing the slotted advection problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print(my_data.__class__)
        msg.fail("ERROR: patch invalid in slotted.py")

    offset = rp.get_param("slotted.offset")
    omega = rp.get_param("slotted.omega")

    myg = my_data.grid

    xctr_dens = 0.5*(myg.xmin + myg.xmax)
    yctr_dens = 0.5*(myg.ymin + myg.ymax) + offset

    # setting initial condition for density
    dens = my_data.get_var("density")
    dens[:, :] = 0.0

    R = 0.15
    slot_width = 0.05

    inside = (myg.x2d - xctr_dens)**2 + (myg.y2d - yctr_dens)**2 < R**2

    slot_x = np.logical_and(myg.x2d > (xctr_dens - slot_width*0.5),
                            myg.x2d < (xctr_dens + slot_width*0.5))
    slot_y = np.logical_and(myg.y2d > (yctr_dens - R),
                            myg.y2d < (yctr_dens))
    slot = np.logical_and(slot_x, slot_y)

    dens[inside] = 1.0
    dens[slot] = 0.0

    # setting initial condition for velocity
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    u[:, :] = omega*(myg.y2d - xctr_dens)
    v[:, :] = -omega*(myg.x2d - (yctr_dens-offset))

    print("extrema: ", np.amax(u), np.amin(u))

def init_data(my_data, rp):
    """ initialize the incompressible shear problem """

    msg.bold("initializing the incompressible shear problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print(my_data.__class__)
        msg.fail("ERROR: patch invalid in shear.py")


    # get the necessary runtime parameters
    rho_s = rp.get_param("shear.rho_s")
    delta_s = rp.get_param("shear.delta_s")


    # get the velocities
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    myg = my_data.grid

    if (myg.xmin != 0 or myg.xmax != 1 or
        myg.ymin != 0 or myg.ymax != 1):
        msg.fail("ERROR: domain should be a unit square")

    y_half = 0.5*(myg.ymin + myg.ymax)

    print('y_half = ', y_half)
    print('delta_s = ', delta_s)
    print('rho_s = ', rho_s)

    idx = myg.y2d <= y_half
    u.d[idx] = np.tanh(rho_s*(myg.y2d[idx] - 0.25))

    idx = myg.y2d > y_half
    u.d[idx] = np.tanh(rho_s*(0.75 - myg.y2d[idx]))

    v.d[:,:] = delta_s*np.sin(2.0*math.pi*myg.x2d)

    print("extrema: ", u.min(), u.max())

def init_data(my_data, rp):
    """ initialize the incompressible converge problem """

    msg.bold("initializing the incompressible converge problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print(my_data.__class__)
        msg.fail("ERROR: patch invalid in converge.py")

    # get the velocities
    u = my_data.get_var("x-velocity")
    v = my_data.get_var("y-velocity")

    myg = my_data.grid

    if (myg.xmin != 0 or myg.xmax != 1 or
        myg.ymin != 0 or myg.ymax != 1):
        msg.fail("ERROR: domain should be a unit square")

    u[:, :] = 1.0 - 2.0*np.cos(2.0*math.pi*myg.x2d)*np.sin(2.0*math.pi*myg.y2d)
    v[:, :] = 1.0 + 2.0*np.sin(2.0*math.pi*myg.x2d)*np.cos(2.0*math.pi*myg.y2d)

def initData(my_data):
    """ initialize the Gaussian diffusion problem """

    msg.bold("initializing the Gaussian diffusion problem...")

    rp = my_data.rp

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print "ERROR: patch invalid in diffuse.py"
        print my_data.__class__
        sys.exit()

    phi = my_data.get_var("phi")

    xmin = my_data.grid.xmin
    xmax = my_data.grid.xmax

    ymin = my_data.grid.ymin
    ymax = my_data.grid.ymax

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)
    
    k = rp.get_param("diffusion.k")
    t_0 = rp.get_param("gaussian.t_0")
    phi_max = rp.get_param("gaussian.phi_max")
    phi_0   = rp.get_param("gaussian.phi_0")

    dist = numpy.sqrt((my_data.grid.x2d - xctr)**2 +
                      (my_data.grid.y2d - yctr)**2)
    
    phi[:,:] = phi_analytic(dist, 0.0, t_0, k, phi_0, phi_max)

    # for later interpretation / analysis, store some auxillary data
    my_data.set_aux("k", k)
    my_data.set_aux("t_0", t_0)
    my_data.set_aux("phi_0", phi_0)
    my_data.set_aux("phi_max", phi_max)

def init_data(my_data, rp):
    """ initialize the smooth advection problem """

    msg.bold("initializing the smooth FV advection problem...")

    # make sure that we are passed a valid patch object
    # if not isinstance(my_data, patch.FV2d):
    #    print("ERROR: patch invalid in smooth.py")
    #    print(my_data.__class__)
    #    sys.exit()

    xmin = my_data.grid.xmin
    xmax = my_data.grid.xmax

    ymin = my_data.grid.ymin
    ymax = my_data.grid.ymax

    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)

    # we need to initialize the cell-averages, so we will create
    # a finer grid, initialize it, and then average down
    mgf = my_data.grid.fine_like(4)

    # since restrict operates in the data class, we need to
    # create a FV2d object here
    fine_data = fv.FV2d(mgf)
    fine_data.register_var("density", my_data.BCs["density"])
    fine_data.create()

    dens_fine = fine_data.get_var("density")

    dens_fine[:, :] = 1.0 + numpy.exp(-60.0*((mgf.x2d-xctr)**2 +
                                             (mgf.y2d-yctr)**2))

    dens = my_data.get_var("density")
    dens[:, :] = fine_data.restrict("density", N=4)

def initData(my_data):
    """ initialize the sod problem """

    msg.bold("initializing the sod problem...")

    rp = my_data.rp

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print "ERROR: patch invalid in sod.py"
        print my_data.__class__
        sys.exit()


    # get the sod parameters
    dens_left = rp.get_param("sod.dens_left")
    dens_right = rp.get_param("sod.dens_right")

    u_left = rp.get_param("sod.u_left")
    u_right = rp.get_param("sod.u_right")

    p_left = rp.get_param("sod.p_left")
    p_right = rp.get_param("sod.p_right")
    

    # 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)
    xmin = rp.get_param("mesh.xmin")
    xmax = rp.get_param("mesh.xmax")

    ymin = rp.get_param("mesh.ymin")
    ymax = rp.get_param("mesh.ymax")

    gamma = rp.get_param("eos.gamma")

    direction = rp.get_param("sod.direction")
    
    xctr = 0.5*(xmin + xmax)
    yctr = 0.5*(ymin + ymax)

    myg = my_data.grid
    
    # 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 pressue 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.
    if direction == "x":

        i = myg.ilo
        while i <= myg.ihi:

            j = myg.jlo
            while j <= myg.jhi:

                if myg.x[i] <= xctr:
                    dens[i,j] = dens_left
                    xmom[i,j] = dens_left*u_left
                    ymom[i,j] = 0.0
                    ener[i,j] = p_left/(gamma - 1.0) + 0.5*xmom[i,j]*u_left
                
                else:
                    dens[i,j] = dens_right
                    xmom[i,j] = dens_right*u_right
                    ymom[i,j] = 0.0
                    ener[i,j] = p_right/(gamma - 1.0) + 0.5*xmom[i,j]*u_right
                    
                j += 1
            i += 1

    else:
        i = myg.ilo
        while i <= myg.ihi:

            j = myg.jlo
            while j <= myg.jhi:

                if myg.y[j] <= yctr:
                    dens[i,j] = dens_left
                    xmom[i,j] = 0.0
                    ymom[i,j] = dens_left*u_left
                    ener[i,j] = p_left/(gamma - 1.0) + 0.5*ymom[i,j]*u_left
                
                else:
                    dens[i,j] = dens_right
                    xmom[i,j] = 0.0
                    ymom[i,j] = dens_right*u_right
                    ener[i,j] = p_right/(gamma - 1.0) + 0.5*ymom[i,j]*u_right
                    
                j += 1
            i += 1

def doit(solver_name, problem_name, param_file,
         other_commands=None,
         comp_bench=False, make_bench=False):

    msg.bold('pyro ...')

    tc = profile.TimerCollection()

    tm_main = tc.timer("main")
    tm_main.begin()

    # import desired solver under "solver" namespace
    solver = importlib.import_module(solver_name)

    #-------------------------------------------------------------------------
    # runtime parameters
    #-------------------------------------------------------------------------

    # parameter defaults
    rp = runparams.RuntimeParameters()
    rp.load_params("_defaults")
    rp.load_params(solver_name + "/_defaults")

    # problem-specific runtime parameters
    rp.load_params(solver_name + "/problems/_" + problem_name + ".defaults")

    # now read in the inputs file
    if not os.path.isfile(param_file):
        # check if the param file lives in the solver's problems directory
        param_file = solver_name + "/problems/" + param_file
        if not os.path.isfile(param_file):
            msg.fail("ERROR: inputs file does not exist")

    rp.load_params(param_file, no_new=1)

    # and any commandline overrides
    if not other_commands == None:
        rp.command_line_params(other_commands)

    # write out the inputs.auto
    rp.print_paramfile()


    #-------------------------------------------------------------------------
    # initialization
    #-------------------------------------------------------------------------

    # initialize the Simulation object -- this will hold the grid and
    # data and know about the runtime parameters and which problem we
    # are running
    sim = solver.Simulation(solver_name, problem_name, rp, timers=tc)

    sim.initialize()
    sim.preevolve()


    #-------------------------------------------------------------------------
    # evolve
    #-------------------------------------------------------------------------
    init_tstep_factor = rp.get_param("driver.init_tstep_factor")
    max_dt_change = rp.get_param("driver.max_dt_change")
    fix_dt = rp.get_param("driver.fix_dt")

    verbose = rp.get_param("driver.verbose")

    plt.ion()

    sim.cc_data.t = 0.0

    # output the 0th data
    basename = rp.get_param("io.basename")
    sim.cc_data.write("{}{:04d}".format(basename, sim.n))

    dovis = rp.get_param("vis.dovis")
    if dovis:
        plt.figure(num=1, figsize=(8,6), dpi=100, facecolor='w')
        sim.dovis()

    while not sim.finished():

        # fill boundary conditions
        sim.cc_data.fill_BC_all()

        # get the timestep
        if fix_dt > 0.0:
            sim.dt = fix_dt
        else:
            sim.compute_timestep()
            if sim.n == 0:
                sim.dt = init_tstep_factor*sim.dt
            else:
                sim.dt = min(max_dt_change*dt_old, sim.dt)
            dt_old = sim.dt

        if sim.cc_data.t + sim.dt > sim.tmax:
            sim.dt = sim.tmax - sim.cc_data.t

        # evolve for a single timestep
        sim.evolve()

#.........这里部分代码省略.........

def init_data(my_data, rp):
    """ initialize the quadrant problem """

    msg.bold("initializing the quadrant problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print("ERROR: patch invalid in quad.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")

    # initialize the components, remember, that ener here is
    # rho*eint + 0.5*rho*v**2, where eint is the specific
    # internal energy (erg/g)
    r1 = rp.get_param("quadrant.rho1")
    u1 = rp.get_param("quadrant.u1")
    v1 = rp.get_param("quadrant.v1")
    p1 = rp.get_param("quadrant.p1")

    r2 = rp.get_param("quadrant.rho2")
    u2 = rp.get_param("quadrant.u2")
    v2 = rp.get_param("quadrant.v2")
    p2 = rp.get_param("quadrant.p2")

    r3 = rp.get_param("quadrant.rho3")
    u3 = rp.get_param("quadrant.u3")
    v3 = rp.get_param("quadrant.v3")
    p3 = rp.get_param("quadrant.p3")

    r4 = rp.get_param("quadrant.rho4")
    u4 = rp.get_param("quadrant.u4")
    v4 = rp.get_param("quadrant.v4")
    p4 = rp.get_param("quadrant.p4")

    cx = rp.get_param("quadrant.cx")
    cy = rp.get_param("quadrant.cy")
    
    gamma = rp.get_param("eos.gamma")
    
    # 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 pressue 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

    iq1 = np.logical_and(myg.x2d >= cx, myg.y2d >= cy)
    iq2 = np.logical_and(myg.x2d < cx,  myg.y2d >= cy)
    iq3 = np.logical_and(myg.x2d < cx,  myg.y2d < cy)
    iq4 = np.logical_and(myg.x2d >= cx, myg.y2d < cy)    
    
    # quadrant 1
    dens.d[iq1] = r1
    xmom.d[iq1] = r1*u1
    ymom.d[iq1] = r1*v1
    ener.d[iq1] = p1/(gamma - 1.0) + 0.5*r1*(u1*u1 + v1*v1)
                
    # quadrant 2
    dens.d[iq2] = r2
    xmom.d[iq2] = r2*u2
    ymom.d[iq2] = r2*v2
    ener.d[iq2] = p2/(gamma - 1.0) + 0.5*r2*(u2*u2 + v2*v2)

    # quadrant 3
    dens.d[iq3] = r3
    xmom.d[iq3] = r3*u3
    ymom.d[iq3] = r3*v3
    ener.d[iq3] = p3/(gamma - 1.0) + 0.5*r3*(u3*u3 + v3*v3)

    # quadrant 4
    dens.d[iq4] = r4
    xmom.d[iq4] = r4*u4
    ymom.d[iq4] = r4*v4
    ener.d[iq4] = p4/(gamma - 1.0) + 0.5*r4*(u4*u4 + v4*v4)

                               to the stored benchmark for this problem
                               (looking in the solver's tests/ sub-
                               directory).

      [runtime parameters] override any of the runtime defaults of
      parameters specified in the inputs file.  For instance, to turn
      off runtime visualization, add:

         vis.dovis=0

      to the end of the commandline.

      """


msg.bold('pyro ...')

tc = profile.TimerCollection()

tm_main = tc.timer("main")
tm_main.begin()

#-----------------------------------------------------------------------------
# command line arguments / solver setup
#-----------------------------------------------------------------------------

# parse the runtime arguments.  We specify a solver (which we import
# locally under the namespace 'solver', the problem name, and the
# input file name

if len(sys.argv) == 1:

def init_data(my_data, rp):
    """ initialize the bubble problem """

    msg.bold("initializing the bubble problem...")

    # make sure that we are passed a valid patch object
    if not isinstance(my_data, patch.CellCenterData2d):
        print("ERROR: patch invalid in bubble.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")

    scale_height = rp.get_param("bubble.scale_height")
    dens_base = rp.get_param("bubble.dens_base")
    dens_cutoff = rp.get_param("bubble.dens_cutoff")

    x_pert = rp.get_param("bubble.x_pert")
    y_pert = rp.get_param("bubble.y_pert")
    r_pert = rp.get_param("bubble.r_pert")
    pert_amplitude_factor = rp.get_param("bubble.pert_amplitude_factor")

    # 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[:, :] = dens_cutoff

    # set the density to be stratified in the y-direction
    myg = my_data.grid

    p = myg.scratch_array()

    cs2 = scale_height*abs(grav)

    for j in range(myg.jlo, myg.jhi+1):
        dens[:, j] = max(dens_base*np.exp(-myg.y[j]/scale_height),
                        dens_cutoff)
        if j == myg.jlo:
            p[:, j] = dens[:, j]*cs2
        else:
            p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav

    # set the energy (P = cs2*dens)
    ener[:, :] = p[:, :]/(gamma - 1.0) + \
                0.5*(xmom[:, :]**2 + ymom[:, :]**2)/dens[:, :]

    r = np.sqrt((myg.x2d - x_pert)**2 + (myg.y2d - y_pert)**2)
    idx = r <= r_pert

    # boost the specific internal energy, keeping the pressure
    # constant, by dropping the density
    eint = (ener[idx] - 0.5*(xmom[idx]**2 - ymom[idx]**2)/dens[idx])/dens[idx]

    pres = dens[idx]*eint*(gamma - 1.0)

    eint = eint*pert_amplitude_factor
    dens[idx] = pres/(eint*(gamma - 1.0))

    ener[idx] = dens[idx]*eint + 0.5*(xmom[idx]**2 + ymom[idx]**2)/dens[idx]