VL+CT support gravity source terms

doc/source/param/initial.incl

@@ -266,12 +266,12 @@ Initial:cloud:perturb_standard_deviation is positive.`
 inclined_wave
 -------------
 
-The :p:`inclined_wave` Initial subgroup is used to setup a HD or
-MHD wave at an angle inclined to the simulation domain for testing
+The :p:`inclined_wave` Initial subgroup is used to setup a HD, MHD,
+or Jeans wave at an angle inclined to the simulation domain for testing
 HD/MHD integrators. If applicable, magnetic fields fields are set
 to zero when a HD wave is initialized.
 
-The procedure to addopted from
+The initialization procedure was adopted from
 `Gardiner & Stone (2008) <http://adsabs.harvard.edu/abs/2008JCoPh.227.4123G>`_ .
 Specifically, a coordinate system "x0", "x1", "x2" is defined and the
 wave is initialized to travel along "x0". The transformation between
@@ -313,9 +313,16 @@ vector potential to ensure that they are divergence-free.
 :Default: :d:`alfven`
 :Scope:   :z:`Enzo`
 
+:e:`This value specifies the type of wave to initialize. We have
+provided more details about each option down below. Note, when using an MHD
+solver with a non-MHD wave, the mangetic fields are uniformly initialized
+to zero.`
+
+Hydro Waves
+~~~~~~~~~~~
 :e:`The values used to initialize hydrodynamical linear waves are taken from
 the columns of the matrix given in equation B3 of`
-`Stone et al. (2008) <http://adsabs.harvard.edu/abs/2008ApJS..178..137S>`_
+`Stone et al. (2008) <https://adsabs.harvard.edu/abs/2008ApJS..178..137S>`_
 :e:`. Valid hydrodynamical waves include:`
 
   * ``"sound"`` :e:`A linear sound wave.`
@@ -331,8 +338,10 @@ the columns of the matrix given in equation B3 of`
     perturbations in velocity component v2 (transverse to the direction of
     bulk motion).`
 
-:e:`Each of the valid MHD waves are described in`
-`Gardiner & Stone (2008) <http://adsabs.harvard.edu/abs/2008JCoPh.227.4123G>`_
+MHD Waves
+~~~~~~~~~
+:e:`Each of the valid MHD waves are described in` `Gardiner & Stone
+(2008) <https://adsabs.harvard.edu/abs/2008JCoPh.227.4123G>`__
 :e:`. Valid MHD wave types include:`
 
   * ``"alfven"`` :e:`A linear Alfven wave with perturbations to the
@@ -343,6 +352,41 @@ the columns of the matrix given in equation B3 of`
   * ``"fast"`` :e:`A linear fast magnetosonic wave.`
   * ``"slow"`` :e:`A linear slow magnetosonic wave.`
 
+Jeans Wave
+~~~~~~~~~~
+
+:e:`To initialize a Jeans wave, set this parameter to` ``"jeans"``.
+:e:`We use equations consistent with what Athena and (earlier
+versions of) Athena++ use. In detail, we use:`
+
+.. math::
+   \rho     &= \rho_{\rm bkg}
+          \left( 1 + A \sin\left(\frac{2\pi}{\lambda} x_0 \right) \right) \\
+   \rho v_0 &= \rho_{\rm bkg} \frac{\sqrt{|\omega^2|}}{2\pi / \lambda} A
+          \begin{cases}
+                0                       & \omega^2 > 0 \\
+                \cos\left(\frac{2\pi}{\lambda} x_0 \right) & \omega^2 < 0 \\
+          \end{cases} \\
+   v_1      &= 0 \\
+   v_2      &= 0 \\
+   \rho E   &= \frac{P_{\rm bkg}}{\gamma - 1} + \gamma A \sin\left(\frac{2\pi}{\lambda} x_0\right) \\
+
+:e:`in which:`
+
+  * :math:`E` :e:`specifies the specific internal energy`
+  * :math:`\gamma` :e:`is the adiabatic index.`
+  * :math:`A` :e:`is the amplitude, specified by
+    Initial:inclined_wave:amplitude and` :math:`\lambda` :e:`is the
+    wavelength, specified by Initial:inclined_wave:lambda`
+  * :math:`\rho_{\rm bkg}=1` :e:`and` :math:`P_{\rm bkg}=1/\gamma`
+    :e:`in the appropriate code units`
+  * :math:`\omega^2=(2\pi c_{s,{\rm bkg}}/ \lambda)^2 (1 -
+    (\lambda/\lambda_J)^2)` :e:`is the dispersion relation. In this
+    equation,` :math:`c_{s,{\rm bkg}}^2=\gamma P_{\rm bkg}/\rho_{\rm
+    bkg}` :e:`and` :math:`\lambda_J = c_{s,{\rm bkg}} \sqrt{\pi / (G
+    \rho_{\rm bkg})}`. :e:`Note that the value of` :math:`G` :e:`is
+    directly set by Method:gravity:grav_const`.
+
 ----
 
 :Parameter:  :p:`Initial` : :p:`inclined_wave` : :p:`amplitude`
@@ -372,9 +416,9 @@ polarized Alfven wave (for that case, amplitude is fixed at 0.1).`
 :Default: :d:`true`
 :Scope:   :z:`Enzo`
 
-:e:`This has no effect for the circularly polarized Alfven wave. This is
-ignored for linear HD entropy waves when Initial:inclined_wave:parallel_vel
-is specified.`
+:e:`Do not specify this parameter when initializing a circularly
+polarized Alfven wave or a Jeans wave. This is ignored for linear HD
+entropy waves when Initial:inclined_wave:parallel_vel is specified.`
 
 ----
 
@@ -384,8 +428,9 @@ is specified.`
 :Default: :d:`none`
 :Scope:   Enzo
 
-:e:`This can be used to specify a background velocity along v0 for
-HD linear waves.`
+:e:`This can be used to specify a background velocity along v0 for HD
+linear waves. At present, this parameter should only be specified for
+the hydrodynamic waves.`
 
 merge_sinks_test
 ----------------

doc/source/tests/gravity.rst

@@ -2,15 +2,58 @@
 Gravity Tests
 -------------
 
-Tests for gravity methods. Runs gravity method using cg solver. Outputs mesh, phi, rho, x acceleration, y acceleration image files.
+Tests for gravity methods.
+
+Completion-Time Tests
+=====================
+
+Runs gravity method using cg solver. Outputs mesh, phi, rho, x
+acceleration, y acceleration image files. Tests succeed or fail based
+on the final completion time of the simulation.
 
 method_gravity_cg-1
-===================
+~~~~~~~~~~~~~~~~~~~
 
 Tests ``EnzoMethodGravityCg`` at P=1 in 2D.
 
 
 method_gravity_cg-8
-===================
+~~~~~~~~~~~~~~~~~~~
 
 Tests ``EnzoMethodGravityCg`` at P=8 in 2D.
+
+Analytic-Tests
+==============
+
+These tests are associated with a dedicated python-script that runs a
+test problem with a known analytic solution.
+
+stable_jeans_wave
+~~~~~~~~~~~~~~~~~
+
+The basic setup for this test is to initialize a 1D stable Jeans wave
+in a 3D box that is inclined with respect all 3 Cartesian axes. Enzo-E
+writes the initial snapshot to disk, evolves the wave over a complete
+period, and then writes the final snapshot result to disk. Then
+the python script computes the :math:`L_1` error norm computed by
+comparing values in the initial and final snapshots (the analytic
+solution expects them to be identical) and compares it with a
+precomputed value.
+
+For this test to pass, the gravity solver and hydro solver's gravity
+source terms must each be implemented correctly. At this time, this
+test has only been defined to use the VL+CT solver (without magnetic
+fields) and the cg-solver. The python script invokes the test at 2
+resolutions (16 cells per wavelength and 32 cells per wavelength).
+
+To invoke this test, invoke the following command from the root
+directory of the repository:
+
+.. code-block:: bash
+
+   python input/Gravity/run_stable_jeans_wave_test.py \
+     --launch_cmd <path-to-enzoe-binary>
+
+where ``<path-to-enzoe-binary>`` is replaced with the path to the Enzo-E
+binary that you want to test. The test prints a summary of the results to
+stdout and an exit code of of 0 indicates that all subtests passed.

input/Gravity/jeans_wave/initial_jeans.incl

@@ -0,0 +1,121 @@
+# Initial conditions for inclined stable Jeans Wave
+
+   include "input/vlct/MHD_linear_wave/domain-lw.incl"
+
+   Adapt {
+      max_level = 0;
+   }
+
+   Method {
+
+      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.
+         # The formula for Jeans wavelength is (in code units):
+         #      lambda_J = sqrt( pi * cs^2 / (G * rho)),
+         # where G is the gravitational constant. The way we set up the problem,
+         # rho and cs^2 are always equal to 1 in our code units.
+         #
+         # If we just adjust code mass, then our formula becomes:
+         #      lambda_J = sqrt(pi / G)
+         #      G = pi / lambda_J^2
+         # let's pick lambda_J = 2.0
+         grav_const = 0.25 * pi;
+      }
+   }
+
+   Initial{
+      list = ["inclined_wave"];
+
+      inclined_wave{
+         # For coordinates x0, x1, x2 - a wave is initialized along x0
+	 # To transform x, y, z to x0, x1, x2:
+	 #    1. Rotate x-y axis about the z-axis by angle -beta (rotate
+	 #       clockwise). The rotated y-axis is now the x1-axis
+	 #    2. Rotate the z-axis and the rotated x axis about the x1-axis
+	 #       by angle -alpha (rotate clock-wise).
+         alpha = 0.7297276562269663; # sin(alpha) = 2/3
+	 beta  = 1.1071487177940904; # sin(beta) =  2/sqrt(5)
+
+         lambda = 1.; #wavelength
+
+         # amplitude of the linear wave (has no effect on circularly
+	 # polarized alfven wave)
+	 amplitude = 1.e-6;
+
+	 # set the sign of the wave speed. default is true
+	 # (doesn't accomplish anything for circularly polarized Alfven waves
+         # or Jeans waves)
+	 positive_vel = true;
+
+         # the dispersion relationship is given by
+         #     omega^2 = (2*pi/lambda)^2 * cs_0^2 * (1 - (lambda / lambda_J)^2)
+         #     lambda_J = sqrt(pi * cs_0^2 / (G * rho_0))
+         # see the notes preceeding Method:gravity:grav_const for details about
+         # adjusting these parameters
+         wave_type = "jeans";
+      }
+   }
+
+   Field {
+      gamma = 1.6666666666666667;
+
+      alignment = 8;
+      ghost_depth = 4;
+      history = 1;
+      list = ["density", "velocity_x", "velocity_y", "velocity_z",
+              "total_energy", "pressure", "internal_energy",
+              "acceleration_x", "acceleration_y", "acceleration_z",
+              "density_gas", "density_particle", "density_total",
+              "density_particle_accumulate", "potential", "potential_temp",
+              "potential_copy", "X", "B", "X_copy", "B_copy", "particle_mass"];
+     padding = 0;
+   }
+
+
+
+   Solver {
+      list = ["cg"];
+      cg {
+         type = "cg";
+         iter_max = 1000;
+         res_tol  = 1e-14;
+         monitor_iter = 25;
+         solve_type = "block";
+      }
+   }
+
+   Stopping{
+      # the final stopping time is 2 * pi / omega, where is omega is angular
+      # frequecy. Omega is given by:
+      #    omega^2 = (2*pi/lambda)^2 * cs_0^2 * (1 - (lambda / lambda_J)^2)
+      # The final time is:
+      #     lambda / (cs_0 * sqrt(1 - (lambda / lambda_J)^2))
+      # in this setup: c_s = 1, lambda = 1, lambda_J = 2
+      time = 1.1547005383792517; # = 1.0 / sqrt(0.75)
+   }
+
+   Output {
+      list = ["data"];
+      data {
+         type = "data";
+	 field_list = ["density", "velocity_x", "velocity_y", "velocity_z",
+                       "total_energy"];
+         name = ["data-%03d.h5", "proc"];
+         #dir = ["method_ppm-cg-1-jeansN16_%.4f","time"];
+         schedule {
+             var = "time";
+             list = [0.0, 1.1547005383792517];
+         }
+      }
+   }

input/Gravity/jeans_wave/initial_jeans_ppm.incl

@@ -0,0 +1,14 @@
+include "input/Gravity/jeans_wave/initial_jeans.incl"
+
+   Method {
+
+      list = ["pm_deposit", "gravity", "ppm"];
+
+      ppm {
+         courant = 0.8;
+         diffusion = false;
+         dual_energy = false;
+         flattening = 0;
+         steepening = false;
+      }
+   }

input/Gravity/jeans_wave/initial_jeans_vlct.incl

@@ -0,0 +1,4 @@
+include "input/vlct/vl.incl"
+include "input/Gravity/jeans_wave/initial_jeans.incl"
+
+Method { list = ["pm_deposit", "gravity", "mhd_vlct"]; }

input/Gravity/jeans_wave/ppm_stable_jeansN16.in

@@ -0,0 +1,11 @@
+   include "input/Gravity/jeans_wave/initial_jeans_ppm.incl"
+
+   Mesh {
+      root_rank = 3; # 3D
+      root_blocks = [1,1,1];
+      root_size = [32,16,16]; # number of cells per axis
+   }
+
+   Output {
+      data { dir = ["method_ppm-cg-1-inclined-jeansN16_%.4f","time"]; }
+   }

input/Gravity/jeans_wave/ppm_stable_jeansN32.in

@@ -0,0 +1,11 @@
+   include "input/Gravity/jeans_wave/initial_jeans_ppm.incl"
+
+   Mesh {
+      root_rank = 3; # 3D
+      root_blocks = [1,1,1];
+      root_size = [64,32,32]; # number of cells per axis
+   }
+
+   Output {
+      data { dir = ["method_ppm-cg-1-inclined-jeansN32_%.4f","time"]; }
+   }

input/Gravity/jeans_wave/vlct_stable_jeansN16.in

@@ -0,0 +1,11 @@
+   include "input/Gravity/jeans_wave/initial_jeans_vlct.incl"
+
+   Mesh {
+      root_rank = 3; # 3D
+      root_blocks = [1,1,1];
+      root_size = [32,16,16]; # number of cells per axis
+   }
+
+   Output {
+      data { dir = ["method_vlct-cg-1-inclined-jeansN16_%.4f","time"]; }
+   }

input/Gravity/jeans_wave/vlct_stable_jeansN32.in

@@ -0,0 +1,11 @@
+   include "input/Gravity/jeans_wave/initial_jeans_vlct.incl"
+
+   Mesh {
+      root_rank = 3; # 3D
+      root_blocks = [1,1,1];
+      root_size = [64,32,32]; # number of cells per axis
+   }
+
+   Output {
+      data { dir = ["method_vlct-cg-1-inclined-jeansN32_%.4f","time"]; }
+   }