Fixed bug relating to drift-time for gas-deposition in EnzoMethodPmDe…

…posit

Previously the drift time for the gas density was set to the value of `alpha` (taken from `Method:pm_deposit:alpha`). As @stefanarridge pointed out in PR #89 this didn't make a lot of sense. I also ran into issues with introducing a Stable Jeans Wave regression test in PR #186 that required me to set `Method:pm_deposit:alpha` to 0.

For comparison, the drift time for the gravitating mass particles is set to `alpha * dt / cosmo_a`, in which:
- `alpha` again comes from `Method:pm_deposit:alpha`
- `dt` is the time-step during the current cycle
- `cosmo_a` is the scale factor computed at `time + alpha * dt` (where `time` is the time at the start of the current cycle)

To investigate this issue I checked the original code from enzo-dev in `Grid_DepositBaryons.C`. I concluded that enzo-dev sets the gas density's drift-timestep should get set in exactly the same way as the Particle's drift timestep. However, when using the PPM and Zeus Hydro-solvers, the drift timestep is set to 0 (this is consistent with a comment @johnwise made during his review of PR #189). Since our 2 primary hydro-solvers in Enzo-E (Ppm and VL+CT) currently handle these source terms with the same temporal order as these solvers, we now force the gas-deposition drift timestep to be 0 in `EnzoMethodPmDeposit`.

We will need to revisit this exact behavior when we eventually modify the VL+CT solver to handle self-gravity in a higher temporal-order manner.

input/Gravity/jeans_wave/initial_jeans.incl

@@ -10,14 +10,6 @@
 
       list = ["pm_deposit", "gravity"];
 
-      pm_deposit {
-         # every cycle, this method tries to "drift" the density field forward
-         # (like it's a grid of particles) by alpha*timestep.
-
-         alpha = 0.0; # this test-problem encounters issues when this is
-                      # non-zero
-     }
-
       gravity {
          solver = "cg";
          # choose mass units to adjust the value of the jeans wavelength is 0.5.

src/Enzo/enzo_EnzoMethodPmDeposit.cpp

@@ -434,7 +434,7 @@ namespace {
   /// @param[in, out] density_tot_arr The array where density gets accumulated
   /// @param[in]      field Contains the field data to use for accumulation
   /// @param[in]      dt Length of time to "drift" the density field before
-  ///     deposition.
+  ///     deposition divided by the scale factor at the deposition time
   /// @param[in]      hx_prop,hy_prop,hz_prop The width of cell along
   ///     each axis. These specify the proper lengths at the time that we
   ///     deposit the density (after any drift).
@@ -443,7 +443,7 @@ namespace {
   /// @param[in]      gx,gy,gz Specifies the number of cells in the ghost zone
   ///     for each dimensions
   void deposit_gas_(const CelloArray<enzo_float, 3>& density_tot_arr,
-                    Field& field, enzo_float dt,
+                    Field& field, enzo_float dt_div_cosmoa,
                     enzo_float hx_prop, enzo_float hy_prop, enzo_float hz_prop,
                     int mx, int my, int mz,
                     int gx, int gy, int gz){
@@ -489,7 +489,7 @@ namespace {
 
     FORTRAN_NAME(dep_grid_cic)(de, deposited_gas_density.data(), temp.data(),
                                vxf, vyf, vzf,
-                               &dt, rfield.data(), &rank,
+                               &dt_div_cosmoa, rfield.data(), &rank,
                                &hx_prop,&hy_prop,&hz_prop,
                                &mx,&my,&mz,
                                &gxi,&gyi,&gzi,
@@ -562,6 +562,7 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
       cosmology->compute_expansion_factor(&cosmo_a,&cosmo_dadt,
                                           block->time() + alpha_*block->dt());
     }
+    const double dt_div_cosmoa = alpha_ * block->dt() / cosmo_a;
 
     double hx,hy,hz;
     block->cell_width(&hx,&hy,&hz);
@@ -576,8 +577,6 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
       if (rank >= 2) inv_vol /= hy;
       if (rank >= 3) inv_vol /= hz;
 
-      const double dt_div_cosmoa = alpha_ * block->dt() / cosmo_a;
-
       deposit_particles_(density_particle_arr, block, dt_div_cosmoa, inv_vol,
                          mx, my, mz,
                          gx, gy, gz);
@@ -590,12 +589,18 @@ void EnzoMethodPmDeposit::compute ( Block * block) throw()
 
     // add mass from gas
     {
-      double gas_dt = alpha_; // probably a typo
+      // Grid_DepositBaryons.C from enzo-dev overrides the drift time for the
+      // gas density to be zero when using the PPM and Zeus solvers.
+      //
+      // The way that the Gravity source terms are currently included with
+      // EnzoMethodPpm and EnzoMethodMHDVlct is currently very similar (in
+      // terms of temporal integration), so we also set this to zero.
+      const double gas_dt_div_cosmoa = 0.0;
 
       // The use of "proper" cell-widths was carried over for consistency with
       // earlier versions of the code. Based on Grid_DepositBaryons.C from
       // enzo-dev, it seems like this may not be correct.
-      deposit_gas_(density_tot_arr, field, gas_dt,
+      deposit_gas_(density_tot_arr, field, gas_dt_div_cosmoa,
                    hx*cosmo_a, hy*cosmo_a, hz*cosmo_a,
                    mx, my, mz,
                    gx, gy, gz);