Merge branch 'master' of https://github.com/zingale/pyro2

Author: Michael Zingale

compressible/BC.py

@@ -16,7 +16,8 @@ def user(bcName, bcEdge, variable, my_data):
     ymom = my_data.get_var("y-momentum")
     ener = my_data.get_var("energy")
 
-    grav = my_data.rp.get_param("compressible.grav")
+    grav = my_data.get_aux("grav")
+    gamma = my_data.get_aux("gamma")
 
     myg = my_data.grid
 
@@ -52,15 +53,15 @@ def user(bcName, bcEdge, variable, my_data):
                     dens[:,myg.jlo]
 
                 eint_base = (ener[:,myg.jlo] - ke_base)/dens[:,myg.jlo]
-                pres_base = eos.pres(dens_base, eint_base)
+                pres_base = eos.pres(gamma, dens_base, eint_base)
 
                 # we are assuming that the density is constant in this
                 # formulation of HSE, so the pressure comes simply from
                 # differencing the HSE equation
                 j = myg.jlo-1
                 while (j >= 0):
                     pres_below = pres_base - grav*dens_base*myg.dy
-                    rhoe = eos.rhoe(pres_below)
+                    rhoe = eos.rhoe(gamma, pres_below)
 
                     ener[:,j] = rhoe + ke_base
 
@@ -102,15 +103,15 @@ def user(bcName, bcEdge, variable, my_data):
                     dens[:,myg.jhi]
 
                 eint_base = (ener[:,myg.jhi] - ke_base)/dens[:,myg.jhi]
-                pres_base = eos.pres(dens_base, eint_base)
+                pres_base = eos.pres(gamma, dens_base, eint_base)
 
                 # we are assuming that the density is constant in this
                 # formulation of HSE, so the pressure comes simply from
                 # differencing the HSE equation
                 j = myg.jhi+1
                 while (j <= myg.jhi+myg.ng):
                     pres_above = pres_base + grav*dens_base*myg.dy
-                    rhoe = eos.rhoe(pres_above)
+                    rhoe = eos.rhoe(gamma, pres_above)
 
                     ener[:,j] = rhoe + ke_base
 

compressible/eos.py

@@ -1,36 +1,24 @@
-gamma = None
+"""
+This is a gamma-law equation of state.  
+"""
 
-def init(gamma_constant):
-    global gamma
-    gamma = gamma_constant
-
-
-def pres(dens, eint):
+def pres(gamma, dens, eint):
     """ given the density and the specific internal energy, return the
         pressure """
-
-    global gamma
     p = dens*eint*(gamma - 1.0)
-
     return p
 
 
-def dens(pres, eint):
+def dens(gamma, pres, eint):
     """ given the pressure and the specific internal energy, return
         the density """
-
-    global gamma
     dens = pres/(eint*(gamma - 1.0))
-
     return dens
 
 
-def rhoe(pres):
+def rhoe(gamma, pres):
     """ given the pressure, return (rho * e) """
-
-    global gamma
     rhoe = pres/(gamma - 1.0)
-
     return rhoe
 
 

compressible/simulation.py

@@ -113,12 +113,10 @@ def initialize(self):
 
 
         # store the EOS gamma as an auxillary quantity so we can have a
-        # self-contained object stored in output files to make plots
-        gamma = self.rp.get_param("eos.gamma")
-        my_data.set_aux("gamma", gamma)
-
-        # initialize the EOS gamma
-        eos.init(gamma)
+        # self-contained object stored in output files to make plots.
+        # store grav because we'll need that in some BCs
+        my_data.set_aux("gamma", self.rp.get_param("eos.gamma"))
+        my_data.set_aux("grav", self.rp.get_param("compressible.grav"))
 
         my_data.create()
 
@@ -160,11 +158,11 @@ def timestep(self):
 
         e = (ener - 0.5*dens*(u*u + v*v))/dens
 
-        p = eos.pres(dens, e)
-
-        # compute the sounds speed
         gamma = self.rp.get_param("eos.gamma")
 
+        p = eos.pres(gamma, dens, e)
+
+        # compute the sounds speed
         cs = numpy.sqrt(gamma*p/dens)
 
 
@@ -246,12 +244,13 @@ def dovis(self):
 
         magvel = numpy.sqrt(magvel)
 
-        # access gamma from the object instead of using the EOS so we can
-        # use dovis outside of a running simulation.
+        e = rhoe/dens
+
+        # access gamma from the cc_data object so we can use dovis 
+        # outside of a running simulation.
         gamma = self.cc_data.get_aux("gamma")
-        p = rhoe*(gamma - 1.0)
 
-        e = rhoe/dens
+        p = eos.pres(gamma, dens, e)
 
         myg = self.cc_data.grid
 

compressible/unsplitFluxes.py

@@ -146,6 +146,8 @@ def unsplitFluxes(my_data, rp, vars, tc, dt):
 
     myg = my_data.grid
 
+    gamma = rp.get_param("eos.gamma")
+
     #=========================================================================
     # compute the primitive variables
     #=========================================================================
@@ -165,7 +167,7 @@ def unsplitFluxes(my_data, rp, vars, tc, dt):
     # get the pressure
     e = (ener - 0.5*(xmom**2 + ymom**2)/dens)/dens
 
-    p = eos.pres(dens, e)
+    p = eos.pres(gamma, dens, e)
 
     smallp = 1.e-10
     p = p.clip(smallp)   # apply a floor to the pressure
@@ -218,8 +220,6 @@ def unsplitFluxes(my_data, rp, vars, tc, dt):
     tm_states = tc.timer("interfaceStates")
     tm_states.begin()
 
-    gamma = rp.get_param("eos.gamma")
-
     V_l = myg.scratch_array(vars.nvar)
     V_r = myg.scratch_array(vars.nvar)
 
@@ -239,13 +239,13 @@ def unsplitFluxes(my_data, rp, vars, tc, dt):
     U_xl[:,:,vars.idens] = V_l[:,:,vars.irho]
     U_xl[:,:,vars.ixmom] = V_l[:,:,vars.irho]*V_l[:,:,vars.iu]
     U_xl[:,:,vars.iymom] = V_l[:,:,vars.irho]*V_l[:,:,vars.iv]
-    U_xl[:,:,vars.iener] = eos.rhoe(V_l[:,:,vars.ip]) + \
+    U_xl[:,:,vars.iener] = eos.rhoe(gamma, V_l[:,:,vars.ip]) + \
         0.5*V_l[:,:,vars.irho]*(V_l[:,:,vars.iu]**2 + V_l[:,:,vars.iv]**2)
 
     U_xr[:,:,vars.idens] = V_r[:,:,vars.irho]
     U_xr[:,:,vars.ixmom] = V_r[:,:,vars.irho]*V_r[:,:,vars.iu]
     U_xr[:,:,vars.iymom] = V_r[:,:,vars.irho]*V_r[:,:,vars.iv]
-    U_xr[:,:,vars.iener] = eos.rhoe(V_r[:,:,vars.ip]) + \
+    U_xr[:,:,vars.iener] = eos.rhoe(gamma, V_r[:,:,vars.ip]) + \
         0.5*V_r[:,:,vars.irho]*(V_r[:,:,vars.iu]**2 + V_r[:,:,vars.iv]**2)
 
 
@@ -271,7 +271,7 @@ def unsplitFluxes(my_data, rp, vars, tc, dt):
                                   vars.nvar,
                                   gamma,
                                   r, u, v, p,
-                                  ldelta_r, ldelta_u, ldelta_v, ldelta_p)                                    
+                                  ldelta_r, ldelta_u, ldelta_v, ldelta_p)  
 
     tm_states.end()
 
@@ -283,13 +283,13 @@ def unsplitFluxes(my_data, rp, vars, tc, dt):
     U_yl[:,:,vars.idens] = V_l[:,:,vars.irho]
     U_yl[:,:,vars.ixmom] = V_l[:,:,vars.irho]*V_l[:,:,vars.iu]
     U_yl[:,:,vars.iymom] = V_l[:,:,vars.irho]*V_l[:,:,vars.iv]
-    U_yl[:,:,vars.iener] = eos.rhoe(V_l[:,:,vars.ip]) + \
+    U_yl[:,:,vars.iener] = eos.rhoe(gamma, V_l[:,:,vars.ip]) + \
         0.5*V_l[:,:,vars.irho]*(V_l[:,:,vars.iu]**2 + V_l[:,:,vars.iv]**2)
 
     U_yr[:,:,vars.idens] = V_r[:,:,vars.irho]
     U_yr[:,:,vars.ixmom] = V_r[:,:,vars.irho]*V_r[:,:,vars.iu]
     U_yr[:,:,vars.iymom] = V_r[:,:,vars.irho]*V_r[:,:,vars.iv]
-    U_yr[:,:,vars.iener] = eos.rhoe(V_r[:,:,vars.ip]) + \
+    U_yr[:,:,vars.iener] = eos.rhoe(gamma, V_r[:,:,vars.ip]) + \
         0.5*V_r[:,:,vars.irho]*(V_r[:,:,vars.iu]**2 + V_r[:,:,vars.iv]**2)
 
 

diffusion/FLOW

@@ -13,9 +13,12 @@ def makeplot(myd, solverName, outfile):
 
     exec 'import ' + solverName + ' as solver'
 
+    sim = solver.Simulation(None, None)
+    sim.cc_data = myd
+
     pylab.figure(num=1, figsize=(8,4.5), dpi=100, facecolor='w')
 
-    solver.dovis(myd, 0)
+    sim.dovis(0)
     pylab.savefig(outfile)
     pylab.show()