Introduce support for using cosmology with VL+CT solver (without bfields)

Strictly speaking, this only required 2 minor tweaks to EnzoMethodMHDVlct:
  - modify the compatability checking of EnzoMethodMHDVlct::post_init_checks_ related to using the VL+CT solver with cosmology.
  - modify all usage of the cell widths in EnzoMethodMHDVlct::compute (which would have defaulted to the cell width in the comoving frame) to instead use the proper cell width.

As part of this commit, occurences of the variables named `cell_width` in the relevant code (i.e. paths executed by EnzoMethodMHDVlct::compute) were renamed so that they are called `proper_cell_width`. The docstrings of these functions were also updated accordingly.

Author: mabruzzo

src/Enzo/enzo_EnzoIntegrationQuanUpdate.cpp

@@ -101,12 +101,12 @@ void EnzoIntegrationQuanUpdate::clear_dUcons_map
 //----------------------------------------------------------------------
 
 void EnzoIntegrationQuanUpdate::accumulate_flux_component
-(int dim, double dt, enzo_float cell_width, const EnzoEFltArrayMap &flux_map,
- EnzoEFltArrayMap &dUcons_map, int stale_depth,
- const str_vec_t &passive_list) const noexcept
+(int dim, double dt, enzo_float proper_cell_width,
+ const EnzoEFltArrayMap &flux_map, EnzoEFltArrayMap &dUcons_map,
+ int stale_depth, const str_vec_t &passive_list) const noexcept
 {
   EnzoPermutedCoordinates coord(dim);
-  enzo_float dtdx_i = dt/cell_width;
+  enzo_float dtdx_i = dt/proper_cell_width;
 
   auto accumulate = [dtdx_i,stale_depth,coord,
                      &flux_map,&dUcons_map](const std::string& key)

src/Enzo/enzo_EnzoIntegrationQuanUpdate.hpp

@@ -66,7 +66,8 @@ class EnzoIntegrationQuanUpdate : public PUP::able
   /// @param[in]     dim The dimension along which to compute the flux
   ///     divergence.
   /// @param[in]     dt The current timestep.
-  /// @param[in]     cell_width The cell width along dimension `dim`.
+  /// @param[in]     proper_cell_width The cell width along dimension `dim`.
+  ///     This should always be a proper distance (even in cosmological sims).
   /// @param[in]     flux_map Map of arrays holding the fluxes computed for
   ///     the current timestep. The values of these fields should be stored
   ///     on the cell faces along the `dim` dimension.
@@ -77,7 +78,12 @@ class EnzoIntegrationQuanUpdate : public PUP::able
   ///     call. This should match the stale depth at the time the fluxes were
   ///     computed.
   /// @param[in]     passive_list A list of keys for passive scalars.
-  void accumulate_flux_component(int dim, double dt, enzo_float cell_width,
+  ///
+  /// @note
+  /// By always using the cell width measured in the proper frame, this is
+  /// compatible with both cosmological and non-cosmological simulations.
+  void accumulate_flux_component(int dim, double dt,
+                                 enzo_float proper_cell_width,
                                  const EnzoEFltArrayMap &flux_map,
                                  EnzoEFltArrayMap &dUcons_map, int stale_depth,
                                  const str_vec_t &passive_list) const noexcept;

src/Enzo/enzo_EnzoMethodMHDVlct.cpp

@@ -246,8 +246,8 @@ EnzoVlctScratchSpace* EnzoMethodMHDVlct::get_scratch_ptr_
 //----------------------------------------------------------------------
 
 void EnzoMethodMHDVlct::save_fluxes_for_corrections_
-(Block * block, const EnzoEFltArrayMap &flux_map, int dim, double cell_width,
- double dt) const noexcept
+(Block * block, const EnzoEFltArrayMap &flux_map, int dim,
+ double proper_cell_width, double dt) const noexcept
 {
 
   Field field = block->data()->field();
@@ -259,7 +259,7 @@ void EnzoMethodMHDVlct::save_fluxes_for_corrections_
   int gx, gy, gz;
   field.ghost_depth(density_field_id, &gx, &gy, &gz);
 
-  double dt_dxi = dt/cell_width;
+  double dt_dxi = dt/proper_cell_width;
 
   FluxData * flux_data = block->data()->flux_data();
   const int nf = flux_data->num_fields();
@@ -415,9 +415,24 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw()
       bfield_method_->register_target_block(block);
     }
 
-    const enzo_float* const cell_widths = enzo::block(block)->CellWidth;
-
-    double dt = block->dt();
+    const double dt = block->dt();
+
+    // compute the proper_cell_widths (independent of whether we're using
+    // comoving coordinates or not). For a pure hydro simulation, this is all
+    // we need to do to be compatible with cosmology. Cosmology is not
+    // currently supported with magnetic fields (an error is raised elsewhere).
+    // - In the future we might want to compute different values of cosmo_a for
+    //   the partial and the full timestep
+    double cosmo_a = 1.0;
+    if (enzo::cosmology() != nullptr) {
+      double cosmo_dadt = 0.0;
+      enzo::cosmology()->compute_expansion_factor
+        (&cosmo_a, &cosmo_dadt, block->time() + 0.5*dt);
+    }
+    const std::array<enzo_float,3> proper_cell_widths =
+      {enzo::block(block)->CellWidth[0] * cosmo_a,
+       enzo::block(block)->CellWidth[1] * cosmo_a,
+       enzo::block(block)->CellWidth[2] * cosmo_a};
 
     // stale_depth indicates the number of field entries from the outermost
     // field value that the region including "stale" values (need to be
@@ -426,7 +441,7 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw()
 
     // repeat the following loop twice (for half time-step and full time-step)
     for (int i=0;i<2;i++){
-      double cur_dt = (i == 0) ? dt/2. : dt;
+      const double cur_dt = (i == 0) ? dt/2. : dt;
       EnzoEFltArrayMap& cur_integration_map =
         (i == 0) ? external_integration_map : temp_integration_map;
       EnzoEFltArrayMap& out_integration_map =
@@ -480,7 +495,7 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw()
           interface_vel_arr_ptr = nullptr;
         }
 
-        compute_flux_(dim, cur_dt, cell_widths[dim], primitive_map,
+        compute_flux_(dim, cur_dt, proper_cell_widths[dim], primitive_map,
                       pl_map, pr_map, *(flux_maps[dim]), dUcons_map,
                       interface_vel_arr_ptr, *reconstructor, bfield_method_,
                       stale_depth, passive_list);
@@ -492,11 +507,11 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw()
         //   different if using SMR/AMR. This means that the internal energy
         //   source terms won't be fully self-consistent along the edges. This
         //   same effect is also present in the Ppm Solver
-        save_fluxes_for_corrections_(block, xflux_map, 0, cell_widths[0],
+        save_fluxes_for_corrections_(block, xflux_map, 0, proper_cell_widths[0],
 				     cur_dt);
-        save_fluxes_for_corrections_(block, yflux_map, 1, cell_widths[1],
+        save_fluxes_for_corrections_(block, yflux_map, 1, proper_cell_widths[1],
 				     cur_dt);
-        save_fluxes_for_corrections_(block, zflux_map, 2, cell_widths[2],
+        save_fluxes_for_corrections_(block, zflux_map, 2, proper_cell_widths[2],
 				     cur_dt);
       }
 
@@ -540,12 +555,18 @@ void EnzoMethodMHDVlct::compute ( Block * block) throw()
 
 void EnzoMethodMHDVlct::post_init_checks_() const noexcept
 {
-  ASSERT("EnzoMethodMHDVlct::post_init_checks_",
-         "This solver isn't currently compatible with cosmological sims",
-         enzo::cosmology() == nullptr);
-
   const Problem* problem = enzo::problem();
 
+  if (enzo::cosmology() != nullptr){
+    ASSERT("EnzoMethodMHDVlct::post_init_checks_",
+           ("This solver isn't currently compatible with cosmological sims "
+            "when using bfields"), mhd_choice_ == bfield_choice::no_bfield);
+    ASSERT("EnzoMethodMHDVlct::post_init_checks_",
+           "when using cosmology, this method MUST precede the "
+           "comoving_expansion method",
+           problem->method_precedes("mhd_vlct","comoving_expansion"));
+  }
+
   // problems would arise relating to particle-mesh deposition (relating to
   // particle drift before deposition) and in cosmological simulations if the
   // the VL+CT method were to precede the gravity method.
@@ -574,7 +595,7 @@ void EnzoMethodMHDVlct::post_init_checks_() const noexcept
 //----------------------------------------------------------------------
 
 void EnzoMethodMHDVlct::compute_flux_
-(const int dim, const double cur_dt, const enzo_float cell_width,
+(const int dim, const double cur_dt, const enzo_float proper_cell_width,
  EnzoEFltArrayMap &primitive_map,
  EnzoEFltArrayMap &priml_map, EnzoEFltArrayMap &primr_map,
  EnzoEFltArrayMap &flux_map, EnzoEFltArrayMap &dUcons_map,
@@ -609,7 +630,8 @@ void EnzoMethodMHDVlct::compute_flux_
   // from the fluxes and use these values to update dUcons_map (which is used
   // to accumulate the total change in these quantities over the current
   // [partial] timestep)
-  integration_quan_updater_->accumulate_flux_component(dim, cur_dt, cell_width,
+  integration_quan_updater_->accumulate_flux_component(dim, cur_dt,
+                                                       proper_cell_width,
                                                        flux_map, dUcons_map,
                                                        cur_stale_depth,
                                                        passive_list);
@@ -618,7 +640,7 @@ void EnzoMethodMHDVlct::compute_flux_
   // energy source term for this dim (and update dUcons_map).
   if (eos_->uses_dual_energy_formalism()){
     EnzoSourceInternalEnergy eint_src;
-    eint_src.calculate_source(dim, cur_dt, cell_width, primitive_map,
+    eint_src.calculate_source(dim, cur_dt, proper_cell_width, primitive_map,
                               dUcons_map, *interface_velocity_arr_ptr, eos_,
                               cur_stale_depth);
   }

src/Enzo/enzo_EnzoMethodMHDVlct.hpp

@@ -33,6 +33,9 @@
 ///        - bfield          (both integration and primitive)
 ///      When using the dual energy formalism there is also:
 ///        - internal_energy (just integration)
+///    - When using comoving coordinates, the integration and primitives are
+///      all defined in the comoving frame. Source terms related to comoving
+///      expansion are all handled externally by `EnzoMethodComovingExpansion`
 ///
 ///    EnzoEFltArrayMap Objects
 ///    ------------------------
@@ -172,13 +175,16 @@ class EnzoMethodMHDVlct : public Method {
   ///
   /// If using the dual energy formalism, this also computes a part of the
   /// internal energy density source term,
-  ///    `dt * pressure * (dvx/dx + dvy/dy + dvz/dz)`
-  /// (`vx`, `vy`, `vz` are velocity components & scale factor dependence is
-  /// omitted), and adds it to the 'internal_energy' entry in `dUcons_map`.
-  /// More specifically, it handles the dimensionally split part of the term
-  /// involving the derivative along `dim`. The velocity component along `dim`
-  /// at the cell-interfaces (estimated by the Riemann Solver) to compute the
-  /// derivatives.
+  /// `dt * pressure * (dvx/dx + dvy/dy + dvz/dz)`. In detail, this computes
+  ///   - `-dt * pressure * dvx/dx`, when `dim == 0`
+  ///   - `-dt * pressure * dvy/dy`, when `dim == 1`
+  ///   - `-dt * pressure * dvz/dz`, when `dim == 2`
+  /// and adds the results to the 'internal_energy' entry in `dUcons_map`. This
+  /// estimates the velocity derivative using the velocity component along
+  /// `dim` at the cell-interfaces (estimated by the Riemann Solver) and the
+  /// `proper_cell_width`. This is valid for both cosmological and
+  /// non-cosmological simulations. In cosmological simulations, the other
+  /// cosmological expansion term is handled in a separate `Method` object.
   ///
   /// This function should NOT be modified to directly compute any other source
   /// terms unless they similarly have dependence on dimensional quantites
@@ -188,7 +194,8 @@ class EnzoMethodMHDVlct : public Method {
   /// @param[in]     dim Dimension along which to compute fluxes. Values of 0,
   ///     1, and 2 correspond to the x, y, and z directions, respectively.
   /// @param[in]     dt The current timestep.
-  /// @param[in]     cell_width The cell width along dimension `dim`.
+  /// @param[in]     proper_cell_width The cell width along dimension `dim`.
+  ///     This should always be a proper distance (even in cosmological sims).
   /// @param[in]     primitive_map Map of arrays holding cell-centered
   ///     primitive quantities that are to be reconstructed (This includes
   ///     specific passive scalars).
@@ -219,7 +226,7 @@ class EnzoMethodMHDVlct : public Method {
   ///     performing reconstruction)
   /// @param[in]     passive_list A list of keys for passively advected scalars.
   void compute_flux_
-  (const int dim, const double cur_dt, const enzo_float cell_width,
+  (const int dim, const double cur_dt, const enzo_float proper_cell_width,
    EnzoEFltArrayMap &primitive_map,
    EnzoEFltArrayMap &priml_map, EnzoEFltArrayMap &primr_map,
    EnzoEFltArrayMap &flux_map, EnzoEFltArrayMap &dUcons_map,
@@ -283,11 +290,12 @@ class EnzoMethodMHDVlct : public Method {
   /// @param[in]  dim indicates the dimension that the fluxes in `flux_map`
   ///     were computed along. Values of 0, 1, and 2 correspond to the x, y,
   ///     and z directions, respectively.
-  /// @param[in]  cell_width is the width of a cell along the dimension `dim`
+  /// @param[in]  proper_cell_width The cell width along dimension `dim`. This
+  ///     should always be a proper distance (even in cosmological sims).
   /// @param[in]  dt is the value of the current timestep
   void save_fluxes_for_corrections_
-  (Block * block, const EnzoEFltArrayMap &flux_map, int dim, double cell_width,
-   double dt) const noexcept;
+  (Block * block, const EnzoEFltArrayMap &flux_map, int dim,
+   double proper_cell_width, double dt) const noexcept;
 
   /// Returns a pointer to the scratch space struct. If the scratch space has
   /// not already been allocated, it will be allocated now.

src/Enzo/enzo_EnzoSourceInternalEnergy.cpp

@@ -11,7 +11,7 @@
 //----------------------------------------------------------------------
 
 void EnzoSourceInternalEnergy::calculate_source
-(const int dim, const double dt, const enzo_float cell_width,
+(const int dim, const double dt, const enzo_float proper_cell_width,
  const EnzoEFltArrayMap &prim_map, EnzoEFltArrayMap &dUcons_map,
  const CelloArray<const enzo_float,3> &interface_velocity,
  const EnzoEquationOfState *eos,
@@ -26,7 +26,7 @@ void EnzoSourceInternalEnergy::calculate_source
 	 "The EOS can't be barotropic and use the dual energy formalism.",
 	 !(eos->is_barotropic()) );
 
-  const enzo_float dtdx = dt/cell_width;
+  const enzo_float dtdx = dt/proper_cell_width;
 
   EnzoPermutedCoordinates coord(dim);
 
@@ -37,9 +37,10 @@ void EnzoSourceInternalEnergy::calculate_source
   //
   // Note  - deint_dens holds the accumulated change in the internal energy
   //         density over the current time-step
-  //       - we compute pressure from the cell-centered density ane specific
-  //         internal from the start of the time-step (adopting convention from
-  //         flux_hll.F - we also could just use the precomputed field)
+  //       - we compute pressure from the cell-centered density and specific
+  //         internal energy that are taken from the start of the time-step
+  //         (if we adopted the convention from flux_hll.F, we could
+  //         alternatively use the precomputed pressure field)
 
   CSlice full_ax(nullptr, nullptr);
 

src/Enzo/enzo_EnzoSourceInternalEnergy.hpp

@@ -34,19 +34,32 @@ class EnzoSourceInternalEnergy
   /// @brief    [\ref Enzo] Encapsulates the calculation of the internal energy
   ///           source term along a given dimension.
   ///
-  /// The source term in question for the internal energy density is formally
-  /// defined as follows. We define the jth component of the interface velocity
-  /// as vbar and grid width along this dimension as dx. For all cell
-  /// locations (i,j,k) there is a source term, over a timestep dt, of:
-  ///     -1 * dt* pressure_j * (vbar_{j+1/2} - vbar_{j+1/2}) / (a * dx_j)
-  /// a corresponding term is added for each dimension (although this
-  /// implementation only handles a single dimension at a time).
+  /// For convenience, we split the encapsulated source term up by into 3
+  /// pieces; one for each dimension (this class needs to be called separately
+  /// for each dimension). We define the portion of the source term for the ith
+  /// dimension that get's applied at every cell locations as:
+  /// @code{.unparsed}
+  ///    -1 * dt * pressure[i] * (vbar[i+1/2] - vbar[i-1/2]) / (a * cellwidth_i)
+  /// @endcode
+  /// In the above equation:
+  ///   * the bracketted variables indicate the relative positions of variables
+  ///     on the grid along the ith dimension (values of `i` map to a
+  ///     cell-centered position). You can assume that the omitted indices
+  ///     along the other dimensions are always cell-centered.
+  ///   * `vbar` is the ith component of the velocity at the cell interface
+  ///   * `a` is the current scale-factor.
+  ///   * `cellwidth_i` is the comoving cell width along the ith dimension. In
+  ///     practice, this calculation always expects to be passed the proper
+  ///     cell-width, which is always equal to `a * cellwidth_i`.
   ///
-  /// As in Enzo's PPM method, the cell-centerred pressure is computed from the
+  /// This source-term is compatible with both non-cosmological and
+  /// cosmological simulations. In cosmological calculations the internal
+  /// energy density has another source term, but that get's handled separately
+  /// by `EnzoMethodComovingExpansion`.
+  ///
+  /// As in Enzo's PPM method, the cell-centered pressure is computed from the
   /// internal energy (although the internal energy and total energy should
   /// have been synchronized at the start of the timestep)
-  ///
-  /// The current implementation ignores the scale factor.
 
 public:
 
@@ -58,6 +71,8 @@ class EnzoSourceInternalEnergy
   ///     internal Energy source term. Values of 0, 1, and 2 correspond to the
   ///     the x, y, and z directions, respectively.
   /// @param[in]  dt The time time-step overwhich to apply the source term
+  /// @param[in]  proper_cell_width The cell width along dimension `dim`. This
+  ///     should always be a proper distance (even in cosmological sims).
   /// @param[in]  prim_map Map holding the values of the cell-centered
   ///     primitives from the start of the (partial) timestep (but after the
   ///     synchronization of the internal energy with the total energy).
@@ -80,7 +95,7 @@ class EnzoSourceInternalEnergy
   ///     reconstructor's delayed_staling_rate should be applied at some
   ///     time after this function call.
   void calculate_source(const int dim, const double dt,
-			const enzo_float cell_width,
+			const enzo_float proper_cell_width,
                         const EnzoEFltArrayMap &prim_map,
                         EnzoEFltArrayMap &dUcons_map,
                         const CelloArray<const enzo_float,3> &interface_velocity,