Updates to EnzoRiemann related to the shift to using primitives to co…

…mpute the fluxes (and the new "integration quantity" terminology)

Other changes include:
- Adjusted the max asymmetry residual in the Dual Energy Cloud Test for the HLLC Riemann Solver.
- I transitioned from registering the integration quantities with the Riemann Solver (in the Riemann Solver factory method) to having the Riemann Solver directly specify its expected primitives and integrations quantities. This seems to make much more sense than independently coming up with lists of integration and primitive quanities in EnzoMethodMHDVlct and making sure that they are consistent with the Riemann Solver's requirements. This change includes tweaks to EnzoRiemann::construct_riemann and the introduction of the EnzoRiemann::integration_quantities() and EnzoRiemann::primitive_quantities() interface methods.
- This change made the EnzoMethodMHDVlct::determine_quantities_ helper method unnecessary, so I deleted it.
- I also deleted the now unnecessary method from EnzoEquationOfState, eint_from_primitive.

Finally, I updated the website documentation for EnzoRiemann to reflect these changes

input/vlct/run_dual_energy_cloud_test.py

@@ -83,7 +83,7 @@ def analyze_tests():
     r += check_cloud_asym('hlld_cloud_0.0625/hlld_cloud_0.0625.block_list',
                           'hlld_cloud', 5.5e-13)
     r += check_cloud_asym('hllc_cloud_0.0625/hllc_cloud_0.0625.block_list',
-                          'hllc_cloud', 3.6e-13)
+                          'hllc_cloud', 3.8e-13)
     r += check_cloud_asym('hlle_cloud_0.0625/hlle_cloud_0.0625.block_list',
                           'hlle_cloud', 3.e-13)
     n_passed = np.sum(r)

src/Enzo/enzo_EnzoEOSIdeal.cpp

@@ -186,34 +186,6 @@ void EnzoEOSIdeal::pressure_from_integration
 
 //----------------------------------------------------------------------
 
-void EnzoEOSIdeal::eint_from_primitive(EnzoEFltArrayMap &primitive,
-                                       EFlt3DArray &eint,
-                                       int stale_depth) const
-{
-  if (grackle_variable_gamma_()){
-    // we would need to model a field with the
-    ERROR("EnzoEOSIdeal::eint_from_primitive",
-	  "Doesn't currently support spatial variations in gamma");
-  }
-
-  EFlt3DArray density = primitive.get("density", stale_depth);
-  EFlt3DArray pressure = primitive.get("pressure", stale_depth);
-
-  CSlice unstaled(stale_depth,-stale_depth);
-  EFlt3DArray e = eint.subarray(unstaled, unstaled, unstaled);
-  enzo_float inv_gm1 = 1./(get_gamma()-1.);
-
-  for (int iz=0; iz<e.shape(0); iz++) {
-    for (int iy=0; iy<e.shape(1); iy++) {
-      for (int ix=0; ix<e.shape(2); ix++) {
-        e(iz,iy,ix) = pressure(iz,iy,ix) * inv_gm1 / density(iz,iy,ix);
-      }
-    }
-  }
-}
-
-//----------------------------------------------------------------------
-
 // based on the enzo's hydro_rk implementation of synchronization (found in the
 // Grid_UpdateMHD.C file)
 void EnzoEOSIdeal::apply_floor_to_energy_and_sync

src/Enzo/enzo_EnzoEOSIdeal.hpp

@@ -58,9 +58,6 @@ class EnzoEOSIdeal : public EnzoEquationOfState
   (EnzoEFltArrayMap &integration_map, const EFlt3DArray &pressure,
    int stale_depth) const;
 
-  void eint_from_primitive(EnzoEFltArrayMap &primitive, EFlt3DArray &eint,
-                           int stale_depth) const;
-
   inline enzo_float get_density_floor() const { return density_floor_; }
 
   enzo_float get_pressure_floor() const { return pressure_floor_; }

src/Enzo/enzo_EnzoEquationOfState.hpp

@@ -129,21 +129,6 @@ class EnzoEquationOfState : public PUP::able
   (EnzoEFltArrayMap &integration_map, const EFlt3DArray &pressure,
    int stale_depth) const = 0;
 
-  /// Computes specific internal energy from primitive quantities (nominally
-  /// after reconstruction)
-  ///
-  /// @param[in]  primitive Map holding primitive values that are used to
-  ///     compute the specific internal energy
-  /// @param[out] pressure Array where the specific internal energy is to be
-  ///     stored.
-  /// @param[in]  stale_depth indicates the current stale_depth for the
-  ///     supplied cell-centered quantities
-  ///
-  /// For a barotropic EOS, this might not do anything
-  virtual void eint_from_primitive(EnzoEFltArrayMap &primitive,
-                                   EFlt3DArray &eint,
-                                   int stale_depth) const = 0;
-
   /// returns the density floor
   virtual enzo_float get_density_floor() const = 0;
 

src/Enzo/enzo_EnzoMethodMHDVlct.cpp

@@ -59,22 +59,35 @@ EnzoMethodMHDVlct::EnzoMethodMHDVlct (std::string rsolver,
   eos_ = new EnzoEOSIdeal(gamma, density_floor, pressure_floor,
 			  dual_energy_formalism, dual_energy_formalism_eta);
 
+#ifdef CONFIG_USE_GRACKLE
+  if (enzo::config()->method_grackle_use_grackle){
+    // we can remove the following once EnzoMethodGrackle no longer requires
+    // the internal_energy to be a permanent field
+    ASSERT("EnzoMethodMHDVlct::determine_quantities_",
+           ("Grackle cannot currently be used alongside this integrator "
+            "unless the dual-energy formalism is in use"),
+           eos_->uses_dual_energy_formalism());
+  }
+#endif /* CONFIG_USE_GRACKLE */
+
   // Determine whether magnetic fields are to be used
   mhd_choice_ = parse_bfield_choice_(mhd_choice);
 
+  riemann_solver_ = EnzoRiemann::construct_riemann
+    (rsolver, mhd_choice_ != bfield_choice::no_bfield,
+     eos_->uses_dual_energy_formalism());
+
   // determine integration and primitive quantities
   std::vector<std::string> integration_quantities, primitive_quantities;
-  EnzoMethodMHDVlct::determine_quantities_(eos_, mhd_choice_,
-                                           integration_quantities,
-                                           primitive_quantities);
+  integration_quantities = riemann_solver_->integration_quantities();
+  primitive_quantities = riemann_solver_->primitive_quantities();
 
   // Initialize the remaining component objects
   half_dt_recon_ = EnzoReconstructor::construct_reconstructor
     (primitive_quantities, half_recon_name, (enzo_float)theta_limiter);
   full_dt_recon_ = EnzoReconstructor::construct_reconstructor
     (primitive_quantities, full_recon_name, (enzo_float)theta_limiter);
-  riemann_solver_ = EnzoRiemann::construct_riemann(integration_quantities,
-                                                   rsolver);
+
   integration_quan_updater_ =
     new EnzoIntegrationQuanUpdate(integration_quantities, true);
 
@@ -131,51 +144,6 @@ EnzoMethodMHDVlct::bfield_choice EnzoMethodMHDVlct::parse_bfield_choice_
 
 //----------------------------------------------------------------------
 
-void EnzoMethodMHDVlct::determine_quantities_
-(const EnzoEquationOfState *eos, EnzoMethodMHDVlct::bfield_choice mhd_choice,
- std::vector<std::string> &integration_quantities,
- std::vector<std::string> &primitive_quantities) noexcept
-{
-#ifdef CONFIG_USE_GRACKLE
-  if (enzo::config()->method_grackle_use_grackle){
-    // we can remove the following once EnzoMethodGrackle no longer requires
-    // the internal_energy to be a permanent field
-    ASSERT("EnzoMethodMHDVlct::determine_quantities_",
-           ("Grackle cannot currently be used alongside this integrator "
-            "unless the dual-energy formalism is in use"),
-           eos->uses_dual_energy_formalism());
-  }
-#endif /* CONFIG_USE_GRACKLE */
-
-  std::string common[] {"density", "velocity"};
-  for (std::string quantity : common){
-    integration_quantities.push_back(quantity);
-    primitive_quantities.push_back(quantity);
-  }
-
-  if (mhd_choice != bfield_choice::no_bfield){
-    integration_quantities.push_back("bfield");
-    primitive_quantities.push_back("bfield");
-  }
-
-  if (eos->is_barotropic()){
-    ERROR("EnzoMethodMHDVlct::determine_quantities_",
-	  "Not presently equipped to handle barotropic equations of state.");
-  } else {
-    // add specific total energy to integration quantities
-    integration_quantities.push_back("total_energy");
-    // add pressure to primitive quantities
-    primitive_quantities.push_back("pressure");
-    // add specific internal energy to integration quantities (if using the
-    // dual energy formalism)
-    if (eos->uses_dual_energy_formalism()){
-      integration_quantities.push_back("internal_energy");
-    }
-  }
-}
-
-//----------------------------------------------------------------------
-
 EnzoMethodMHDVlct::~EnzoMethodMHDVlct()
 {
   delete eos_;

src/Enzo/enzo_EnzoMethodMHDVlct.hpp

@@ -145,23 +145,6 @@ class EnzoMethodMHDVlct : public Method {
   /// returns the bfield_choice enum that matches the input string
   bfield_choice parse_bfield_choice_(std::string choice) const noexcept;
 
-  /// Determines the integration quantities and primitives from (FIELD_TABLE)
-  /// that are to be integrated and reconstructed
-  ///
-  /// @param[in]  eos Pointer to the fluid's EquationOfState object.
-  /// @param[in]  mhd_choice Encodes how the integrator will handle B-fields
-  /// @param[out] integration_quantities Reference to a vector that get's filled
-  ///     by this function with the integration quantities (matching names in
-  ///     FIELD_TABLE) used by the integrator
-  /// @param[out] primitive_quantities Reference to a vector that get's
-  ///     filled by this function with the names of primitive quantities
-  ///     (matching names in FIELD_TABLE) that are used by the integrator for
-  ///     reconstruction
-  static void determine_quantities_
-  (const EnzoEquationOfState *eos, bfield_choice mhd_choice,
-   std::vector<std::string> &integration_quantities,
-   std::vector<std::string> &primitive_quantities) noexcept;
-
   /// Checks that the mesh size is sufficiently large to handle the given ghost
   /// depth and confirms that the ghost depth is consistent with the
   /// requirements of the reconstructors

src/Enzo/enzo_EnzoRiemann.cpp

@@ -12,30 +12,44 @@
 
 //----------------------------------------------------------------------
 
-EnzoRiemann* EnzoRiemann::construct_riemann
-(std::vector<std::string> integrable_quantities, std::string solver)
+EnzoRiemann* EnzoRiemann::construct_riemann(std::string solver, bool mhd,
+                                            bool internal_energy)
 {
   // determine the type of solver to construct:
   // convert string to lower case (https://stackoverflow.com/a/313990)
   std::string formatted(solver.size(), ' ');
   std::transform(solver.begin(), solver.end(), formatted.begin(),
 		 ::tolower);
-  EnzoRiemann* out;
+  EnzoRiemann* out = nullptr; // set to NULL to suppress compiler warnings
+
+  if (formatted == "hll"){
+    ASSERT("EnzoRiemann::construct_riemann",
+           ("An \"HLL\" Riemann solver without magnetic fields isn't "
+            "currently supported."), mhd);
+    out = new EnzoRiemannHLLMHD(internal_energy);
+
+  } else if (formatted == "hlle"){
+    if (mhd){
+      out = new EnzoRiemannHLLEMHD(internal_energy);
+    } else {
+      ERROR("EnzoRiemann::construct_riemann",
+            "The \"HLLE\" Riemann solver without magnetic fields is untested");
+      out = new EnzoRiemannHLLE(internal_energy);
+    }
 
-  // Eventually we may want to check for non-MHD Riemann solvers
-  if (formatted == std::string("hll")){
-    out = new EnzoRiemannHLLMHD(integrable_quantities);
-  } else if (formatted == std::string("hlle")){
-    out = new EnzoRiemannHLLEMHD(integrable_quantities);
   } else if (formatted == std::string("hllc")){
-    out = new EnzoRiemannHLLC(integrable_quantities);
+    ASSERT("EnzoRiemann::construct_riemann",
+           "The \"HLLC\" Riemann Solver can't support mhd", !mhd);
+    out = new EnzoRiemannHLLC(internal_energy);
+
   } else if (formatted == std::string("hlld")){
-    // could possibly check that MHD fields are included
-    out = new EnzoRiemannHLLD(integrable_quantities);
+    ASSERT("EnzoRiemann::construct_riemann",
+           "The \"HLLD\" Riemann Solver requires magnetic fields", mhd);
+    out = new EnzoRiemannHLLD(internal_energy);
+
   } else {
     ERROR("EnzoRiemann::construct_riemann",
 	  "The only known solvers are HLL, HLLE, HLLC, & HLLD");
-    out = NULL;  // Deals with compiler warning
   }
 
   return out;

src/Enzo/enzo_EnzoRiemann.hpp

@@ -22,15 +22,15 @@ class EnzoRiemann : public PUP::able
   /// Factory method for constructing the EnzoRiemann object. (The signature
   /// may need to be modified as additional physics get added)
   ///
-  /// @param integrable_quantities A vector of integration quantities (which
-  ///     must be listed as advected quantities in FIELD_TABLE). This is used
-  ///     to register the quantities that the Riemann Solver operates on.
   /// @param solver The name of the Riemann solver to use. Valid names include
   ///     "hll", "hlle", and "hlld"
-  static EnzoRiemann* construct_riemann
-    (std::vector<std::string> integrable_quantities, std::string solver);
+  /// @param mhd Indicates whether magnetic fields are present
+  /// @param internal_energy Indicates whether the internal energy is an
+  ///     integration quantity
+  static EnzoRiemann* construct_riemann(std::string solver, bool mhd,
+                                        bool internal_energy);
 
-  EnzoRiemann() throw()
+  EnzoRiemann() noexcept
   {}
 
   /// Virtual destructor
@@ -54,12 +54,16 @@ class EnzoRiemann : public PUP::able
   /// Computes the Riemann Fluxes for each conserved field along a given
   /// dimension, dim
   /// @param[in]     priml_map,primr_map Maps of arrays holding the left/right
-  ///     reconstructed face-centered primitives. This should be
-  ///     face-centered along `dim` (without having values on the exterior
-  ///     faces of the block) and cell-centered along the other dimensions.
-  /// @param[out]    flux_map Holds arrays where the calculated fluxes
-  ///     will be stored. The arrays should be face-centered along `dim`
-  ///     (without having values on the exterior faces of the block)
+  ///     reconstructed face-centered primitives. A list of the expected
+  ///     primitive quantities is provided by the `primitive_quantities` method.
+  ///     These should be face-centered along `dim` (without having values on
+  ///     the exterior faces of the block) and cell-centered along the other
+  ///     dimensions.
+  /// @param[out]    flux_map Holds arrays where the calculated fluxes for the
+  ///     integration quantities will be stored. Fluxes are computed for each
+  ///     integration quantities in the list provided by the
+  ///     ``integration_quantities`` method. The arrays should be face-centered
+  ///     along `dim` (without having values on the exterior faces of the block)
   /// @param[in]     dim Dimension along which to compute Riemann fluxes.
   ///     Values of 0, 1, and 2 correspond to the x, y, and z directions.
   /// @param[in]     eos Instance of the fluid's EnzoEquationOfState object
@@ -85,6 +89,15 @@ class EnzoRiemann : public PUP::able
    int stale_depth, const str_vec_t &passive_list,
    EFlt3DArray *interface_velocity) const = 0;
 
+  /// Return the names of the actively advected integration quantities for
+  /// which fluxes will be computed.
+  virtual const std::vector<std::string> integration_quantities()
+    const noexcept = 0;
+
+  /// Return the names of the (non passive scalar) primitives that are required
+  /// to compute the flux.
+  virtual const std::vector<std::string> primitive_quantities()
+    const noexcept = 0;
 };
 
 #endif /* ENZO_ENZO_RIEMANN_HPP */

src/Enzo/enzo_EnzoRiemannHLL.hpp

@@ -168,7 +168,7 @@ struct EinfeldtWavespeed
 template<class inputLUT>
 struct DavisWavespeed
 {
-  /// @class    DavisHLLEWavespeed
+  /// @class    DavisWavespeed
   /// @ingroup  Enzo
   /// @brief    [\ref Enzo] Encapsulates one of the wavespeed estimators
   ///           described by Davis 1998. It was used in Enzo's Runge-Kutta MHD
@@ -214,9 +214,9 @@ struct DavisWavespeed
 template <class WaveSpeedFunctor>
 struct HLLImpl
 {
-  /// @class    EnzoRiemannHLLE
+  /// @class    HLLImpl
   /// @ingroup  Enzo
-  /// @brief    [\ref Enzo] Encapsulates operations of HLLE approximate Riemann
+  /// @brief    [\ref Enzo] Encapsulates operations of HLL approximate Riemann
   /// Solver. 
 public: