Skip to content

Commit ad46661

Browse files
mabruzzomabruzzo
committed
Introduce an initializer for Zeldovich Pancake test problem.
1 parent dfb8f5e commit ad46661

8 files changed

+328
-9
lines changed

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

+92
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
include "input/vlct/vl.incl"
2+
3+
Adapt {
4+
max_level = 0;
5+
}
6+
7+
Method {
8+
9+
list = ["cosmology", "pm_deposit", "gravity", "mhd_vlct", "comoving_expansion"];
10+
11+
gravity {
12+
solver = "cg";
13+
grav_const = 1.0;
14+
}
15+
}
16+
17+
Initial{ list = ["zeldovich_pancake"]; }
18+
19+
Boundary { type = "periodic"; }
20+
21+
Domain {
22+
lower = [ 0.0, 0.0, 0.0 ];
23+
upper = [ 1.0, 0.0625, 0.0625 ];
24+
}
25+
26+
Mesh {
27+
root_rank = 3;
28+
root_size = [ 64, 4, 4 ];
29+
root_blocks = [1, 1, 1];
30+
}
31+
32+
Physics {
33+
list = [ "cosmology" ];
34+
cosmology {
35+
omega_baryon_now = 1.0;
36+
omega_cdm_now = 0.0;
37+
omega_matter_now = 1.0;
38+
omega_lambda_now = 0.0;
39+
hubble_constant_now = 0.5;
40+
max_expansion_rate = 0.01;
41+
initial_redshift = 20.0;
42+
final_redshift = 0.0;
43+
comoving_box_size = 64.0; # specifies the code length unit in terms
44+
# of comoving Mpc / h
45+
}
46+
}
47+
48+
Field {
49+
gamma = 1.6666666666666667;
50+
51+
alignment = 8;
52+
ghost_depth = 4;
53+
history = 1;
54+
list += ["internal_energy", "acceleration_x", "acceleration_y", "acceleration_z",
55+
"density_gas", "density_particle", "density_total",
56+
"density_particle_accumulate", "potential", "potential_temp",
57+
"potential_copy", "X", "B", "X_copy", "B_copy", "particle_mass"];
58+
padding = 0;
59+
}
60+
61+
62+
63+
Solver {
64+
list = ["cg"];
65+
cg {
66+
type = "cg";
67+
iter_max = 1000;
68+
res_tol = 1e-14;
69+
monitor_iter = 25;
70+
solve_type = "block";
71+
}
72+
}
73+
74+
Stopping{
75+
#close to redshift= 1
76+
time = 2.782174164581e+01;
77+
}
78+
79+
Output {
80+
list = ["data"];
81+
data {
82+
type = "data";
83+
field_list = ["density", "velocity_x", "velocity_y", "velocity_z",
84+
"total_energy"];
85+
name = ["data-%03d.h5", "proc"];
86+
dir = ["method-vlct-x-aligned-pancake_%04d","cycle"];
87+
schedule {
88+
var = "cycle";
89+
list = [0, 1, 2, 3, 4, 5, 6, 7, 8]; # 237];
90+
}
91+
}
92+
}

‎src/Enzo/_enzo.hpp

+1
Original file line numberDiff line numberDiff line change
@@ -201,6 +201,7 @@ extern "C" {
201201
#include "enzo_EnzoInitialIsolatedGalaxy.hpp"
202202
#include "enzo_EnzoInitialBurkertBodenheimer.hpp"
203203
#include "enzo_EnzoInitialMergeSinksTest.hpp"
204+
#include "enzo_EnzoInitialZeldovichPancake.hpp"
204205

205206
#include "enzo_EnzoRefineShock.hpp"
206207
#include "enzo_EnzoRefineParticleMass.hpp"

‎src/Enzo/enzo.ci

+1
Original file line numberDiff line numberDiff line change
@@ -67,6 +67,7 @@ module enzo {
6767
PUPable EnzoInitialTurbulence;
6868
PUPable EnzoInitialIsolatedGalaxy;
6969
PUPable EnzoInitialMergeSinksTest;
70+
PUPable EnzoInitialZeldovichPancake;
7071

7172
PUPable IoEnzoBlock;
7273

‎src/Enzo/enzo_EnzoInitialCosmology.cpp

+20-7
Original file line numberDiff line numberDiff line change
@@ -26,7 +26,6 @@ EnzoInitialCosmology::EnzoInitialCosmology
2626
void EnzoInitialCosmology::enforce_block
2727
( Block * block, const Hierarchy * hierarchy ) throw()
2828
{
29-
EnzoUnits * units = enzo::units();
3029

3130
// "If temperature is left unset, set it assuming that T=550 K at z=200"
3231
EnzoPhysicsCosmology * cosmology = enzo::cosmology();
@@ -36,16 +35,33 @@ void EnzoInitialCosmology::enforce_block
3635
pow((1.0 + cosmology->initial_redshift())/(1.0 + 200.00), 2.0);
3736
}
3837

38+
EnzoInitialCosmology::init_cosmology(block, temperature_, gamma_);
39+
40+
block->initial_done();
41+
42+
}
43+
44+
//----------------------------------------------------------------------
45+
46+
void EnzoInitialCosmology::init_cosmology
47+
(Block * block, double temperature, double gamma) throw()
48+
{
49+
EnzoPhysicsCosmology * cosmology = enzo::cosmology();
50+
3951
// Set initial time based on initial redshift
4052
double r0 = cosmology->initial_redshift();
4153
double t0 = cosmology->time_from_redshift(r0);
4254
block->set_time(t0);
4355
enzo::simulation()->set_time(t0);
44-
56+
57+
if (temperature < 0){ return; }
58+
59+
EnzoUnits * units = enzo::units();
60+
// initialize energy fields
4561
const double default_mu = 0.6;
4662

47-
const double internal_energy = temperature_/units->temperature()
48-
/default_mu/(gamma_-1.0);
63+
const double internal_energy = temperature/units->temperature()
64+
/default_mu/(gamma-1.0);
4965

5066
Field field = block->data()->field();
5167

@@ -70,7 +86,4 @@ void EnzoInitialCosmology::enforce_block
7086
}
7187
}
7288
}
73-
74-
block->initial_done();
75-
7689
}

‎src/Enzo/enzo_EnzoInitialCosmology.hpp

+9
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,15 @@ class EnzoInitialCosmology : public Initial {
5050
virtual void enforce_block
5151
( Block * block, const Hierarchy * hierarchy ) throw();
5252

53+
/// Performs all of the heavy lifting related to initializing cosmology and
54+
/// initializing the energy fields from a temperature
55+
///
56+
/// This is a public static can be called by other initializers. (There's
57+
/// little harm in this function getting called multiple times for a single
58+
/// block)
59+
static void init_cosmology
60+
( Block * block, double temperature, double gamma ) throw();
61+
5362
private: // functions
5463

5564

‎src/Enzo/enzo_EnzoInitialZeldovichPancake.cpp

+142
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,142 @@
1+
// See LICENSE_CELLO file for license and copyright information
2+
3+
/// @file enzo_EnzoInitialZeldovichPancake.cpp
4+
/// @author Matthew Abruzzo ([email protected])
5+
/// @date Sun May 22 2022
6+
/// @brief Implementation of Zeldovich Pancake initializer
7+
8+
#include "cello.hpp"
9+
#include "enzo.hpp"
10+
11+
//----------------------------------------------------------------------
12+
13+
namespace{
14+
struct coordinate_pair{ double x_eulerian; double x_lagrangian; };
15+
}
16+
17+
//----------------------------------------------------------------------
18+
19+
void EnzoInitialZeldovichPancake::enforce_block
20+
( Block * block, const Hierarchy * hierarchy ) throw()
21+
{
22+
// this function tries to port the initialization routine from
23+
// Grid_ZeldovichPancakeInitializeGrid.C
24+
// The initial conditions seem to come from Ryu+ (1993) - they are further
25+
// discussed in Collins+ (2010) and Li+ (2008). Unfortunately, it's not
26+
// obvious to me what the analytic solution for linear collapse is
27+
const double _TOL = 1e-6;
28+
29+
EnzoPhysicsCosmology * cosmology = enzo::cosmology();
30+
31+
ASSERT("EnzoInitialZeldovichPancake::enforce_block",
32+
"Cosmology must be in use", enzo::cosmology() != nullptr);
33+
// todo: add more flexibility or at least be more explicit about the problem
34+
ASSERT("EnzoInitialZeldovichPancake::enforce_block",
35+
"invalid cosmological parameters",
36+
(cosmology->omega_matter_now() == 1.0) &
37+
(cosmology->omega_baryon_now() == 1.0) &
38+
(cosmology->omega_cdm_now() == 0.0) &
39+
(cosmology->omega_lambda_now() == 0.0) );
40+
41+
const double init_z = cosmology->initial_redshift();
42+
const double omega_baryon_now = cosmology->omega_baryon_now();
43+
44+
// todo: make some/most of the following configurable in the future
45+
const double init_temperature_K = 100.0; // in units of kelvin
46+
const double pancake_central_offset = 0.0;
47+
const double collapse_z = 1.0;
48+
49+
double x_domain_min, x_domain_max;
50+
cello::hierarchy()->lower(&x_domain_min, nullptr, nullptr);
51+
cello::hierarchy()->upper(&x_domain_max, nullptr, nullptr);
52+
const double lambda = x_domain_max - x_domain_min;
53+
54+
// 1. precompute some relevant quantities:
55+
const double dbl_pi = 2.0 * cello::pi;
56+
const double kx = dbl_pi / lambda;
57+
const double amplitude = (1.0 + collapse_z) / (1.0 + init_z);
58+
const double amplitude_vel =
59+
-std::sqrt(2.0/3.0) * (1.0 + collapse_z) / ((1.0 + init_z) * dbl_pi);
60+
61+
// 2. Construct a function that returns lagrangian and eularian positions
62+
63+
Data* data = block->data();
64+
Field field = data->field();
65+
int mx,my,mz;
66+
int gx,gy,gz;
67+
field.dimensions (field.field_id("density"),&mx,&my,&mz);
68+
field.ghost_depth(field.field_id("density"),&gx,&gy,&gz);
69+
70+
std::vector<double> xc(mx);
71+
std::vector<double> dummy(std::max(my, mz));
72+
data->field_cells (xc.data(), dummy.data(), dummy.data(), gx, gy, gz);
73+
const double pancake_center = (0.5*(x_domain_min + x_domain_max) +
74+
pancake_central_offset);
75+
76+
auto get_pos_pair = [&](int ix)
77+
{
78+
// I don't know why we aren't starting with eulerian coordinates and then
79+
// computing lagrangian coordinates... (we're following what Enzo does)
80+
double x_lagrange = xc[ix] - pancake_center;
81+
82+
double x_euler = x_lagrange;
83+
double x_euler_old = std::numeric_limits<double>::max();
84+
while (fabs((x_euler-x_euler_old)/x_euler) > _TOL) {
85+
x_euler_old = x_euler;
86+
x_euler += (x_lagrange - x_euler + amplitude*std::sin(kx*x_euler)/kx)
87+
/(1 - amplitude*std::cos(kx*x_euler) );
88+
}
89+
coordinate_pair out = {x_euler, x_lagrange};
90+
return out;
91+
};
92+
93+
// todo: use CelloArrays (we aren't using them now to avoid breaking stuff
94+
// when we merge in another PR)
95+
enzo_float * density = (enzo_float *) field.values("density");
96+
enzo_float * e_int = (enzo_float *) field.values("internal_energy");
97+
const bool idual = (e_int != nullptr);
98+
enzo_float * e_tot = (enzo_float *) field.values("total_energy");
99+
enzo_float * vx = (enzo_float *) field.values("velocity_x");
100+
enzo_float * vy = (enzo_float *) field.values("velocity_y");
101+
enzo_float * vz = (enzo_float *) field.values("velocity_z");
102+
103+
EnzoUnits * units = enzo::units();
104+
105+
const double gm1 = gamma_ - 1.0;
106+
// convert to problem units
107+
const double bulkv = 0.0;
108+
const double temperature_units = units->temperature();
109+
110+
for (int iz=gz; iz<mz-gz; iz++) {
111+
for (int iy=gy; iy<my-gy; iy++) {
112+
for (int ix=gx; ix<mx-gx; ix++) {
113+
coordinate_pair coord_pair = get_pos_pair(ix);
114+
double x_euler = coord_pair.x_eulerian;
115+
116+
int i = ix + mx*(iy + my*iz);
117+
118+
double cur_density =
119+
omega_baryon_now / (1 - amplitude*std::cos(kx*x_euler));
120+
double cur_vx = amplitude_vel * sin(kx*x_euler) + bulkv;
121+
double cur_eint =
122+
((init_temperature_K / temperature_units) *
123+
std::pow(cur_density / omega_baryon_now, gm1) / gm1);
124+
double cur_etot = cur_eint + 0.5 * cur_vx * cur_vx;
125+
126+
density[i] = (enzo_float)cur_density;
127+
vx[i] = cur_vx;
128+
vy[i] = 0.0;
129+
vz[i] = 0.0;
130+
e_tot[i] = cur_etot;
131+
if (idual) { e_int[i] = cur_eint; }
132+
133+
}
134+
}
135+
}
136+
137+
// set initial time!
138+
EnzoInitialCosmology::init_cosmology(block, -1, -1);
139+
140+
block->initial_done();
141+
}
142+

‎src/Enzo/enzo_EnzoInitialZeldovichPancake.hpp

+59
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
// See LICENSE_CELLO file for license and copyright information
2+
3+
/// @file enzo_EnzoInitialZeldovichPancake.hpp
4+
/// @author Matthew Abruzzo ([email protected])
5+
/// @date Sun May 22 2022
6+
/// @brief [\ref Enzo] Initialization routine ZeldovichPancake cosmology
7+
/// test problem.
8+
9+
#ifndef ENZO_ENZO_INITIAL_ZELDOVICH_PANCAKE_HPP
10+
#define ENZO_ENZO_INITIAL_ZELDOVICH_PANCAKE_HPP
11+
12+
class EnzoInitialZeldovichPancake : public Initial {
13+
/// @class EnzoInitialZeldovichPancake
14+
/// @ingroup Enzo
15+
/// @brief [\ref Enzo] Initializer for the axis-aligned Zeldovich Pancake
16+
/// test problem.
17+
///
18+
/// In the future, we might want to combine with the EnzoInitialInclinedWave
19+
/// initializer or reuse some of the machinery so that we can incline this
20+
/// problem.
21+
22+
public: // interface
23+
24+
/// Constructor
25+
EnzoInitialZeldovichPancake(double gamma, int cycle, double time)
26+
: Initial(cycle, time), gamma_(gamma)
27+
{ }
28+
29+
/// CHARM++ PUP::able declaration
30+
PUPable_decl(EnzoInitialZeldovichPancake);
31+
32+
/// CHARM++ migration constructor
33+
EnzoInitialZeldovichPancake(CkMigrateMessage *m)
34+
: Initial (m), gamma_(0.0)
35+
{ }
36+
37+
/// Destructor
38+
virtual ~EnzoInitialZeldovichPancake() throw()
39+
{ }
40+
41+
/// CHARM++ Pack / Unpack function
42+
void pup (PUP::er &p) {
43+
// NOTE: update whenever attributes change
44+
TRACEPUP;
45+
Initial::pup(p);
46+
p | gamma_;
47+
}
48+
49+
/// Initialize the block
50+
virtual void enforce_block
51+
( Block * block, const Hierarchy * hierarchy ) throw();
52+
53+
private: // attributes
54+
55+
/// adiabatic index
56+
double gamma_;
57+
};
58+
59+
#endif /* ENZO_ENZO_INITIAL_ZELDOVICH_PANCAKE_HPP */

‎src/Enzo/enzo_EnzoProblem.cpp

+4-2
Original file line numberDiff line numberDiff line change
@@ -254,8 +254,10 @@ Initial * EnzoProblem::create_initial_
254254
initial = new EnzoInitialIsolatedGalaxy (enzo_config);
255255
} else if (type == "merge_sinks_test") {
256256
initial = new EnzoInitialMergeSinksTest (enzo_config);
257-
}
258-
else {
257+
} else if (type == "zeldovich_pancake") {
258+
initial = new EnzoInitialZeldovichPancake (enzo_config->field_gamma,
259+
cycle, time);
260+
} else {
259261
initial = Problem::create_initial_
260262
(type,index,config,parameters);
261263
}

0 commit comments

Comments
 (0)