mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
commit
bfce9dc8c7
7 changed files with 230 additions and 125 deletions
120
example.conf
120
example.conf
|
@ -1,95 +1,112 @@
|
||||||
#########################################################################################
|
#########################################################################################
|
||||||
# Example conf file for MUSIC2 - monofonIC single resolution simulation ICs
|
# Example config file for MUSIC2 - monofonIC single resolution simulation ICs
|
||||||
|
# version 1 from 2020/08/23
|
||||||
#########################################################################################
|
#########################################################################################
|
||||||
|
|
||||||
#########################################################################################
|
#########################################################################################
|
||||||
[setup]
|
[setup]
|
||||||
# number of grid cells per linear dimension for calculations = particles for sc initial load
|
|
||||||
GridRes = 128
|
|
||||||
|
|
||||||
# length of the box in Mpc/h
|
GridRes = 128 # number of grid cells per linear dimension for calculations
|
||||||
BoxLength = 300
|
# = particles for sc initial load
|
||||||
|
BoxLength = 300 # length of the box in Mpc/h
|
||||||
|
zstart = 24.0 # starting redshift
|
||||||
|
|
||||||
# starting redshift
|
LPTorder = 3 # order of the LPT to be used (1,2 or 3)
|
||||||
zstart = 24.0
|
|
||||||
|
|
||||||
# order of the LPT to be used (1,2 or 3)
|
DoBaryons = no # also do baryon ICs?
|
||||||
LPTorder = 3
|
|
||||||
|
|
||||||
# also do baryon ICs?
|
DoFixing = yes # do mode fixing à la Angulo&Pontzen (https://arxiv.org/abs/1603.05253)
|
||||||
DoBaryons = no
|
DoInversion = no # invert phases (for paired simulations)
|
||||||
|
|
||||||
# do mode fixing à la Angulo&Pontzen
|
ParticleLoad = sc # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x)
|
||||||
DoFixing = yes
|
# (increases number of particles by given factor!), or 'glass'
|
||||||
|
|
||||||
# invert phases (for paired simulations)
|
|
||||||
DoInversion = no
|
|
||||||
|
|
||||||
# particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x) (increases number of particles by factor!) or 'glass'
|
|
||||||
ParticleLoad = sc
|
|
||||||
|
|
||||||
# if `ParticleLoad = glass' then specify here where to load the glass distribution from
|
|
||||||
#GlassFileName = glass128.hdf5
|
|
||||||
#GlassTiles = 1
|
|
||||||
|
|
||||||
|
## if `ParticleLoad = glass' then specify here where to load the glass distribution from
|
||||||
|
# GlassFileName = glass128.hdf5
|
||||||
|
# GlassTiles = 1
|
||||||
|
|
||||||
#########################################################################################
|
#########################################################################################
|
||||||
[cosmology]
|
[cosmology]
|
||||||
## transfer = ... specifies the Einstein-Boltzmann plugin module
|
## transfer = ... specifies the Einstein-Boltzmann plugin module
|
||||||
# transfer = eisenstein # Eisenstein&Hu fitting formula
|
|
||||||
|
|
||||||
# transfer = file_CAMB # CAMB file to be specified as 'transfer_file = ...'
|
|
||||||
# transfer_file = wmap5_transfer_out_z0.dat
|
|
||||||
|
|
||||||
transfer = CLASS # CLASS module (if enabled in CMake file)
|
|
||||||
|
|
||||||
ztarget = 2.5 # target redshift for CLASS module, output at ztarget will be back-scaled to zstart
|
|
||||||
|
|
||||||
|
# main cosmological parameters
|
||||||
Omega_m = 0.302
|
Omega_m = 0.302
|
||||||
Omega_b = 0.045
|
Omega_b = 0.045
|
||||||
Omega_L = 0.698
|
Omega_L = 0.698
|
||||||
H0 = 70.3
|
H0 = 70.3
|
||||||
nspec = 0.961
|
nspec = 0.961
|
||||||
sigma_8 = 0.811
|
sigma_8 = 0.811
|
||||||
# A_s = 2.148752e-09 # can use instead of sigma_8
|
|
||||||
|
|
||||||
ZeroRadiation = true # For Back-scaling: set to false if your simulation code can deal with Omega_r!=0
|
ZeroRadiation = false # For Back-scaling only: set to true if your simulation code
|
||||||
|
# can deal with Omega_r!=0 in its background FLRW model
|
||||||
|
|
||||||
# Additional cosmological parameters (set by default to the given values)
|
## Additional cosmological parameters (set by default to the given values)
|
||||||
# w0 = -1.0
|
# w0 = -1.0
|
||||||
# wa = 0.0
|
# wa = 0.0
|
||||||
# Tcmb = 2.7255
|
# Tcmb = 2.7255
|
||||||
# Neff = 3.046
|
# Neff = 3.046
|
||||||
|
|
||||||
# anisotropic large scale tidal field
|
## Use below for anisotropic large scale tidal field ICs up to 2LPT
|
||||||
# see Stuecker+2020
|
## see Stuecker+2020 (https://arxiv.org/abs/2003.06427)
|
||||||
# LSS_aniso_lx = +0.1
|
# LSS_aniso_lx = +0.1
|
||||||
# LSS_aniso_ly = +0.1
|
# LSS_aniso_ly = +0.1
|
||||||
# LSS_aniso_lz = -0.2
|
# LSS_aniso_lz = -0.2
|
||||||
|
|
||||||
|
### MODULES: ###
|
||||||
|
|
||||||
|
##> Eisenstein & Hu (1997) fitting formulae
|
||||||
|
## this is fast, but not too accurate. Also baryons trace CDM here.
|
||||||
|
## see https://arxiv.org/abs/astro-ph/9709112
|
||||||
|
|
||||||
|
# transfer = eisenstein
|
||||||
|
|
||||||
|
##> CAMB transfer function file module
|
||||||
|
## This should be transfer function output with CAMB (https://camb.info)
|
||||||
|
## at the *target* redshift
|
||||||
|
|
||||||
|
# transfer = file_CAMB # CAMB file to be specified as 'transfer_file = ...'
|
||||||
|
# transfer_file = wmap5_transfer_out_z0.dat
|
||||||
|
|
||||||
|
##> CLASS module, which links to the actual CLASS C-code.
|
||||||
|
## note that CLASS needs to be cloned as a git submodule and enabled in CMake file
|
||||||
|
|
||||||
|
transfer = CLASS
|
||||||
|
ztarget = 2.5 # target redshift for CLASS module, output at ztarget will be back-scaled to zstart
|
||||||
|
# A_s = 2.148752e-09 # can use A_s instead of sigma_8 when using CLASS
|
||||||
|
|
||||||
|
|
||||||
#########################################################################################
|
#########################################################################################
|
||||||
[random]
|
[random]
|
||||||
## generator = ... specifies the random field generator plugin module
|
## generator = ... specifies the random field generator plugin module
|
||||||
generator = NGENIC
|
|
||||||
seed = 9001
|
|
||||||
|
|
||||||
# generator = PANPHASIA
|
##> NGenIC compatible random number generator module compatible with V. Springel's original code
|
||||||
# descriptor = [Panph1,L10,(800,224,576),S9,CH1564365824,MXXL]
|
## (https://www.h-its.org/2014/11/05/ngenic-code/) as well as the 2LPT code by Pueblas&Scoccmiarro
|
||||||
|
## (https://cosmo.nyu.edu/roman/2LPT/)
|
||||||
|
generator = NGENIC
|
||||||
|
seed = 12345
|
||||||
|
|
||||||
|
##> The PANPHASIA generator uses a plugin based on original code by A. Jenkins
|
||||||
|
## Warning: Before using this module, please make sure you read and agree to the distinct license
|
||||||
|
## requirements by registering on the website http://icc.dur.ac.uk/Panphasia.php
|
||||||
|
|
||||||
|
# generator = PANPHASIA
|
||||||
|
# descriptor = [Panph1,L10,(800,224,576),S9,CH1564365824,MXXL]
|
||||||
|
|
||||||
|
##> The MUSIC1 multi-scale random number generator is provided for convenience
|
||||||
|
## warning: MUSIC1 generator is not MPI parallel (yet) (memory is needed for full field on each task)
|
||||||
|
# generator = MUSIC1
|
||||||
|
# seed[7] = 12345
|
||||||
|
# seed[8] = 23456
|
||||||
|
# seed[9] = 34567
|
||||||
|
|
||||||
# Add a possible constraint field here:
|
# Add a possible constraint field here:
|
||||||
# ConstraintFieldFile = initial_conditions.hdf5
|
# ConstraintFieldFile = initial_conditions.hdf5
|
||||||
# ConstraintFieldName = ic_white_noise
|
# ConstraintFieldName = ic_white_noise
|
||||||
|
|
||||||
|
|
||||||
#########################################################################################
|
|
||||||
[testing]
|
|
||||||
# enables diagnostic output
|
|
||||||
# can be 'none' (default), 'potentials_and_densities', 'velocity_displacement_symmetries', or 'convergence'
|
|
||||||
test = none
|
|
||||||
|
|
||||||
|
|
||||||
#########################################################################################
|
#########################################################################################
|
||||||
[execution]
|
[execution]
|
||||||
|
# Specify the number of threads / task
|
||||||
NumThreads = 8
|
NumThreads = 8
|
||||||
|
|
||||||
|
|
||||||
|
@ -97,20 +114,29 @@ NumThreads = 8
|
||||||
[output]
|
[output]
|
||||||
## format = .... specifies the output plugin module
|
## format = .... specifies the output plugin module
|
||||||
|
|
||||||
|
##> RAMSES / GRAFIC2 compatible format
|
||||||
# format = grafic2
|
# format = grafic2
|
||||||
# filename = ics_ramses
|
# filename = ics_ramses
|
||||||
# grafic_use_SPT = no # if no then uses PPT, otherwise linear SPT
|
# grafic_use_SPT = no # if no then uses PPT, otherwise linear SPT
|
||||||
|
|
||||||
|
##> Gadget-2/3 'fortran unformatted binary'-style format
|
||||||
# format = gadget2
|
# format = gadget2
|
||||||
# filename = ics_gadget.dat
|
# filename = ics_gadget.dat
|
||||||
# UseLongids = false
|
# UseLongids = false
|
||||||
|
|
||||||
|
##> Gadget-2/3 HDF5 format
|
||||||
# format = gadget_hdf5
|
# format = gadget_hdf5
|
||||||
# filename = ics_gadget.hdf5
|
# filename = ics_gadget.hdf5
|
||||||
|
|
||||||
|
##> Arepo HDF5 format (virtually identical to gadget_hdf5)
|
||||||
# format = AREPO
|
# format = AREPO
|
||||||
# filename = ics_arepo.hdf5
|
# filename = ics_arepo.hdf5
|
||||||
|
|
||||||
|
##> HACC compatible generic-io format
|
||||||
|
# format = genericio
|
||||||
|
# filename = ics_hacc
|
||||||
|
|
||||||
|
##> Generic HDF5 output format for testing or PT-based calculations
|
||||||
# format = generic
|
# format = generic
|
||||||
# filename = debug.hdf5
|
# filename = debug.hdf5
|
||||||
# generic_out_eulerian = yes # if yes then uses PPT for output
|
# generic_out_eulerian = yes # if yes then uses PPT for output
|
||||||
|
|
35
external/panphasia/LICENSE
vendored
Normal file
35
external/panphasia/LICENSE
vendored
Normal file
|
@ -0,0 +1,35 @@
|
||||||
|
The code in this subdirectory is part of Adrian Jenkins' PANPHASIA packet,
|
||||||
|
obtained from here http://icc.dur.ac.uk/Panphasia.php
|
||||||
|
|
||||||
|
PANPHASIA is not published under the GPL but has its own proprietary license,
|
||||||
|
make sure to visit the website before using the PANPHASIA functionality of
|
||||||
|
MUSIC2 and register your name.
|
||||||
|
|
||||||
|
We reproduce the licensing requirements for PANPHASIA from the above website
|
||||||
|
as retrieved on 2020/08/23:
|
||||||
|
|
||||||
|
We make our software available for free but with a licence that includes the
|
||||||
|
condition that users make sure the phases of any new simulation volumes set up
|
||||||
|
using Panphasia are published.
|
||||||
|
|
||||||
|
We are happy to collaborate with others on improving the software and providing
|
||||||
|
support for languages other than fortran. Contact: A.R.Jenkins@durham.ac.uk
|
||||||
|
|
||||||
|
LICENCE:
|
||||||
|
|
||||||
|
You are licensed to use this software free of charge on condition that:
|
||||||
|
|
||||||
|
- you will publish the phase descriptors and reference Jenkins (13) for any new
|
||||||
|
simulations that use Panphasia phases. You will pass on this condition to others
|
||||||
|
for any software or data you make available publically or privately that makes
|
||||||
|
use of Panphasia.
|
||||||
|
- that you will ensure any publications using results derived from Panphasia will
|
||||||
|
be submitted as a final version to arXiv prior to or coincident with publication
|
||||||
|
in a journal.
|
||||||
|
- that you report any bugs in this software as soon as confirmed to
|
||||||
|
A.R.Jenkins@durham.ac.uk
|
||||||
|
- that you understand that the software comes with no warranty and that is your
|
||||||
|
responsibility to ensure that it is suitable for the purpose that you intend.
|
||||||
|
- that you agree to having your name and email address stored for an indefinite
|
||||||
|
period in the future electronically in a database as a record that you agreed
|
||||||
|
the licence conditions.
|
|
@ -24,6 +24,7 @@ class RNG_music : public RNG_plugin
|
||||||
protected:
|
protected:
|
||||||
std::vector<long> rngseeds_;
|
std::vector<long> rngseeds_;
|
||||||
std::vector<std::string> rngfnames_;
|
std::vector<std::string> rngfnames_;
|
||||||
|
std::vector<rng *> randc_;
|
||||||
unsigned ran_cube_size_;
|
unsigned ran_cube_size_;
|
||||||
|
|
||||||
int levelmin_, levelmax_, levelmin_seed_;
|
int levelmin_, levelmax_, levelmin_seed_;
|
||||||
|
@ -49,21 +50,73 @@ protected:
|
||||||
//! store the white noise fields in memory or on disk
|
//! store the white noise fields in memory or on disk
|
||||||
//void store_rnd(int ilevel, rng *prng);
|
//void store_rnd(int ilevel, rng *prng);
|
||||||
|
|
||||||
public:
|
bool is_power_of_two(size_t x) const
|
||||||
explicit RNG_music(config_file &cf) : RNG_plugin(cf), initialized_(false) {}
|
{
|
||||||
|
return (x & (x - 1)) == 0;
|
||||||
|
}
|
||||||
|
|
||||||
~RNG_music() {}
|
unsigned int_log2( size_t v ) const
|
||||||
|
{
|
||||||
|
unsigned r{0}; // r will be lg(v)
|
||||||
|
while (v >>= 1) r++;
|
||||||
|
return r;
|
||||||
|
}
|
||||||
|
|
||||||
|
public:
|
||||||
|
explicit RNG_music(config_file &cf) : RNG_plugin(cf), initialized_(false)
|
||||||
|
{
|
||||||
|
// we need to make sure that the chosen resolution is a power of 2 resolution
|
||||||
|
size_t res = pcf_->get_value<size_t>("setup", "GridRes");
|
||||||
|
if( !is_power_of_two(res) ){
|
||||||
|
std::string errmsg("MUSIC random number plugin requires [setup]/GridRes to be a power of 2!");
|
||||||
|
music::flog << errmsg << std::endl;
|
||||||
|
throw std::runtime_error(errmsg.c_str());
|
||||||
|
}
|
||||||
|
levelmin_ = int_log2(res);
|
||||||
|
levelmax_ = levelmin_;
|
||||||
|
music::ilog << "MUSIC1 RNG plugin: setting levelmin = levelmax = " << levelmin_ << std::endl;
|
||||||
|
|
||||||
|
this->initialize_for_grid_structure( );
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
~RNG_music()
|
||||||
|
{
|
||||||
|
for( auto rr : randc_ ){
|
||||||
|
if( rr != nullptr ) delete rr;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
bool isMultiscale() const { return true; }
|
bool isMultiscale() const { return true; }
|
||||||
|
|
||||||
void Fill_Grid( Grid_FFT<real_t>& g ) {} //const { }
|
void Fill_Grid( Grid_FFT<real_t>& g )
|
||||||
|
|
||||||
void initialize_for_grid_structure()//const refinement_hierarchy &refh)
|
|
||||||
{
|
{
|
||||||
//prefh_ = &refh;
|
// determine extent of grid to be filled (can be a slab with MPI)
|
||||||
levelmin_ = pcf_->get_value<unsigned>("setup", "levelmin");
|
const size_t i0 = g.local_0_start_, j0{0}, k0{0};
|
||||||
levelmax_ = pcf_->get_value<unsigned>("setup", "levelmax");
|
const size_t Ni = g.rsize(0), Nj = g.rsize(1), Nk = g.rsize(2);
|
||||||
|
|
||||||
|
// make sure we're in real space
|
||||||
|
g.FourierTransformBackward();
|
||||||
|
|
||||||
|
// copy over
|
||||||
|
#pragma omp parallel for
|
||||||
|
for( auto i = i0; i<Ni; ++i )
|
||||||
|
{
|
||||||
|
size_t ip = i-i0; // index in g
|
||||||
|
for( auto j = j0; j<Nj; ++j )
|
||||||
|
{
|
||||||
|
auto jp = j-j0; // index in g
|
||||||
|
for( auto k = k0; k<Nk; ++k )
|
||||||
|
{
|
||||||
|
auto kp = k-k0; // index in g
|
||||||
|
g.relem(ip,jp,kp) = (*randc_[levelmin_])(i,j,k);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
void initialize_for_grid_structure()
|
||||||
|
{
|
||||||
ran_cube_size_ = pcf_->get_value_safe<unsigned>("random", "cubesize", DEF_RAN_CUBE_SIZE);
|
ran_cube_size_ = pcf_->get_value_safe<unsigned>("random", "cubesize", DEF_RAN_CUBE_SIZE);
|
||||||
disk_cached_ = pcf_->get_value_safe<bool>("random", "disk_cached", true);
|
disk_cached_ = pcf_->get_value_safe<bool>("random", "disk_cached", true);
|
||||||
restart_ = pcf_->get_value_safe<bool>("random", "restart", false);
|
restart_ = pcf_->get_value_safe<bool>("random", "restart", false);
|
||||||
|
@ -159,17 +212,18 @@ void RNG_music::compute_random_numbers(void)
|
||||||
{
|
{
|
||||||
bool rndsign = pcf_->get_value_safe<bool>("random", "grafic_sign", false);
|
bool rndsign = pcf_->get_value_safe<bool>("random", "grafic_sign", false);
|
||||||
|
|
||||||
std::vector<rng *> randc(std::max(levelmax_, levelmin_seed_) + 1, (rng *)NULL);
|
|
||||||
|
|
||||||
//--- FILL ALL WHITE NOISE ARRAYS FOR WHICH WE NEED THE FULL FIELD ---//
|
//--- FILL ALL WHITE NOISE ARRAYS FOR WHICH WE NEED THE FULL FIELD ---//
|
||||||
|
|
||||||
|
randc_.assign(std::max(levelmax_, levelmin_seed_) + 1, nullptr);
|
||||||
|
|
||||||
|
|
||||||
//... seeds are given for a level coarser than levelmin
|
//... seeds are given for a level coarser than levelmin
|
||||||
if (levelmin_seed_ < levelmin_)
|
if (levelmin_seed_ < levelmin_)
|
||||||
{
|
{
|
||||||
if (rngfnames_[levelmin_seed_].size() > 0)
|
if (rngfnames_[levelmin_seed_].size() > 0)
|
||||||
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
|
randc_[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
|
||||||
else
|
else
|
||||||
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_], true);
|
randc_[levelmin_seed_] = new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_], true);
|
||||||
|
|
||||||
for (int i = levelmin_seed_ + 1; i <= levelmin_; ++i)
|
for (int i = levelmin_seed_ + 1; i <= levelmin_; ++i)
|
||||||
{
|
{
|
||||||
|
@ -178,9 +232,9 @@ void RNG_music::compute_random_numbers(void)
|
||||||
if (rngfnames_[i].size() > 0)
|
if (rngfnames_[i].size() > 0)
|
||||||
music::ilog.Print("Warning: Cannot use filenames for higher levels currently! Ignoring!");
|
music::ilog.Print("Warning: Cannot use filenames for higher levels currently! Ignoring!");
|
||||||
|
|
||||||
randc[i] = new rng(*randc[i - 1], ran_cube_size_, rngseeds_[i], true);
|
randc_[i] = new rng(*randc_[i - 1], ran_cube_size_, rngseeds_[i], true);
|
||||||
delete randc[i - 1];
|
delete randc_[i - 1];
|
||||||
randc[i - 1] = NULL;
|
randc_[i - 1] = NULL;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -188,9 +242,9 @@ void RNG_music::compute_random_numbers(void)
|
||||||
if (levelmin_seed_ > levelmin_)
|
if (levelmin_seed_ > levelmin_)
|
||||||
{
|
{
|
||||||
if (rngfnames_[levelmin_seed_].size() > 0)
|
if (rngfnames_[levelmin_seed_].size() > 0)
|
||||||
randc[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
|
randc_[levelmin_seed_] = new rng(1 << levelmin_seed_, rngfnames_[levelmin_seed_], rndsign);
|
||||||
else
|
else
|
||||||
randc[levelmin_seed_] =
|
randc_[levelmin_seed_] =
|
||||||
new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_], true); //, x0, lx );
|
new rng(1 << levelmin_seed_, ran_cube_size_, rngseeds_[levelmin_seed_], true); //, x0, lx );
|
||||||
|
|
||||||
for (int ilevel = levelmin_seed_ - 1; ilevel >= (int)levelmin_; --ilevel)
|
for (int ilevel = levelmin_seed_ - 1; ilevel >= (int)levelmin_; --ilevel)
|
||||||
|
@ -201,14 +255,14 @@ void RNG_music::compute_random_numbers(void)
|
||||||
ilevel, levelmin_seed_);
|
ilevel, levelmin_seed_);
|
||||||
|
|
||||||
// if( brealspace_tf && ilevel < levelmax_ )
|
// if( brealspace_tf && ilevel < levelmax_ )
|
||||||
// randc[ilevel] = new rng( *randc[ilevel+1], false );
|
// randc_[ilevel] = new rng( *randc_[ilevel+1], false );
|
||||||
// else // do k-space averaging
|
// else // do k-space averaging
|
||||||
randc[ilevel] = new rng(*randc[ilevel + 1], true);
|
randc_[ilevel] = new rng(*randc_[ilevel + 1], true);
|
||||||
|
|
||||||
if (ilevel + 1 > levelmax_)
|
if (ilevel + 1 > levelmax_)
|
||||||
{
|
{
|
||||||
delete randc[ilevel + 1];
|
delete randc_[ilevel + 1];
|
||||||
randc[ilevel + 1] = NULL;
|
randc_[ilevel + 1] = NULL;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -216,12 +270,12 @@ void RNG_music::compute_random_numbers(void)
|
||||||
//--- GENERATE AND STORE ALL LEVELS, INCLUDING REFINEMENTS ---//
|
//--- GENERATE AND STORE ALL LEVELS, INCLUDING REFINEMENTS ---//
|
||||||
|
|
||||||
//... levelmin
|
//... levelmin
|
||||||
if (randc[levelmin_] == NULL)
|
if (randc_[levelmin_] == NULL)
|
||||||
{
|
{
|
||||||
if (rngfnames_[levelmin_].size() > 0)
|
if (rngfnames_[levelmin_].size() > 0)
|
||||||
randc[levelmin_] = new rng(1 << levelmin_, rngfnames_[levelmin_], rndsign);
|
randc_[levelmin_] = new rng(1 << levelmin_, rngfnames_[levelmin_], rndsign);
|
||||||
else
|
else
|
||||||
randc[levelmin_] = new rng(1 << levelmin_, ran_cube_size_, rngseeds_[levelmin_], true);
|
randc_[levelmin_] = new rng(1 << levelmin_, ran_cube_size_, rngseeds_[levelmin_], true);
|
||||||
}
|
}
|
||||||
|
|
||||||
// if( levelmax_ == levelmin_ )
|
// if( levelmax_ == levelmin_ )
|
||||||
|
@ -232,11 +286,11 @@ void RNG_music::compute_random_numbers(void)
|
||||||
//... constraints.apply will return without doing anything
|
//... constraints.apply will return without doing anything
|
||||||
int x0[3] = { 0, 0, 0 };
|
int x0[3] = { 0, 0, 0 };
|
||||||
int lx[3] = { 1<<levelmin_, 1<<levelmin_, 1<<levelmin_ };
|
int lx[3] = { 1<<levelmin_, 1<<levelmin_, 1<<levelmin_ };
|
||||||
constraints.apply( levelmin_, x0, lx, randc[levelmin_] );
|
constraints.apply( levelmin_, x0, lx, randc_[levelmin_] );
|
||||||
}
|
}
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
// store_rnd(levelmin_, randc[levelmin_]);
|
// store_rnd(levelmin_, randc_[levelmin_]);
|
||||||
|
|
||||||
//... refinement levels
|
//... refinement levels
|
||||||
// for (int ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel)
|
// for (int ilevel = levelmin_ + 1; ilevel <= levelmax_; ++ilevel)
|
||||||
|
@ -258,22 +312,22 @@ void RNG_music::compute_random_numbers(void)
|
||||||
// x0[1] = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - lx[1] / 4;
|
// x0[1] = prefh_->offset_abs(ilevel, 1) - lfac * shift[1] - lx[1] / 4;
|
||||||
// x0[2] = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - lx[2] / 4;
|
// x0[2] = prefh_->offset_abs(ilevel, 2) - lfac * shift[2] - lx[2] / 4;
|
||||||
|
|
||||||
// if (randc[ilevel] == NULL)
|
// if (randc_[ilevel] == NULL)
|
||||||
// randc[ilevel] =
|
// randc_[ilevel] =
|
||||||
// new rng(*randc[ilevel - 1], ran_cube_size_, rngseeds_[ilevel], kavg, ilevel == levelmin_ + 1, x0, lx);
|
// new rng(*randc_[ilevel - 1], ran_cube_size_, rngseeds_[ilevel], kavg, ilevel == levelmin_ + 1, x0, lx);
|
||||||
// delete randc[ilevel - 1];
|
// delete randc_[ilevel - 1];
|
||||||
// randc[ilevel - 1] = NULL;
|
// randc_[ilevel - 1] = NULL;
|
||||||
|
|
||||||
// //... apply constraints to this level, if any
|
// //... apply constraints to this level, if any
|
||||||
// // if( ilevel == levelmax_ )
|
// // if( ilevel == levelmax_ )
|
||||||
// // constraints.apply( ilevel, x0, lx, randc[ilevel] );
|
// // constraints.apply( ilevel, x0, lx, randc_[ilevel] );
|
||||||
|
|
||||||
// //... store numbers
|
// //... store numbers
|
||||||
// store_rnd(ilevel, randc[ilevel]);
|
// store_rnd(ilevel, randc_[ilevel]);
|
||||||
// }
|
// }
|
||||||
|
|
||||||
delete randc[levelmax_];
|
// delete randc_[levelmax_];
|
||||||
randc[levelmax_] = NULL;
|
// randc_[levelmax_] = NULL;
|
||||||
|
|
||||||
//... make sure that the coarse grid contains oct averages where it overlaps with a fine grid
|
//... make sure that the coarse grid contains oct averages where it overlaps with a fine grid
|
||||||
//... this also ensures that constraints enforced on fine grids are carried to the coarser grids
|
//... this also ensures that constraints enforced on fine grids are carried to the coarser grids
|
||||||
|
@ -287,9 +341,9 @@ void RNG_music::compute_random_numbers(void)
|
||||||
/*if( levelmax_rand_ >= (int)levelmin_ )
|
/*if( levelmax_rand_ >= (int)levelmin_ )
|
||||||
{
|
{
|
||||||
std::cerr << "lmaxread >= (int)levelmin\n";
|
std::cerr << "lmaxread >= (int)levelmin\n";
|
||||||
randc[levelmax_rand_] = new rng( (unsigned)pow(2,levelmax_rand_), rngfnames_[levelmax_rand_] );
|
randc_[levelmax_rand_] = new rng( (unsigned)pow(2,levelmax_rand_), rngfnames_[levelmax_rand_] );
|
||||||
for( int ilevel = levelmax_rand_-1; ilevel >= (int)levelmin_; --ilevel )
|
for( int ilevel = levelmax_rand_-1; ilevel >= (int)levelmin_; --ilevel )
|
||||||
randc[ilevel] = new rng( *randc[ilevel+1] );
|
randc_[ilevel] = new rng( *randc_[ilevel+1] );
|
||||||
}*/
|
}*/
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -297,5 +351,5 @@ void RNG_music::compute_random_numbers(void)
|
||||||
|
|
||||||
namespace
|
namespace
|
||||||
{
|
{
|
||||||
RNG_plugin_creator_concrete<RNG_music> creator("MUSIC");
|
RNG_plugin_creator_concrete<RNG_music> creator("MUSIC1");
|
||||||
}
|
}
|
|
@ -41,38 +41,8 @@ music_wnoise_generator<T>::music_wnoise_generator(unsigned res, unsigned cubesiz
|
||||||
double mean = 0.0;
|
double mean = 0.0;
|
||||||
size_t res_l = res;
|
size_t res_l = res;
|
||||||
|
|
||||||
bool musicnoise = true;
|
|
||||||
if (!musicnoise)
|
|
||||||
cubesize_ = res_;
|
|
||||||
|
|
||||||
if (!musicnoise)
|
|
||||||
music::elog.Print("This currently breaks compatibility. Need to disable by hand! Make sure to not check into repo");
|
|
||||||
|
|
||||||
initialize();
|
initialize();
|
||||||
|
mean = fill_all();
|
||||||
if (musicnoise)
|
|
||||||
mean = fill_all();
|
|
||||||
else
|
|
||||||
{
|
|
||||||
rnums_.push_back(new Meshvar<T>(res, 0, 0, 0));
|
|
||||||
cubemap_[0] = 0; // create dummy map index
|
|
||||||
register_cube(0, 0, 0);
|
|
||||||
//rapid_proto_ngenic_rng( res_, baseseed_, *this );
|
|
||||||
}
|
|
||||||
|
|
||||||
/*
|
|
||||||
|
|
||||||
if( musicnoise )
|
|
||||||
mean = fill_all();
|
|
||||||
else
|
|
||||||
{
|
|
||||||
rnums_.push_back( new Meshvar<T>( res, 0, 0, 0 ) );
|
|
||||||
cubemap_[0] = 0; // create dummy map index
|
|
||||||
register_cube(0,0,0);
|
|
||||||
rapid_proto_ngenic_rng( res_, baseseed_, *this );
|
|
||||||
}
|
|
||||||
|
|
||||||
*/
|
|
||||||
|
|
||||||
if (zeromean)
|
if (zeromean)
|
||||||
{
|
{
|
||||||
|
|
|
@ -117,7 +117,7 @@ public:
|
||||||
if (i == nres_ / 2 || j == nres_ / 2 || k == nres_ / 2) continue;
|
if (i == nres_ / 2 || j == nres_ / 2 || k == nres_ / 2) continue;
|
||||||
if (i == 0 && j == 0 && k == 0) continue;
|
if (i == 0 && j == 0 && k == 0) continue;
|
||||||
|
|
||||||
ampl = -std::sqrt(-std::log(ampl));
|
ampl = std::sqrt(-std::log(ampl));
|
||||||
ccomplex_t zrand(ampl*std::cos(phase),ampl*std::sin(phase));
|
ccomplex_t zrand(ampl*std::cos(phase),ampl*std::sin(phase));
|
||||||
|
|
||||||
if (k > 0) {
|
if (k > 0) {
|
||||||
|
@ -143,7 +143,7 @@ public:
|
||||||
}
|
}
|
||||||
} else if (i < nres_ / 2) {
|
} else if (i < nres_ / 2) {
|
||||||
if(i_in_range) g.kelem(ip,j,k) = zrand;
|
if(i_in_range) g.kelem(ip,j,k) = zrand;
|
||||||
if (ii_in_range) g.kelem(iip,jj,k) = std::conj(zrand);
|
if(ii_in_range) g.kelem(iip,jj,k) = std::conj(zrand);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
|
@ -1,6 +1,7 @@
|
||||||
// This file is part of monofonIC (MUSIC2)
|
// This file is part of monofonIC (MUSIC2)
|
||||||
// A software package to generate ICs for cosmological simulations
|
// A software package to generate ICs for cosmological simulations
|
||||||
// Copyright (C) 2020 by Oliver Hahn and Adrian Jenkins (this file)
|
// Copyright (C) 2020 by Oliver Hahn and Adrian Jenkins (this file)
|
||||||
|
// but see distinct licensing for PANPHASIA below
|
||||||
//
|
//
|
||||||
// monofonIC is free software: you can redistribute it and/or modify
|
// monofonIC is free software: you can redistribute it and/or modify
|
||||||
// it under the terms of the GNU General Public License as published by
|
// it under the terms of the GNU General Public License as published by
|
||||||
|
@ -14,6 +15,13 @@
|
||||||
//
|
//
|
||||||
// You should have received a copy of the GNU General Public License
|
// You should have received a copy of the GNU General Public License
|
||||||
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
// along with this program. If not, see <http://www.gnu.org/licenses/>.
|
||||||
|
//
|
||||||
|
// IMPORTANT NOTICE:
|
||||||
|
// Note that PANPHASIA itself is not released under the GPL. Make sure
|
||||||
|
// to read and agree to its distinct licensing before you use or modify
|
||||||
|
// the code below or in the /external/panphasia directory which can be
|
||||||
|
// found here: http://icc.dur.ac.uk/Panphasia.php
|
||||||
|
// NOTE THAT PANPHASIA REQUIRES REGISTRATION ON THIS WEBSITE PRIOR TO USE
|
||||||
|
|
||||||
#if defined(USE_PANPHASIA)
|
#if defined(USE_PANPHASIA)
|
||||||
|
|
||||||
|
@ -282,14 +290,15 @@ public:
|
||||||
{
|
{
|
||||||
panphasia_descriptor d(descriptor_string_);
|
panphasia_descriptor d(descriptor_string_);
|
||||||
|
|
||||||
int lextra = (log10((double)ngrid_ / (double)d.i_base) + 0.001) / log10(2.0);
|
const int lextra = (log10((double)ngrid_ / (double)d.i_base) + 0.001) / log10(2.0);
|
||||||
int level_p = d.wn_level_base + lextra;
|
int level_p = d.wn_level_base + lextra;
|
||||||
int ratio = 1 << lextra;
|
const int ratio = 1 << lextra;
|
||||||
|
|
||||||
lstate[mythread].layer_min = 0;
|
lstate[mythread].layer_min = 0;
|
||||||
lstate[mythread].layer_max = level_p;
|
lstate[mythread].layer_max = level_p;
|
||||||
lstate[mythread].indep_field = 1;
|
lstate[mythread].indep_field = 1;
|
||||||
|
|
||||||
|
_unused(ratio);
|
||||||
assert(ngrid_ == ratio * d.i_base);
|
assert(ngrid_ == ratio * d.i_base);
|
||||||
|
|
||||||
long long ix_rel[3];
|
long long ix_rel[3];
|
||||||
|
@ -419,14 +428,15 @@ public:
|
||||||
{
|
{
|
||||||
panphasia_descriptor d(descriptor_string_);
|
panphasia_descriptor d(descriptor_string_);
|
||||||
|
|
||||||
int lextra = (log10((double)ngrid_ / (double)d.i_base) + 0.001) / log10(2.0);
|
const int lextra = (log10((double)ngrid_ / (double)d.i_base) + 0.001) / log10(2.0);
|
||||||
int level_p = d.wn_level_base + lextra;
|
int level_p = d.wn_level_base + lextra;
|
||||||
int ratio = 1 << lextra;
|
const int ratio = 1 << lextra;
|
||||||
|
|
||||||
lstate[mythread].layer_min = 0;
|
lstate[mythread].layer_min = 0;
|
||||||
lstate[mythread].layer_max = level_p;
|
lstate[mythread].layer_max = level_p;
|
||||||
lstate[mythread].indep_field = 1;
|
lstate[mythread].indep_field = 1;
|
||||||
|
|
||||||
|
_unused(ratio);
|
||||||
assert(ngrid_ == ratio * d.i_base);
|
assert(ngrid_ == ratio * d.i_base);
|
||||||
|
|
||||||
long long ix_rel[3];
|
long long ix_rel[3];
|
||||||
|
|
|
@ -184,6 +184,8 @@ public:
|
||||||
transfer_CAMB_file_plugin(config_file &cf)
|
transfer_CAMB_file_plugin(config_file &cf)
|
||||||
: TransferFunction_plugin(cf)
|
: TransferFunction_plugin(cf)
|
||||||
{
|
{
|
||||||
|
music::wlog << "The CAMB file plugin is not well tested! Proceed with checks of correctness of output before running a simulation!" << std::endl;
|
||||||
|
|
||||||
m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file");
|
m_filename_Tk = pcf_->get_value<std::string>("cosmology", "transfer_file");
|
||||||
m_Omega_m = cf.get_value<double>("cosmology", "Omega_m"); //MvD
|
m_Omega_m = cf.get_value<double>("cosmology", "Omega_m"); //MvD
|
||||||
m_Omega_b = cf.get_value<double>("cosmology", "Omega_b"); //MvD
|
m_Omega_b = cf.get_value<double>("cosmology", "Omega_b"); //MvD
|
||||||
|
@ -249,6 +251,7 @@ public:
|
||||||
double lk = log10(k);
|
double lk = log10(k);
|
||||||
double dk = m_tab_k[n] - m_tab_k[n1];
|
double dk = m_tab_k[n] - m_tab_k[n1];
|
||||||
double delk = lk - m_tab_k[n];
|
double delk = lk - m_tab_k[n];
|
||||||
|
double dc{0.0}, db{0.0};
|
||||||
|
|
||||||
switch (type)
|
switch (type)
|
||||||
{
|
{
|
||||||
|
@ -280,7 +283,14 @@ public:
|
||||||
v1 = m_tab_Tk_tot[n1];
|
v1 = m_tab_Tk_tot[n1];
|
||||||
v2 = m_tab_Tk_tot[n];
|
v2 = m_tab_Tk_tot[n];
|
||||||
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
return pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||||
#warning TODO: add deltabc here
|
case deltabc:
|
||||||
|
v1 = m_tab_Tk_cdm[n1];
|
||||||
|
v2 = m_tab_Tk_cdm[n];
|
||||||
|
dc = pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||||
|
v1 = m_tab_Tk_baryon[n1];
|
||||||
|
v2 = m_tab_Tk_baryon[n];
|
||||||
|
db = pow(10.0, (v2 - v1) / dk * (delk) + v2);
|
||||||
|
return db-dc;
|
||||||
default:
|
default:
|
||||||
throw std::runtime_error(
|
throw std::runtime_error(
|
||||||
"Invalid type requested in transfer function evaluation");
|
"Invalid type requested in transfer function evaluation");
|
||||||
|
|
Loading…
Reference in a new issue