Introduce an initializer for Zeldovich Pancake test problem.

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

@@ -0,0 +1,92 @@
+   include "input/vlct/vl.incl"
+
+   Adapt {
+      max_level = 0;
+   }
+
+   Method {
+
+      list = ["cosmology", "pm_deposit", "gravity", "mhd_vlct", "comoving_expansion"];
+
+      gravity {
+         solver = "cg";
+         grav_const = 1.0;
+      }
+   }
+
+   Initial{ list = ["zeldovich_pancake"]; }
+
+   Boundary { type = "periodic"; }
+
+   Domain {
+      lower = [ 0.0, 0.0, 0.0 ];
+      upper = [ 1.0, 0.0625, 0.0625 ];
+   }
+
+   Mesh {
+      root_rank = 3;
+      root_size = [ 64, 4, 4 ];
+      root_blocks = [1, 1, 1];
+   }
+
+   Physics {
+      list = [ "cosmology" ];
+      cosmology {
+         omega_baryon_now    = 1.0;
+         omega_cdm_now       = 0.0;
+         omega_matter_now    = 1.0;
+         omega_lambda_now    = 0.0;
+         hubble_constant_now = 0.5;
+         max_expansion_rate  = 0.01;
+         initial_redshift    = 20.0;
+         final_redshift      = 0.0;
+	 comoving_box_size   = 64.0; # specifies the code length unit in terms
+                                     # of comoving Mpc / h
+      }
+   }
+
+   Field {
+      gamma = 1.6666666666666667;
+
+      alignment = 8;
+      ghost_depth = 4;
+      history = 1;
+      list += ["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{
+      #close to redshift= 1
+      time = 2.782174164581e+01;
+   }
+
+   Output {
+      list = ["data"];
+      data {
+         type = "data";
+	 field_list = ["density", "velocity_x", "velocity_y", "velocity_z",
+                       "total_energy"];
+         name = ["data-%03d.h5", "proc"];
+         dir = ["method-vlct-x-aligned-pancake_%04d","cycle"];
+         schedule {
+             var = "cycle";
+             list = [0, 1, 2, 3, 4, 5, 6, 7, 8]; # 237];
+         }
+      }
+   }

src/Enzo/_enzo.hpp

@@ -201,6 +201,7 @@ extern "C" {
 #include "enzo_EnzoInitialIsolatedGalaxy.hpp"
 #include "enzo_EnzoInitialBurkertBodenheimer.hpp"
 #include "enzo_EnzoInitialMergeSinksTest.hpp"
+#include "enzo_EnzoInitialZeldovichPancake.hpp"
 
 #include "enzo_EnzoRefineShock.hpp"
 #include "enzo_EnzoRefineParticleMass.hpp"

src/Enzo/enzo.ci

@@ -67,6 +67,7 @@ module enzo {
   PUPable EnzoInitialTurbulence;
   PUPable EnzoInitialIsolatedGalaxy;
   PUPable EnzoInitialMergeSinksTest;
+  PUPable EnzoInitialZeldovichPancake;
 
   PUPable IoEnzoBlock;
 

src/Enzo/enzo_EnzoInitialCosmology.cpp

@@ -26,7 +26,6 @@ EnzoInitialCosmology::EnzoInitialCosmology
 void EnzoInitialCosmology::enforce_block
 ( Block * block, const Hierarchy * hierarchy ) throw()
 {
-  EnzoUnits * units = enzo::units();
 
   // "If temperature is left unset, set it assuming that T=550 K at z=200"
   EnzoPhysicsCosmology * cosmology = enzo::cosmology();
@@ -36,16 +35,33 @@ void EnzoInitialCosmology::enforce_block
       pow((1.0 + cosmology->initial_redshift())/(1.0 + 200.00), 2.0);
   }
 
+  EnzoInitialCosmology::init_cosmology(block, temperature_, gamma_);
+
+  block->initial_done();
+
+}
+
+//----------------------------------------------------------------------
+
+void EnzoInitialCosmology::init_cosmology
+(Block * block, double temperature, double gamma) throw()
+{
+  EnzoPhysicsCosmology * cosmology = enzo::cosmology();
+
   // Set initial time based on initial redshift
   double r0 = cosmology->initial_redshift();
   double t0 = cosmology->time_from_redshift(r0);
   block->set_time(t0);
   enzo::simulation()->set_time(t0);
-
+
+  if (temperature < 0){ return; }
+
+  EnzoUnits * units = enzo::units();
+  // initialize energy fields
   const double default_mu = 0.6;
 
-  const double internal_energy = temperature_/units->temperature()
-    /default_mu/(gamma_-1.0);
+  const double internal_energy = temperature/units->temperature()
+    /default_mu/(gamma-1.0);
 
   Field field = block->data()->field();
 
@@ -70,7 +86,4 @@ void EnzoInitialCosmology::enforce_block
       }
     }
   }
-
-  block->initial_done();
-
 }

src/Enzo/enzo_EnzoInitialCosmology.hpp

@@ -50,6 +50,15 @@ class EnzoInitialCosmology : public Initial {
   virtual void enforce_block
   ( Block * block, const Hierarchy * hierarchy ) throw();
 
+  /// Performs all of the heavy lifting related to initializing cosmology and
+  /// initializing the energy fields from a temperature
+  ///
+  /// This is a public static can be called by other initializers. (There's
+  /// little harm in this function getting called multiple times for a single
+  /// block)
+  static void init_cosmology
+  ( Block * block, double temperature, double gamma ) throw();
+
 private: // functions
 
 

src/Enzo/enzo_EnzoInitialZeldovichPancake.cpp

@@ -0,0 +1,142 @@
+// See LICENSE_CELLO file for license and copyright information
+
+/// @file     enzo_EnzoInitialZeldovichPancake.cpp
+/// @author   Matthew Abruzzo ([email protected])
+/// @date     Sun May 22 2022
+/// @brief    Implementation of Zeldovich Pancake initializer
+
+#include "cello.hpp"
+#include "enzo.hpp"
+
+//----------------------------------------------------------------------
+
+namespace{
+struct coordinate_pair{ double x_eulerian; double x_lagrangian; };
+}
+
+//----------------------------------------------------------------------
+
+void EnzoInitialZeldovichPancake::enforce_block
+( Block * block, const Hierarchy  * hierarchy ) throw()
+{
+  // this function tries to port the initialization routine from
+  //   Grid_ZeldovichPancakeInitializeGrid.C
+  // The initial conditions seem to come from Ryu+ (1993) - they are further
+  // discussed in Collins+ (2010) and Li+ (2008). Unfortunately, it's not
+  // obvious to me what the analytic solution for linear collapse is
+  const double _TOL = 1e-6;
+
+  EnzoPhysicsCosmology * cosmology = enzo::cosmology();
+
+  ASSERT("EnzoInitialZeldovichPancake::enforce_block",
+         "Cosmology must be in use", enzo::cosmology() != nullptr);
+  // todo: add more flexibility or at least be more explicit about the problem
+  ASSERT("EnzoInitialZeldovichPancake::enforce_block",
+         "invalid cosmological parameters",
+         (cosmology->omega_matter_now() == 1.0) &
+         (cosmology->omega_baryon_now() == 1.0) &
+         (cosmology->omega_cdm_now() == 0.0) &
+         (cosmology->omega_lambda_now() == 0.0) );
+
+  const double init_z = cosmology->initial_redshift();
+  const double omega_baryon_now = cosmology->omega_baryon_now();
+
+  // todo: make some/most of the following configurable in the future
+  const double init_temperature_K = 100.0; // in units of kelvin
+  const double pancake_central_offset = 0.0;
+  const double collapse_z = 1.0;
+
+  double x_domain_min, x_domain_max;
+  cello::hierarchy()->lower(&x_domain_min, nullptr, nullptr);
+  cello::hierarchy()->upper(&x_domain_max, nullptr, nullptr);
+  const double lambda = x_domain_max - x_domain_min;
+
+  // 1. precompute some relevant quantities:
+  const double dbl_pi = 2.0 * cello::pi;
+  const double kx = dbl_pi / lambda;
+  const double amplitude = (1.0 + collapse_z) / (1.0 + init_z);
+  const double amplitude_vel =
+    -std::sqrt(2.0/3.0) * (1.0 + collapse_z) / ((1.0 + init_z) * dbl_pi);
+
+  // 2. Construct a function that returns lagrangian and eularian positions
+
+  Data* data = block->data();
+  Field field = data->field();
+  int mx,my,mz;
+  int gx,gy,gz;
+  field.dimensions (field.field_id("density"),&mx,&my,&mz);
+  field.ghost_depth(field.field_id("density"),&gx,&gy,&gz);
+
+  std::vector<double> xc(mx);
+  std::vector<double> dummy(std::max(my, mz));
+  data->field_cells (xc.data(), dummy.data(), dummy.data(), gx, gy, gz);
+  const double pancake_center = (0.5*(x_domain_min + x_domain_max) +
+                                 pancake_central_offset);
+
+  auto get_pos_pair = [&](int ix)
+    {
+      // I don't know why we aren't starting with eulerian coordinates and then
+      // computing lagrangian coordinates... (we're following what Enzo does)
+      double x_lagrange = xc[ix] - pancake_center;
+
+      double x_euler = x_lagrange;
+      double x_euler_old = std::numeric_limits<double>::max();
+      while (fabs((x_euler-x_euler_old)/x_euler) > _TOL) {
+        x_euler_old = x_euler;
+        x_euler += (x_lagrange - x_euler + amplitude*std::sin(kx*x_euler)/kx)
+                  /(1                    - amplitude*std::cos(kx*x_euler)   );
+      }
+      coordinate_pair out = {x_euler, x_lagrange};
+      return out;
+    };
+
+  // todo: use CelloArrays (we aren't using them now to avoid breaking stuff
+  // when we merge in another PR)
+  enzo_float * density = (enzo_float *) field.values("density");
+  enzo_float * e_int = (enzo_float *) field.values("internal_energy");
+  const bool idual = (e_int != nullptr);
+  enzo_float * e_tot = (enzo_float *) field.values("total_energy");
+  enzo_float * vx = (enzo_float *) field.values("velocity_x");
+  enzo_float * vy = (enzo_float *) field.values("velocity_y");
+  enzo_float * vz = (enzo_float *) field.values("velocity_z");
+
+  EnzoUnits * units = enzo::units();
+
+  const double gm1 = gamma_ - 1.0;
+  // convert to problem units
+  const double bulkv = 0.0;
+  const double temperature_units = units->temperature();
+
+  for (int iz=gz; iz<mz-gz; iz++) {
+    for (int iy=gy; iy<my-gy; iy++) {
+      for (int ix=gx; ix<mx-gx; ix++) {
+        coordinate_pair coord_pair = get_pos_pair(ix);
+        double x_euler = coord_pair.x_eulerian;
+
+	int i = ix + mx*(iy + my*iz);
+
+        double cur_density =
+          omega_baryon_now / (1 - amplitude*std::cos(kx*x_euler));
+        double cur_vx = amplitude_vel * sin(kx*x_euler) + bulkv;
+        double cur_eint =
+          ((init_temperature_K / temperature_units) *
+           std::pow(cur_density / omega_baryon_now, gm1) / gm1);
+        double cur_etot = cur_eint + 0.5 * cur_vx * cur_vx;
+
+        density[i] = (enzo_float)cur_density;
+        vx[i] = cur_vx;
+        vy[i] = 0.0;
+        vz[i] = 0.0;
+        e_tot[i] = cur_etot;
+        if (idual) { e_int[i] = cur_eint; }
+
+      }
+    }
+  }
+
+  // set initial time!
+  EnzoInitialCosmology::init_cosmology(block, -1, -1);
+
+  block->initial_done();
+}
+

src/Enzo/enzo_EnzoInitialZeldovichPancake.hpp

@@ -0,0 +1,59 @@
+// See LICENSE_CELLO file for license and copyright information
+
+/// @file     enzo_EnzoInitialZeldovichPancake.hpp
+/// @author   Matthew Abruzzo ([email protected])
+/// @date     Sun May 22 2022
+/// @brief    [\ref Enzo] Initialization routine ZeldovichPancake cosmology
+///           test problem.
+
+#ifndef ENZO_ENZO_INITIAL_ZELDOVICH_PANCAKE_HPP
+#define ENZO_ENZO_INITIAL_ZELDOVICH_PANCAKE_HPP
+
+class EnzoInitialZeldovichPancake : public Initial {
+  /// @class    EnzoInitialZeldovichPancake
+  /// @ingroup  Enzo
+  /// @brief    [\ref Enzo] Initializer for the axis-aligned Zeldovich Pancake
+  /// test problem.
+  ///
+  /// In the future, we might want to combine with the EnzoInitialInclinedWave
+  /// initializer or reuse some of the machinery so that we can incline this
+  /// problem.
+
+public: // interface
+
+  /// Constructor
+  EnzoInitialZeldovichPancake(double gamma, int cycle, double time)
+    : Initial(cycle, time), gamma_(gamma)
+  { }
+
+  /// CHARM++ PUP::able declaration
+  PUPable_decl(EnzoInitialZeldovichPancake);
+
+  /// CHARM++ migration constructor
+  EnzoInitialZeldovichPancake(CkMigrateMessage *m)
+    : Initial (m), gamma_(0.0)
+  {  }
+
+  /// Destructor
+  virtual ~EnzoInitialZeldovichPancake() throw()
+  {  }
+
+  /// CHARM++ Pack / Unpack function
+  void pup (PUP::er &p) {
+    // NOTE: update whenever attributes change
+    TRACEPUP;
+    Initial::pup(p);
+    p | gamma_;
+  }
+
+  /// Initialize the block
+  virtual void enforce_block
+  ( Block * block, const Hierarchy * hierarchy ) throw();
+
+private: // attributes
+
+  /// adiabatic index
+  double gamma_;
+};
+
+#endif /* ENZO_ENZO_INITIAL_ZELDOVICH_PANCAKE_HPP */

src/Enzo/enzo_EnzoProblem.cpp

@@ -254,8 +254,10 @@ Initial * EnzoProblem::create_initial_
     initial = new EnzoInitialIsolatedGalaxy (enzo_config);
   } else if (type == "merge_sinks_test") {
     initial = new EnzoInitialMergeSinksTest (enzo_config);
-  }
-  else {
+  } else if (type == "zeldovich_pancake") {
+    initial = new EnzoInitialZeldovichPancake (enzo_config->field_gamma,
+                                               cycle, time);
+  } else {
     initial = Problem::create_initial_
       (type,index,config,parameters);
   }