diff --git a/input/Cosmology/ZeldovichPancake/vlct-HD-x-aligned.in b/input/Cosmology/ZeldovichPancake/vlct-HD-x-aligned.in new file mode 100644 index 0000000000..1894e020bd --- /dev/null +++ b/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]; + } + } + } \ No newline at end of file diff --git a/src/Enzo/_enzo.hpp b/src/Enzo/_enzo.hpp index 96d9967d15..3c58b973e9 100644 --- a/src/Enzo/_enzo.hpp +++ b/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" diff --git a/src/Enzo/enzo.ci b/src/Enzo/enzo.ci index a9f89155aa..13b679f58c 100644 --- a/src/Enzo/enzo.ci +++ b/src/Enzo/enzo.ci @@ -67,6 +67,7 @@ module enzo { PUPable EnzoInitialTurbulence; PUPable EnzoInitialIsolatedGalaxy; PUPable EnzoInitialMergeSinksTest; + PUPable EnzoInitialZeldovichPancake; PUPable IoEnzoBlock; diff --git a/src/Enzo/enzo_EnzoInitialCosmology.cpp b/src/Enzo/enzo_EnzoInitialCosmology.cpp index b82ebc4218..68155683d7 100644 --- a/src/Enzo/enzo_EnzoInitialCosmology.cpp +++ b/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(); - } diff --git a/src/Enzo/enzo_EnzoInitialCosmology.hpp b/src/Enzo/enzo_EnzoInitialCosmology.hpp index e6698bad88..6876d32f21 100644 --- a/src/Enzo/enzo_EnzoInitialCosmology.hpp +++ b/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 diff --git a/src/Enzo/enzo_EnzoInitialZeldovichPancake.cpp b/src/Enzo/enzo_EnzoInitialZeldovichPancake.cpp new file mode 100644 index 0000000000..81600261e0 --- /dev/null +++ b/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 (matthewabruzzo@gmail.com) +/// @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 xc(mx); + std::vector 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::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; izinitial_done(); +} + diff --git a/src/Enzo/enzo_EnzoInitialZeldovichPancake.hpp b/src/Enzo/enzo_EnzoInitialZeldovichPancake.hpp new file mode 100644 index 0000000000..6111e9eeca --- /dev/null +++ b/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 (matthewabruzzo@gmail.com) +/// @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 */ diff --git a/src/Enzo/enzo_EnzoProblem.cpp b/src/Enzo/enzo_EnzoProblem.cpp index 035596f6db..634e9cabde 100644 --- a/src/Enzo/enzo_EnzoProblem.cpp +++ b/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); }