Consolidate pressure routines, improve VL+CT & Grackle compatability

doc/source/devel/hydro_infrastructure.rst

@@ -9,8 +9,7 @@ the Enzo Layer that can be optionally used to implement other
 hydro/MHD methods. The infrastructure was used to implement other the
 VL + CT MHD solver.
 
-*Note: Currently, barotropic equations of state and compatibility with
-Grackle (for* ``primordial_chemistry > 0`` *) are not yet implemented
+*Note: Currently, barotropic equations of state are not yet implemented
 within the infrastucture. However they are mentioned throughout this
 guide and slots have been explicitly left open for them to be
 implemented within the framework.*
@@ -369,10 +368,6 @@ Returns whether the dual energy formalism is in use.
 Returns the ratio of the specific heats. This is only required to
 yield a reasonable value if the gas is not barotropic.
 
-*In the future, the interface will need to be revisited once Grackle
-is fully supported and it will be possible for gamma to vary
-spatially.*
-
 .. code-block:: c++
 
    enzo_float get_isothermal_sound_speed();
@@ -407,32 +402,34 @@ energy. If the equation of state is barotropic, this should do nothing.
 
    void pressure_from_integration(const EnzoEFltArrayMap &integration_map,
                                   const CelloArray<enzo_float, 3> &pressure,
-                                  int stale_depth);
+                                  int stale_depth,
+                                  bool ignore_grackle = false);
 
 This method computes the pressure from the integration quantities
 (stored in ``integration_map``) and stores the result in ``pressure``.
-
-*In principle this should wrap* ``EnzoComputePressure``, *but
-currently that is not the case. Some minor refactoring is needed to
-allow EnzoComputePressure to compute Pressure based on arrays
-specified in a* ``EnzoEFltArrayMap`` *object and we are holding off on
-this until we implement full support for Grackle. Currently, when the
-dual-energy_formalism is in use, pressure is simply computed from
-internal energy.*
+This wraps the ``EnzoComputePressure`` object whose default behavior
+is to use the Grackle-supplied routine for computing pressure when the
+simulation is configured to use ``EnzoMethodGrackle``. The
+``ignore_grackle`` parameter can be used to avoid using that routine (the
+parameter is meaningless if the Grackle routine would not otherwise
+get used). This parameter's primary purpose is to provide the option
+to suppress the effects of molecular hydrogen on the adiabatic index
+(when Grackle is configured with ``primordial_chemistry > 1``).
 
 .. code-block:: c++
 
    void primitive_from_integration
      (const EnzoEFltArrayMap &integration_map, EnzoEFltArrayMap &primitive_map,
-      int stale_depth, const std::vector<std::string> &passive_list);
+      int stale_depth, const std::vector<std::string> &passive_list,
+      bool ignore_grackle = false);
 
 This method is responsible for computing the primitive quantities (to
 be held in ``primitive_map``) from the integration quantities (stored
 in ``integration_map``).  Non-passive scalar quantities appearing in
 both ``integration_map`` and ``primitive_map`` are simply deepcopied
 and passive scalar quantities are converted from conserved-form to
 specific form. For a non-barotropic EOS, this also computes pressure
-(by calling ``EnzoEquationOfState::pressure_from_integration``)
+(by calling ``EnzoEquationOfState::pressure_from_integration``).
 
 *In the future, it might be worth considering making this into a subclass
 of Cello's ``Physics`` class. If that is done, it may be advisable to

doc/source/user/problem_method.rst

@@ -548,6 +548,72 @@ A sample Method for implementing forward-euler to solve the heat equation.
 
 Calls methods provided by the external Grackle 3.0 chemistry and
 cooling library.
+
+Compatability with hydro/mhd solvers
+------------------------------------
+
+The ``"grackle"`` method is compatible with both the ``"ppm"`` and the
+``"mhd_vlct"`` methods. The convention is to list the hydro method
+before ``"grackle"`` in the ``Field:list`` parameter.  This
+configuration performs advection and radiative cooling in an
+operator-split manner (*Note: there isn't currently support for
+performing radiative cooling during the predictor step of the
+VL+CT solver*).
+
+Integration with hydro-solvers is self-consistent when
+``Field:Grackle:primordial_chemistry`` has values of ``0`` or ``1``.
+However, the integration is somewhat inconsistent when the parameter
+exceeds ``1``. While users shouldn't be too concerned about this
+latter scenario unless they are simulating conditions where
+:math:`{\rm H}_2` makes up a significant fraction of the gas density,
+we describe the inconsistencies in greater detail below.
+
+When ``Field:Grackle:primordial_chemistry > 1``, the Grackle library
+explicitly models chemistry involving :math:`{\rm H}_2` and how it
+modifies the adiabtic index. Grackle's routines treat
+:math:`\gamma_0`, the "nominal adiabatic index" specified by
+``Field:gamma``, as the adiabatic index for all monatomic species
+(this should be ``5.0/3.0``). To that end, Grackle supplies functions
+that can effectively be represented as :math:`\gamma(e, n_{{\rm H}_2},
+n_{\rm other})` and :math:`p(\rho, e, n_{{\rm H}_2}, n_{\rm
+other})`. In these formulas:
+
+- :math:`p`, :math:`\rho` and :math:`e` correspond to the quantities
+  held by the ``pressure``, ``density`` and ``internal_energy``
+  fields.  *(Note: the* :math:`\gamma` *function's dependence on*
+  :math:`e` *accounts for the dependence of* :math:`\gamma_{{\rm
+  H}_2}` *on temperature)*
+
+- :math:`n_{{\rm H}_2}` specifies the number density of
+  :math:`{\rm H}_2`. :math:`n_{\rm other}` specifies a selection of
+  the other primordial species (that roughly approximate the total
+  number density). In practice, these are computed from passively
+  advected species fields.
+
+There are a handful of locations within the ``"ppm"`` and
+``"mhd_vlct"`` methods where this treatment is relevant:
+
+1. **Computing the timestep:** each hydro/mhd
+   method uses the :math:`p(\rho, e, n_{{\rm H}_2}, n_{\rm other})`
+   function for the pressure values. However, they both use
+   :math:`\gamma_0` in other places.
+
+2. **Pre-reconstruction pressure calculation:** each hydro/mhd
+   solver internally computes the pressure that is to be reconstructed
+   with :math:`p=(\gamma_0 - 1)e\rho`.
+
+3. **Riemann Solver:** in each hydro/mhd solver, the Riemann Solver
+   completely ignore the grackle supplied functions.
+
+4. **VL+CT Energy floor and DE synchronization:** the internal energy
+   floor is computed from the pressure floor using: :math:`e_{\rm
+   floor} = \frac{p_{\rm floor}}{(\gamma_0 - 1)\rho}` (thus,
+   :math:`p_{\rm floor}` may exceed :math:`p(\rho, e_{\rm floor},
+   \ldots)`). Additionally, synchronizing the internal energy with
+   total energy relies on :math:`\gamma_0`.
+
+5. **PPM reconstruction:** uses :math:`\gamma_0`.
+
 
 ``"comoving_expansion"``: comoving expansion
 ============================================

input/Grackle/method_grackle_shock_tube.in

@@ -0,0 +1,74 @@
+# this file is used to just simply check if Enzo-E can run a hydro solver with
+# Grackle without crashing (in the future, this should be replaced with a more
+# meaningful test problem)
+
+include "input/vlct/vlct_de.incl"
+include "input/vlct/dual_energy_shock_tube/initial_sod_shock_tube.incl"
+
+
+   Domain {
+      lower = [0.0, 0.0, 0.0];
+      upper = [1.0, 0.03125, 0.03125];
+   }
+
+   Units {
+      density = 1.6726219E-24;  # m_H in grams
+      time    = 3.15576E13;     # 1 Myr in seconds
+      length  = 3.086E16;       # 0.01 pc in cm - does not actually matter
+   }
+
+   Mesh {
+      root_rank = 3; # 3D
+      root_blocks = [1,1,1];
+      root_size = [128,4,4]; # number of cells per axis
+   }
+
+   Boundary {
+      list = ["two", "three","one"];
+      two{
+         type = "periodic";
+   	 axis = "y";
+      }
+      three{
+         type = "periodic";
+   	 axis = "z";
+      }
+      one{
+         type = "outflow";
+   	 axis = "x";
+      }
+   }
+
+   Output {
+     list = [];
+   }
+
+   Field {
+       gamma = 5./3.;
+   }
+
+      Method {
+      list = [ "mhd_vlct", "grackle"];
+
+      grackle {
+        courant = 0.40; # meaningless unless use_cooling_timestep = true;
+
+        data_file = "CloudyData_UVB=HM2012_shielded.h5";
+
+        with_radiative_cooling = 1;
+        primordial_chemistry   = 3;  # 1, 2, or 3
+        metal_cooling          = 1;  # 0 or 1 (off/on)
+        UVbackground           = 1;  # on or off
+        self_shielding_method  = 0;  # 0 - 3 (0 or 3 recommended)
+
+        HydrogenFractionByMass = 0.73;
+
+        # set this to true to limit the maximum timestep to the product of the
+        # minimum cooling/heating time and courant.
+        use_cooling_timestep = true; # default is false
+      }
+   }
+
+
+   Stopping { cycle = 200; }
+

src/Enzo/_enzo.hpp

@@ -177,6 +177,7 @@ extern "C" {
 #include "enzo_IoEnzoBlock.hpp"
 
 #include "enzo_EnzoEFltArrayMap.hpp"
+#include "enzo_EnzoFieldAdaptor.hpp"
 #include "enzo_EnzoPermutedCoordinates.hpp"
 #include "enzo_EnzoCenteredFieldRegistry.hpp"
 

src/Enzo/enzo.cpp

@@ -43,4 +43,23 @@ namespace enzo {
     return static_cast<EnzoBlock*> (block);
   }
 
+  bool uses_dual_energy_formalism(bool default_ret /* false */)
+  {
+    // TODO(mabruzzo): this is meant to be a short-term solution. My immediate
+    // priority is to create a Physics object for storing EOS properties
+    // (including dual energy formalism parameters). This function will be
+    // modified or deleted at that time.
+
+    EnzoProblem* prob = problem();
+    if (prob->method("ppm")){
+      return config()->ppm_dual_energy;
+    } else if (prob->method("mhd_vlct")){
+      return config()->method_vlct_dual_energy;
+    } else if (prob->method("ppml")){
+      return false;
+    } else {
+      return default_ret;
+    }
+  }
+
 }

src/Enzo/enzo.hpp

@@ -64,6 +64,12 @@ namespace enzo {
   const EnzoConfig * config();
   CProxy_EnzoBlock block_array();
   EnzoBlock * block ( Block * block);
+
+  /// Returns whether the dual energy formalism is in use.
+  ///
+  /// @param default_ret[in] The value to return if no hydro methods are used.
+  ///     The default value is false.
+  bool uses_dual_energy_formalism(bool default_ret = false);
 }
 
 extern CProxy_EnzoSimulation proxy_enzo_simulation;

src/Enzo/enzo_EnzoComputeCoolingTime.cpp

@@ -73,8 +73,8 @@ void EnzoComputeCoolingTime::compute_(Block * block,
          "Grackle must be enabled in order to compute the cooling time",
          enzo::config()->method_grackle_use_grackle );
   const EnzoMethodGrackle* grackle_method = enzo::grackle_method();
-  grackle_method->calculate_cooling_time(block, ct, grackle_units,
-					 grackle_fields, i_hist_);
+  grackle_method->calculate_cooling_time(EnzoFieldAdaptor(block, i_hist_), ct,
+                                         0, grackle_units, grackle_fields);
 }
 
 #endif