compressible: add the ability for a problem-dependent external source (…

…#289)

this also adds a "heating" test problem.

pyro/compressible/problems/heating.py

@@ -0,0 +1,65 @@
+"""A test of the energy sources.  This uses a uniform domain and
+slowly adds heat to the center over time."""
+
+import numpy as np
+
+from pyro.util import msg
+
+DEFAULT_INPUTS = "inputs.heating"
+
+PROBLEM_PARAMS = {"heating.rho_ambient": 1.0,  # ambient density
+                  "heating.p_ambient": 10.0,  # ambient pressure
+                  "heating.r_src": 0.1,  # physical size of the heating src
+                  "heating.e_rate": 0.1}  # energy generation rate (energy / mass / time)
+
+
+def init_data(my_data, rp):
+    """ initialize the heating problem """
+
+    if rp.get_param("driver.verbose"):
+        msg.bold("initializing the heating 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")
+
+    # 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[:, :] = rp.get_param("heating.rho_ambient")
+    xmom[:, :] = 0.0
+    ymom[:, :] = 0.0
+    ener[:, :] = rp.get_param("heating.p_ambient") / (gamma - 1.0)
+
+
+def source_terms(myg, U, ivars, rp):
+    """source terms to be added to the evolution"""
+
+    S = myg.scratch_array(nvar=ivars.nvar)
+
+    xctr = 0.5 * (myg.xmin + myg.xmax)
+    yctr = 0.5 * (myg.ymin + myg.ymax)
+
+    dist = np.sqrt((myg.x2d - xctr)**2 +
+                   (myg.y2d - yctr)**2)
+
+    e_rate = rp.get_param("heating.e_rate")
+    r_src = rp.get_param("heating.r_src")
+
+    S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_src)**2)
+    return S
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """
+
+    print("""
+          The script analysis/sedov_compare.py can be used to analyze these
+          results.  That will perform an average at constant radius and
+          compare the radial profiles to the exact solution.  Sample exact
+          data is provided as analysis/cylindrical-sedov.out
+          """)

pyro/compressible/problems/inputs.heating

@@ -0,0 +1,41 @@
+[driver]
+max_steps = 5000
+tmax = 1.0
+
+
+[compressible]
+limiter = 2
+cvisc = 0.1
+
+
+[io]
+basename = heating_
+dt_out = 0.1
+
+
+[eos]
+gamma = 1.4
+
+
+[mesh]
+nx = 64
+ny = 64
+xmax = 1.0
+ymax = 1.0
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = outflow
+yrboundary = outflow
+
+
+[heating]
+rho_ambient = 1.0
+p_ambient = 10.0
+r_src = 0.05
+e_rate = 0.1
+
+
+[vis]
+dovis = 1

pyro/compressible/problems/inputs.plume

@@ -0,0 +1,39 @@
+# simple inputs files for the four-corner problem.
+
+[driver]
+max_steps = 10000
+tmax = 10.0
+
+
+[io]
+basename = plume_
+n_out = 100
+
+
+[mesh]
+nx = 128
+ny = 256
+xmax = 4.0
+ymax = 8.0
+
+xlboundary = outflow
+xrboundary = outflow
+
+ylboundary = hse
+yrboundary = hse
+
+
+[plume]
+scale_height = 3.0
+dens_base = 1000.0
+
+x_pert = 2.0
+y_pert = 2.0
+r_pert = 0.25
+e_rate = 0.5
+
+
+[compressible]
+grav = -2.0
+
+limiter = 2

pyro/compressible/problems/plume.py

@@ -0,0 +1,90 @@
+"""A heat source at a point creates a plume that buoynantly rises in
+an adiabatically stratified atmosphere."""
+
+import numpy as np
+
+from pyro.util import msg
+
+DEFAULT_INPUTS = "inputs.plume"
+
+PROBLEM_PARAMS = {"plume.dens_base": 10.0,  # density at the base of the atmosphere
+                  "plume.scale_height": 4.0,  # scale height of the isothermal atmosphere
+                  "plume.x_pert": 2.0,
+                  "plume.y_pert": 2.0,
+                  "plume.r_pert": 0.25,
+                  "plume.e_rate": 0.1,
+                  "plume.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("plume.scale_height")
+    dens_base = rp.get_param("plume.dens_base")
+    dens_cutoff = rp.get_param("plume.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
+        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[:, :]
+
+
+def source_terms(myg, U, ivars, rp):
+    """source terms to be added to the evolution"""
+
+    S = myg.scratch_array(nvar=ivars.nvar)
+
+    x_pert = rp.get_param("plume.x_pert")
+    y_pert = rp.get_param("plume.y_pert")
+
+    dist = np.sqrt((myg.x2d - x_pert)**2 +
+                   (myg.y2d - y_pert)**2)
+
+    e_rate = rp.get_param("plume.e_rate")
+    r_pert = rp.get_param("plume.r_pert")
+
+    S[:, :, ivars.iener] = U[:, :, ivars.idens] * e_rate * np.exp(-(dist / r_pert)**2)
+    return S
+
+
+def finalize():
+    """ print out any information to the user at the end of the run """

pyro/compressible/simulation.py

@@ -91,7 +91,7 @@ def prim_to_cons(q, gamma, ivars, myg):
     return U
 
 
-def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None):
+def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None, problem_source=None):
     """compute the external sources, including gravity"""
 
     _ = t  # maybe unused
@@ -142,6 +142,11 @@ def get_external_sources(t, dt, U, ivars, rp, myg, *, U_old=None):
 
             S[:, :, ivars.iener] = ymom_new * grav
 
+    # now add the heating
+    if problem_source:
+        S_heating = problem_source(myg, U, ivars, rp)
+        S[...] += S_heating
+
     return S
 
 
@@ -274,7 +279,8 @@ def evolve(self):
         # Only gravitional source for Cartesian2d
         U_xl, U_xr, U_yl, U_yr = flx.apply_source_terms(U_xl, U_xr, U_yl, U_yr,
                                                         self.cc_data, self.aux_data, self.rp,
-                                                        self.ivars, self.tc, self.dt)
+                                                        self.ivars, self.tc, self.dt,
+                                                        problem_source=self.problem_source)
 
         # Apply transverse corrections.
         U_xl, U_xr, U_yl, U_yr = flx.apply_transverse_flux(U_xl, U_xr, U_yl, U_yr,
@@ -361,7 +367,8 @@ def evolve(self):
         # * correct: U^{n+1} = U^{n+1,*} + dt/2 (S^{n+1} - S^n)
 
         S_old = get_external_sources(self.cc_data.t, self.dt, U_old,
-                                     self.ivars, self.rp, myg)
+                                     self.ivars, self.rp, myg,
+                                     problem_source=self.problem_source)
 
         for n in range(self.ivars.nvar):
             var = self.cc_data.get_var_by_index(n)
@@ -370,7 +377,8 @@ def evolve(self):
         # now get the new time source
 
         S_new = get_external_sources(self.cc_data.t, self.dt, self.cc_data.data,
-                                     self.ivars, self.rp, myg, U_old=U_old)
+                                     self.ivars, self.rp, myg, U_old=U_old,
+                                     problem_source=self.problem_source)
 
         # and correct
         for n in range(self.ivars.nvar):

pyro/compressible/unsplit_fluxes.py

@@ -243,7 +243,8 @@ def interface_states(my_data, rp, ivars, tc, dt):
 
 
 def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
-                       my_data, my_aux, rp, ivars, tc, dt):
+                       my_data, my_aux, rp, ivars, tc, dt, *,
+                       problem_source=None):
     """
     This function applies source terms including external (gravity),
     geometric terms, and pressure terms to the left and right
@@ -270,6 +271,8 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
         The timers we are using to profile
     dt : float
         The timestep we are advancing through.
+    problem_source : function (optional)
+        A problem-specific source function to call
 
     Returns
     -------
@@ -290,7 +293,8 @@ def apply_source_terms(U_xl, U_xr, U_yl, U_yr,
     ymom_src = my_aux.get_var("ymom_src")
     E_src = my_aux.get_var("E_src")
 
-    U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg)
+    U_src = comp.get_external_sources(my_data.t, dt, my_data.data, ivars, rp, myg,
+                                      problem_source=problem_source)
 
     dens_src[:, :] = U_src[:, :, ivars.idens]
     xmom_src[:, :] = U_src[:, :, ivars.ixmom]

pyro/compressible_fv4/simulation.py

@@ -7,9 +7,11 @@
 class Simulation(compressible_rk.Simulation):
 
     def __init__(self, solver_name, problem_name, problem_func, rp, *,
-                 problem_finalize_func=None, timers=None, data_class=fv.FV2d):
+                 problem_finalize_func=None, problem_source_func=None,
+                 timers=None, data_class=fv.FV2d):
         super().__init__(solver_name, problem_name, problem_func, rp,
                          problem_finalize_func=problem_finalize_func,
+                         problem_source_func=problem_source_func,
                          timers=timers, data_class=data_class)
 
     def substep(self, myd):
@@ -29,7 +31,8 @@ def substep(self, myd):
 
         # cell-centered sources
         S = get_external_sources(myd.t, self.dt, U_cc,
-                                 self.ivars, self.rp, myg)
+                                 self.ivars, self.rp, myg,
+                                 problem_source=self.problem_source)
 
         # bring the sources back to averages -- we only care about
         # the interior (no ghost cells)

pyro/compressible_rk/simulation.py

@@ -20,7 +20,8 @@ def substep(self, myd):
         # source terms -- note: this dt is the entire dt, not the
         # stage's dt
         S = compressible.get_external_sources(myd.t, self.dt, myd.data,
-                                              self.ivars, self.rp, myg)
+                                              self.ivars, self.rp, myg,
+                                              problem_source=self.problem_source)
 
         k = myg.scratch_array(nvar=self.ivars.nvar)
 

pyro/lm_atm/simulation.py

@@ -36,10 +36,13 @@ def jp(self, shift, buf=0):
 class Simulation(NullSimulation):
 
     def __init__(self, solver_name, problem_name, problem_func, rp, *,
-                 problem_finalize_func=None, timers=None):
+                 problem_finalize_func=None, problem_source_func=None,
+                 timers=None):
 
         super().__init__(solver_name, problem_name, problem_func, rp,
-                         problem_finalize_func=problem_finalize_func, timers=timers)
+                         problem_finalize_func=problem_finalize_func,
+                         problem_source_func=problem_source_func,
+                         timers=timers)
 
         self.base = {}
         self.aux_data = None