[WIP] add a compressible convection problem (#293)

this also adds an ambient upper boundary condition

pyro/compressible/BC.py

@@ -141,6 +141,45 @@ def user(bc_name, bc_edge, variable, ccdata):
         else:
             msg.fail("error: hse BC not supported for xlb or xrb")
 
+    elif bc_name == "ambient":
+
+        # we'll assume that the problem setup defines these
+        # they will not be available for source terms
+        ambient_rho = ccdata.get_aux("ambient_rho")
+        ambient_u = ccdata.get_aux("ambient_u")
+        ambient_v = ccdata.get_aux("ambient_v")
+        ambient_p = ccdata.get_aux("ambient_p")
+
+        if bc_edge == "yrb":
+
+            # upper y boundary
+
+            # by default, use a zero gradient
+            v = ccdata.get_var(variable)
+            for j in range(myg.jhi+1, myg.jhi+myg.ng+1):
+                v[:, j] = v[:, myg.jhi]
+
+            # overwrite with ambient conditions
+            if variable == "density":
+                v[:, myg.jhi+1:myg.jhi+myg.ng+1] = ambient_rho
+
+            elif variable == "x-momentum":
+                rhou = ambient_rho * ambient_u
+                v[:, myg.jhi+1:myg.jhi+myg.ng+1] = rhou
+
+            elif variable == "y-momentum":
+                rhov = ambient_rho * ambient_v
+                v[:, myg.jhi+1:myg.jhi+myg.ng+1] = rhov
+
+            elif variable == "energy":
+                gamma = ccdata.get_aux("gamma")
+                ke = 0.5 * ambient_rho * (ambient_u**2 + ambient_v**2)
+                rhoE = ambient_p / (gamma - 1.0) + ke
+                v[:, myg.jhi+1:myg.jhi+myg.ng+1] = rhoE
+
+        else:
+            msg.fail("error: ambient BC not supported for xlb, xrb, or ylb")
+
     elif bc_name == "ramp":
         # Boundary conditions for double Mach reflection problem
 

pyro/compressible/problems/convection.py

@@ -0,0 +1,110 @@
+"""A heat source in a layer at some height above the bottom will drive
+convection in an adiabatically stratified atmosphere."""
+
+import numpy as np
+
+from pyro.util import msg
+
+DEFAULT_INPUTS = "inputs.convection"
+
+PROBLEM_PARAMS = {"convection.dens_base": 10.0,  # density at the base of the atmosphere
+                  "convection.scale_height": 4.0,  # scale height of the isothermal atmosphere
+                  "convection.y_height": 2.0,
+                  "convection.thickness": 0.25,
+                  "convection.e_rate": 0.1,
+                  "convection.dens_cutoff": 0.01}
+
+
+def init_data(my_data, rp):
+    """ initialize the bubble problem """
+
+    if rp.get_param("driver.verbose"):
+        msg.bold("initializing the bubble 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")
+
+    gamma = rp.get_param("eos.gamma")
+
+    grav = rp.get_param("compressible.grav")
+
+    scale_height = rp.get_param("convection.scale_height")
+    dens_base = rp.get_param("convection.dens_base")
+    dens_cutoff = rp.get_param("convection.dens_cutoff")
+
+    # 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()
+
+    pres_base = scale_height*dens_base*abs(grav)
+
+    for j in range(myg.jlo, myg.jhi+1):
+        profile = 1.0 - (gamma-1.0)/gamma * myg.y[j]/scale_height
+        if profile > 0.0:
+            dens[:, j] = max(dens_base*(profile)**(1.0/(gamma-1.0)),
+                             dens_cutoff)
+        else:
+            dens[:, j] = dens_cutoff
+
+        if j == myg.jlo:
+            p[:, j] = pres_base
+        elif dens[0, j] <= dens_cutoff + 1.e-30:
+            p[:, j] = p[:, j-1]
+        else:
+            #p[:, j] = p[:, j-1] + 0.5*myg.dy*(dens[:, j] + dens[:, j-1])*grav
+            p[:, j] = pres_base * (dens[:, j] / dens_base)**gamma
+
+    # set the ambient conditions
+    my_data.set_aux("ambient_rho", dens_cutoff)
+    my_data.set_aux("ambient_u", 0.0)
+    my_data.set_aux("ambient_v", 0.0)
+    my_data.set_aux("ambient_p", p.v().min())
+
+    # set the energy (P = cs2*dens) -- assuming zero velocity
+    ener[:, :] = p[:, :]/(gamma - 1.0)
+
+    # pairs of random numbers between [-1, 1]
+    vel_pert = 2.0 * np.random.random_sample((myg.qx, myg.qy, 2)) - 1
+
+    cs = np.sqrt(gamma * p / dens)
+
+    # make vel_pert have M < 0.05
+    vel_pert[:, :, 0] *= 0.05 * cs
+    vel_pert[:, :, 1] *= 0.05 * cs
+
+    idx = dens > 2 * dens_cutoff
+    xmom[idx] = dens[idx] * vel_pert[idx, 0]
+    ymom[idx] = dens[idx] * vel_pert[idx, 1]
+
+    ener[:, :] += 0.5 * (xmom[:, :]**2 + ymom[:, :]**2) / dens[:, :]
+
+
+def source_terms(myg, U, ivars, rp):
+    """source terms to be added to the evolution"""
+
+    S = myg.scratch_array(nvar=ivars.nvar)
+
+    y_height = rp.get_param("convection.y_height")
+
+    dist = np.abs(myg.y2d - y_height)
+
+    e_rate = rp.get_param("convection.e_rate")
+    thick = rp.get_param("convection.thickness")
+
+    S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / thick)**2)
+    return S
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """

pyro/compressible/problems/inputs.convection

@@ -0,0 +1,37 @@
+# simple inputs files for the four-corner problem.
+
+[driver]
+max_steps = 10000
+tmax = 10.0
+
+
+[io]
+basename = convection_
+n_out = 100
+
+
+[mesh]
+nx = 128
+ny = 384
+xmax = 4.0
+ymax = 12.0
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = reflect
+yrboundary = ambient
+
+
+[convection]
+scale_height = 2.0
+dens_base = 1000.0
+dens_cutoff = 1.e-3
+
+e_rate = 0.5
+
+
+[compressible]
+grav = -2.0
+
+limiter = 2

pyro/compressible/simulation.py

@@ -176,6 +176,7 @@ def initialize(self, *, extra_vars=None, ng=4):
 
         # define solver specific boundary condition routines
         bnd.define_bc("hse", BC.user, is_solid=False)
+        bnd.define_bc("ambient", BC.user, is_solid=False)
         bnd.define_bc("ramp", BC.user, is_solid=False)  # for double mach reflection problem
 
         bc, bc_xodd, bc_yodd = bc_setup(self.rp)
@@ -489,3 +490,4 @@ def write_extras(self, f):
 
         # the value here is the value of "is_solid"
         gb.create_dataset("hse", data=False)
+        gb.create_dataset("ambient", data=False)