WIP: Cosmology support for VL+CT

doc/source/param/initial.incl

@@ -28,6 +28,7 @@
    * :t:`"turbulence"` :e:`Initialize fields for driving turbulence, including "driving_[xyz]" fields.`
    * :t:`"value"` :e:`Initialize fields using expressions directly from the parameter file.`
    * :t:`"vlct_bfield"` :e:`Initialize the cell-centered magnetic fields for use by the VL + CT method. For more details, see` :ref:`vlct_bfield_init`
+   * :t:`"zeldovich_pancake"` :e:`Initialize the Zeldovich Pancake test problem. For more details, see` :ref:`the section describing the tests using this initializer <zeldovich-pancake-test>`.
 
 Parameters specific to individual initializers are specified in subgroups.
 
@@ -1269,4 +1270,23 @@ integrators that don't require face-centered B-fields*
    The arguements for this parameter follow the same sets of rules as the
    parameters of Initial:value. If this parameter is not specified, but the
    values of the other components of the magnetic vector potential are, then
-   this component is assumed to be zero everywhere.`
+   this component is assumed to be zero everywhere.`
+
+
+zeldovich_pancake
+-----------------
+
+This subgroup is used to initialize a problem based on the solution to the Zel'dovich Pancake problem.
+We use the density and velocity profiles from the solution for a pressure-free fluid to initialize an ordinary gas.
+
+----
+
+.. par:parameter:: Initial:zeldovich_pancake:aligned_ax
+
+   :Summary: :s:`Specify the axis along which the Zel'dovich Pancake evolves.`
+   :Type:    :par:typefmt:`string`
+   :Default: :d:`x`
+   :Scope:   :z:`Enzo`
+
+   :e:`Allowed values are` ``"x"`` :e:`,`  ``"y"`` :e:`, or` ``"z"`` :e:`.`
+

doc/source/param/method_mhd_vlct.incl

@@ -40,6 +40,31 @@ The Method:mhd_vlct:mhd_choice determines whether the method is used as a pure h
 
 ----
 
+.. par:parameter:: Method:mhd_vlct:time_scheme
+
+   :Summary: :s:`name of the time-integration scheme to use`
+   :Type:   :par:typefmt:`string`
+   :Default: :d:`vl`
+   :Scope:     :z:`Enzo`
+
+   :e:`Name of the time integration scheme to use. The recommended
+   choice is` ``"vl"``, :e:`which corresponds to the default 2-stage
+   predictor-corrector scheme. This should generally be used with a
+   courant factor satisfying` :math:`C \leq 0.5`.
+
+   :e:`At present, the only other option is` ``"euler"``, :e:`which just
+   updates the MHD fields in a single-stage. It is` **ONLY INTENDED FOR
+   TESTING PURPOSES** :e:`and at the time of writing this documentation, it
+   has not been rigorously tested. When using this choice, make sure to pass`
+   ``"nn"`` :e:`to the` :par:param:`~Method:mhd_vlct:reconstruct_method`
+   :e:`parameter. This scheme should generally be compatible with courant
+   factors satisfying` :math:`C \leq 1`.
+
+   :e:`In the future, we may add additional options to this parameter to
+   support higher order Runge-Kutta integration schemes.`
+
+----
+
 .. par:parameter:: Method:mhd_vlct:riemann_solver
 
    :Summary: :s:`name of the Riemann solver to use`
@@ -52,6 +77,39 @@ The Method:mhd_vlct:mhd_choice determines whether the method is used as a pure h
 
 ----
 
+.. par:parameter:: Method:mhd_vlct:reconstruct_method
+
+   :Summary: :s:`name of the reconstruction method`
+   :Type:   :par:typefmt:`string`
+   :Default: :d:`plm`
+   :Scope:     :z:`Enzo`
+
+   :e:`Name of the interpolation method used to reconstruct face-centered
+   primitives for computing the fluxes. For a list of options, see`
+   :ref:`using-vlct-reconstruction`
+
+   .. note::
+
+      When :par:param:`~Method:mhd_vlct:time_scheme` is ``"vl"``, this
+      has no effect on the predictor stage (aka the half-timestep); it
+      only affects the reconstruction method in the second stage. The
+      predictor stage **MUST** use nearest-neighbor reconstruction.
+
+----
+
+.. par:parameter:: Method:mhd_vlct:theta_limiter
+
+   :Summary: :s:`controls the dissipation of certain slope limiters.`
+   :Type:   :par:typefmt:`float`
+   :Default: :d:`1.5`
+   :Scope:     :z:`Enzo`
+
+   :e:`Modifies the disipation of the slope limiter of the`
+   ``"plm"``/``"plm_enzo"`` :e:`piecewise linear reconstruction
+   algorithm.  For more details, see` :ref:`using-vlct-reconstruction`
+
+----
+
 .. par:parameter:: Method:mhd_vlct:half_dt_reconstruct_method
 
    :Summary: :s:`name of the reconstruction method to use for the full timestep`
@@ -65,6 +123,16 @@ The Method:mhd_vlct:mhd_choice determines whether the method is used as a pure h
    :e:`are used). For a list of options, see`
    :ref:`using-vlct-reconstruction`
 
+   .. warning::
+
+      This parameter is deprecated and will be removed in a future version.
+      It only carried meaning while using the original predictor-corrector
+      time integration scheme.
+
+      Furthermore, this parameter only conveyed the illusion of choice. In
+      reality, the integrator ONLY worked when this was set to ``"nn"``. For
+      that reason, this parameter is not replaced.
+
 ----
 
 .. par:parameter:: Method:mhd_vlct:full_dt_reconstruct_method
@@ -78,18 +146,13 @@ The Method:mhd_vlct:mhd_choice determines whether the method is used as a pure h
    primitives for the full timestep. For a list of options, see`
    :ref:`using-vlct-reconstruction`
 
-----
+   
 
-.. par:parameter:: Method:mhd_vlct:theta_limiter
+   .. warning::
 
-   :Summary: :s:`controls the dissipation of certain slope limiters.`
-   :Type:   :par:typefmt:`float`
-   :Default: :d:`1.5`
-   :Scope:     :z:`Enzo`
-
-   :e:`Modifies the disipation of the slope limiter of the`
-   ``"plm"``/``"plm_enzo"`` :e:`piecewise linear reconstruction
-   algorithm.  For more details, see` :ref:`using-vlct-reconstruction`
+      This parameter is deprecated and will be removed in a future version.
+      The :par:param:`~Method:mhd_vlct:reconstruct_method` parameter is a
+      direct replacement.
 
 ----
 

doc/source/tests/existing_tests.rst

@@ -72,3 +72,4 @@ The answer test suite currently covers the following simulations:
 
 .. toctree::
    grackle-pytest
+   zeldovich-pancake

doc/source/tests/zeldovich-pancake.rst

@@ -0,0 +1 @@
+.. include:: ../../../input/Cosmology/ZeldovichPancake/description.rst

doc/source/user/problem_initial.rst

@@ -121,3 +121,8 @@ a 3D array of Sedov blast waves.
    partially initialized ``"total_energy"`` fields with the specific
    magnetic energy computed from the newly computed cell-centered
    bfields and pre-initialized ``"density"`` fields.
+
+``"zeldovich_pancake"``
+
+   Initialize the Zeldovich-Pancake test problem.
+   For more information, see :ref:`the section describing the tests using this initializer <zeldovich-pancake-test>`.

doc/source/user/problem_method.rst

@@ -306,6 +306,7 @@ parameters
        the minimum across all 3 dimensions), at the highest refinement
        level.`
 
+.. _vlct_overview:
 
 ``"mhd_vlct"`` method
 =====================
@@ -356,23 +357,23 @@ parameter) should be less than 0.5.
      - Type
      - Default
      - Description
-   * - ``"mhd_choice"``
+   * - :par:param:`~Method:mhd_vlct:mhd_choice`
      - `string`
      - `none`
      - `Specifies handling of magnetic fields (or lack thereof)`
-   * - ``"riemann_solver"``
+   * - :par:param:`~Method:mhd_vlct:time_scheme`
+     - `string`
+     - `vl`
+     - `Specifies time integration scheme (CURRENTLY JUST FOR TESTING)`
+   * - :par:param:`~Method:mhd_vlct:riemann_solver`
      - `string`
      - `hlld`
      - `name of the Riemann Solver to use`
-   * - ``"half_dt_reconstruct_method"``
-     - `string`
-     - `nn`
-     - `Reconstruction method for half timestep`
-   * - ``"full_dt_reconstruct_method"``
+   * - :par:param:`~Method:mhd_vlct:reconstruct_method`
      - `string`
      - `plm`
-     - `Reconstruction method for full timestep`
-   * - ``"theta_limiter"``
+     - `Reconstruction method`
+   * - :par:param:`~Method:mhd_vlct:theta_limiter`
      - `float`
      - `1.5`
      - `controls dissipation of the "plm"/"plm_enzo" reconstruction
@@ -487,13 +488,13 @@ the stale depth increases. We define the amount by which it
 increases as the "staling rate" (which depends on the choice
 of interpolation method).
 
-For a unigrid simulation, the number of required ghost zones
-is given by the sum of the staling rates for each selected
-reconstruction method.
+For a unigrid simulation, the number of ghost zones is ``1 + staling_rate``
+where ``staling_rate`` depends on the chosen reconstruction method.
 
-We provide the names used to specify each available method in
-the input file, the associated staling depth, and a brief
-description.
+For each available reconstruction method, the following table provides
+the name that must be passed to
+:par:param:`~Method:mhd_vlct:reconstruct_method` in the input file,
+the associated staling depth, and a brief description.
 
 .. list-table:: Available ``mhd_vlct`` reconstructors (and slope
 		limiters)
@@ -531,21 +532,32 @@ description.
 
 We provide a few notes about the choice of interpolator for this algorithm:
 
-   * The recommended choices of reconstruction algorithms are ``"nn"`` for the
-     half-timestep and then piecewise-linear reconstruction for the
-     full-timestep (most test problems have been run using ``plm`` with
-     ``theta_limiter=2``, matching the integrator description in
-     `Stone & Gardiner 2009
-     <https://adsabs.harvard.edu/abs/2009NewA...14..139S>`_ ). Using ``"nn"``
-     both times also works, however tests show that errors arise when
-     piecewise linear reconstruction is used both times.
+   * The recommended choices of reconstruction algorithm is
+     piecewise-linear reconstruction (most test problems have been run
+     using ``plm`` with ``theta_limiter=2``, matching the integrator
+     description in `Stone & Gardiner 2009
+     <https://adsabs.harvard.edu/abs/2009NewA...14..139S>`_ ).
    * It is supposed to be possible to reconstruct the characteristic quantities
      for this method or to use higher order reconstruction in place of ``"plm"``
    * Reconstruction is always performed on the cell-centered magnetic fields.
      After reconstructing values along a given axis, the values of the
      reconstructed magnetic field component for that axis are replaced by the
      face-centered magnetic field values.
 
+.. note::
+
+   Under the hood, the default predictor-corrector temporal integration scheme
+   **always** uses ``"nn"`` reconstruction to compute fluxes during the
+   prediction stage (aka the partial timestep); the scheme breaks for any other
+   choices.
+
+   This is why the number of required ghost zones is ``1 + staling_rate``:
+      * the first term is the staling depth of the ``"nn"`` reconstructor
+      * the second term is the staling depth of the reconstructor specified by
+        the :par:param:`~Method:mhd_vlct:reconstruct_method` parameter; this
+        is used for computing fluxes for the second stage (for the full
+        timestep). 
+
 .. _using-vlct-riemann-solver:
 
 riemann solvers
@@ -657,7 +669,7 @@ Required Fields
 ---------------
 
 The evolved fields are ``"photon_density_i"``, ``"flux_x_i"``, ``"flux_y_i"``, and ``"flux_z_i"``, 
-where photon densities and fluxes are in units of :math:`[L]^{-3}` :e:`and` :math:`[L]^2[T]^{-1}`, 
+where photon densities and fluxes are in units of :math:`[L]^{-3}` and :math:`[L]^2[T]^{-1}`, 
 respectively, and :math:`i` denotes the group number. A set of fields must be defined for each group, and
 fluxes must be defined in the x, y, and z directions regardless of the dimensionality of the simulation.
 For technical reasons, an additional field called ``"photon_density_i_deposit"`` must be defined for each group.
@@ -855,8 +867,6 @@ This method currently ignores all of the floor parameters that are set in the ``
 
 PPML ideal MHD solver
 
-.. _vlct_overview:
-
 ``"sink_maker"`` method
 =======================
 

input/Checkpoint/checkpoint_vlct.in

@@ -24,8 +24,8 @@
          courant = 0.4;
          mhd_choice = "no_bfield";
          riemann_solver = "hllc";
-         half_dt_reconstruct_method = "nn";
-         full_dt_reconstruct_method = "plm";
+         time_scheme = "vl";
+         reconstruct_method = "plm";
          theta_limiter = 2.;
          density_floor = 1.e-30;
          pressure_floor = 1.e-30;

input/Cosmology/ZeldovichPancake/description.rst

@@ -0,0 +1,65 @@
+.. _zeldovich-pancake-test:
+
+------------------
+Zel'dovich Pancake
+------------------
+
+The Zel'dovich Pancake `(Zel'dovich 1970) <https://ui.adsabs.harvard.edu/abs/1970A%2526A.....5...84Z>`__ test problem is a standard 1D test of cosmology.
+
+* In his original work, Zel'dovich applied linear perturbation theory to solve the evolution of a pressure-free fluid (i.e. dark matter) in an isotropic universe, when a sinusoidal perturbation.
+  He arrived at an analytic solution that is exact during the linear phase, up until the point where the caustic form (this is when the pancake forms).
+  You can think of the formation of the caustic as a 1D analog for collapse into filaments and clusters.
+
+* We primarily consider the description of the solution presented in `Anninos et al. (1994) <https://adsabs.harvard.edu/abs/1994ApJ...436...11A>`__.
+  The density and velocity profiles of a pressure-free fluid at some arbitrary redshift :math:`z_0`, that occurs before :math:`z_c` (the redshift of collapse or caustic formation) are given by
+
+  .. math::
+     :nowrap:
+
+     \begin{eqnarray}
+       \rho(x_l) &=& \rho_{\rm bkg} \left[1 - \frac{1+z_c}{1+z_0}\cos(k x_l)\right]^{-1}, \\
+       v(x_l) &=& -H_0 \frac{1+z_c}{\sqrt{1+z_0}}\, \frac{\sin(k x_l)}{k}.
+     \end{eqnarray}
+
+  In these equations :math:`H_0` is the Hubble constant's value at :math:`z=0` and :math:`k=2\pi/\lambda` is the wave-number of the perturbation.
+  These equations are functions of :math:`x_l`, the Lagrangian position of the fluid elements or the positions of the fluid elements at the redshift when the perturbation is imposed.
+  When you consider the position on the grid at :math:`z_0`, you are considering a Eulerian coordinate, :math:`x_e`.
+  The conversion between these two quantities is:
+
+  .. math::
+
+     x_e = x_l - \frac{1 + z_c}{1+z_0} \frac{sin(k x_l)}{k}.
+
+In the test, we start a simulation at some redshift :math:`z_{\rm 0}`, and initialize an ordinary ideal gas such that it obeys the above density and velocity profiles (of the pressure-free fluid).
+The internal energy of this gas is initialized so that entropy is constant everywhere.
+
+This problem is useful because it checks the minimal set of physics required for a cosmological simulation (self-gravity, hydrodynamics, expansion).
+It also serves as a test of the dual energy formalism.
+
+The actual test we run is described as follows:
+
+* We consider a flat universe with :math:`\Omega_b=1` with :math:`h=0.5` in a narrow box that extends over :math:`64\, {\rm Mpc}\, h^{-1}`.
+  We initialize the simulation at :math:`z_0 = 20`, with :math:`z_c = 1`, :math:`\rho_0 = \rho_c` (the critical density), and a background temperature of :math:`100 {\rm K}`.
+
+* While this test technically is just 1D, the VL integrator can only run in 3D.
+  For that reason, we consider a long narrow box that is thin along the other axes (the side of each cell is constant).
+  We effectively run multiple identical versions of this problem in parallel.
+
+**Note on timing:** The ppm integrator will run this test in fewer cycles because of differences in the minimum courant factor.
+
+
+.. COMMENT BLOCK
+
+   We may wish to introduce other variants of this problem. For example we could modify the conditions to introduce Bfields.
+
+     * This is similar to how `Li et al. (2008) <https://adsabs.harvard.edu/abs/2008ApJS..174....1L>`__ and `Collins et al. (2010) <https://adsabs.harvard.edu/abs/2010ApJS..186..308C>`__ extend the conditions from  `Ryu et al. (1993) <https://adsabs.harvard.edu/abs/1993ApJ...414....1R>`__
+
+     * I believe the initializer from Enzo shows us how to do this
+
+   We may also want to run a version of the problem with AMR, and a version with drift velocity.
+
+   May also want to try a dual-pancake version... (for a more
+
+
+
+

input/Cosmology/ZeldovichPancake/ppm-HD-x-aligned.in

@@ -0,0 +1,34 @@
+# Problem: x-axis aligned Zeldovich Pancake problem for ppm integrator
+# Author: Matthew Abruzzo (matthewabruzzo@gmail.com)
+
+   include "input/Cosmology/ZeldovichPancake/x-aligned-HD.incl"
+
+   Physics {
+      list += ["fluid_props"];
+      fluid_props {
+         dual_energy {
+            type = "bryan95";
+            eta = [0.001, 0.1];
+         }
+      }
+   }
+
+   Method {
+      list = ["cosmology",
+              "pm_deposit",
+              "gravity",
+              "ppm",
+              "comoving_expansion"];
+      ppm {
+         courant = 0.8;
+         diffusion = false;
+         flattening = 0;
+         steepening = false;
+      }
+   }
+
+   Output {
+      data {
+         dir = ["method-ppm-HD_x-aligned-pancake_%.3e","time"];
+      }
+   }