mirror/monofonIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00
Code Releases Activity

Merge branch 'master' of bitbucket.org:ohahn/monofonic into develop

Browse source
This commit is contained in:
Oliver Hahn 2021-05-04 22:19:09 +02:00
commit 73c0957b66
4 changed files with 118 additions and 46 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
example.conf
View file

@ -98,6 +98,7 @@ seed = 12345
# generator = PANPHASIA
# descriptor = [Panph1,L10,(800,224,576),S9,CH1564365824,MXXL]
# PanphasiaMinRootResolution = 512 # requires the white noise reallisation to be made at least at that resolution (default is 512)
##> 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)

55
src/plugins/output_genericio.cc
View file

@ -10,7 +10,7 @@ protected:
real_t lunit_, vunit_, munit_;
bool hacc_hydro_;
float hacc_etamax_;
float hh_value_, rho_value_, mu_value_;
float hh_value_, rho_value_, mu_value_, uu_value_;
std::vector<int64_t> ids;
std::vector<float> xx, yy, zz;
std::vector<float> vx, vy, vz;
@ -20,22 +20,43 @@ protected:
public:
//! constructor
explicit genericio_output_plugin(config_file &cf, cosmology::calculator &cc)
: output_plugin(cf, cc, "GenericIO")
explicit genericio_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
: output_plugin(cf, pcc, "GenericIO")
{
real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
const double astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
const double hubble_param = pcc->cosmo_param_["h"];
lunit_ = cf_.get_value<double>("setup", "BoxLength");
vunit_ = lunit_/astart;
hacc_hydro_ = cf_.get_value_safe<bool>("output", "GenericIO_HACCHydro", false);
// initial smoothing length is mean particle seperation
hh_value_ = lunit_ / cf_.get_value<float>("setup", "GridRes");
// neutral value for molecular weight. FIXME: account for ionization?
const float primordial_x = 0.75;
mu_value_ = 4.0 / (1.0 + 3.0*primordial_x);
double rhoc = 27.7519737 * 1e10; // in h^2 M_sol / Mpc^3
rho_value_ = cc_.cosmo_param_["Omega_b"] * rhoc;
munit_ = rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3);
if(hacc_hydro_) {
const double omegab = pcc_->cosmo_param_["Omega_b"];
const double gamma = cf_.get_value_safe<double>("cosmology", "gamma", 5.0 / 3.0);
const double YHe = pcc_->cosmo_param_["YHe"];
const double Tcmb0 = pcc_->cosmo_param_["Tcmb"];
// compute gas internal energy
const double npol = (fabs(1.0 - gamma) > 1e-7) ? 1.0 / (gamma - 1.) : 1.0;
const double unitv = 1e5; // km/s to cm/s
const double adec = 1.0 / (160. * std::pow(omegab * hubble_param * hubble_param / 0.022, 2.0 / 5.0));
const double Tini = astart < adec ? Tcmb0 / astart : Tcmb0 / astart / astart * adec;
const double mu = (Tini > 1.e4) ? 4.0 / (8. - 5. * YHe) : 4.0 / (1. + 3. * (1. - YHe));
const double ceint = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / unitv / unitv;
uu_value_ = ceint/astart/astart; // probably needs a scale factor correction (might be in physical units).
mu_value_ = mu;
music::ilog.Print("HACC : calculated redshift of decoupling: z_dec = %.2f", 1./adec - 1.);
music::ilog.Print("HACC : set initial gas temperature to %.2e K/mu", Tini / mu);
music::ilog.Print("HACC : set initial internal energy to %.2e km^2/s^2", ceint);
rho_value_ = omegab * rhoc;
}
}
output_type write_species_as(const cosmo_species &) const { return output_type::particles; }
@ -65,17 +86,15 @@ public:
vz.reserve(vz.size() + npart);
ids.reserve(ids.size() + npart);
auto _pos = reinterpret_cast<const float*>(pc.get_pos32_ptr());
auto _vel = reinterpret_cast<const float*>(pc.get_vel32_ptr());
auto _ids = reinterpret_cast<const uint64_t*>(pc.get_ids64_ptr());
auto _mass = reinterpret_cast<const float*>(pc.get_mass32_ptr());
for(size_t i=0; i<npart; ++i) {
xx.push_back(_pos[3*i + 0]);
yy.push_back(_pos[3*i + 1]);
zz.push_back(_pos[3*i + 2]);
xx.push_back(fmod(_pos[3*i + 0]+lunit_, lunit_));
yy.push_back(fmod(_pos[3*i + 1]+lunit_, lunit_));
zz.push_back(fmod(_pos[3*i + 2]+lunit_, lunit_));
vx.push_back(_vel[3*i + 0]);
vy.push_back(_vel[3*i + 1]);
vz.push_back(_vel[3*i + 2]);
@ -100,13 +119,13 @@ public:
std::fill(mask.begin() + prev_size, mask.end(), s == cosmo_species::baryon ? 1<<2 : 0);
hh.resize(new_size);
std::fill(hh.begin() + prev_size, hh.end(), s == cosmo_species::baryon ? hh_value_ : 0.0f);
std::fill(hh.begin() + prev_size, hh.end(), hh_value_);
uu.resize(new_size);
std::fill(uu.begin() + prev_size, uu.end(), s == cosmo_species::baryon ? 0.0f : 0.0f);
std::fill(uu.begin() + prev_size, uu.end(), s == cosmo_species::baryon ? uu_value_ : 0.0f);
rho.resize(new_size);
std::fill(rho.begin() + prev_size, rho.end(), s == cosmo_species::baryon ? rho_value_ : 0.0f);
std::fill(rho.begin() + prev_size, rho.end(), rho_value_);
mu.resize(new_size);
std::fill(mu.begin() + prev_size, mu.end(), s == cosmo_species::baryon ? mu_value_ : 0.0f);
@ -141,7 +160,7 @@ public:
namespace
{
output_plugin_creator_concrete<genericio_output_plugin> creator("genericio");
output_plugin_creator_concrete<genericio_output_plugin> creator1("genericio");
} // namespace
#endif // ENABLE_GENERICIO

97
src/plugins/output_swift.cc
View file

@ -22,9 +22,9 @@
#include "HDF_IO.hh"
template <typename T>
std::vector<T> from_6array(const std::array<T,6>& a)
std::vector<T> from_7array(const std::array<T,7>& a)
{
return std::vector<T>{{a[0], a[1], a[2], a[3], a[4], a[5]}};
return std::vector<T>{{a[0], a[1], a[2], a[3], a[4], a[5], a[6]}};
}
template <typename T>
@ -41,12 +41,12 @@ protected:
int num_files_, num_simultaneous_writers_;
real_t lunit_, vunit_, munit_;
real_t boxsize_, hubble_param_, astart_;
real_t boxsize_, hubble_param_, astart_, zstart_;
bool blongids_, bdobaryons_;
std::string this_fname_;
std::array<uint32_t,6> npart_, npartTotal_, npartTotalHighWord_;
std::array<double,6> mass_;
std::array<uint32_t,7> npart_, npartTotal_, npartTotalHighWord_;
std::array<double,7> mass_;
double time_;
public:
@ -63,7 +63,8 @@ public:
const double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
hubble_param_ = pcc->cosmo_param_["h"];
astart_ = 1.0/(1.0+cf_.get_value<double>("setup","zstart"));
zstart_ = cf_.get_value<double>("setup","zstart");
astart_ = 1.0/(1.0 + zstart_);
boxsize_ = cf_.get_value<double>("setup", "BoxLength");
lunit_ = boxsize_ / hubble_param_; // final units will be in Mpc (without h)
@ -73,7 +74,7 @@ public:
blongids_ = cf_.get_value_safe<bool>("output", "UseLongids", false);
bdobaryons_ = cf_.get_value<bool>("setup","DoBaryons");
for (int i = 0; i < 6; ++i)
for (int i = 0; i < 7; ++i)
{
npart_[i] = 0;
npartTotal_[i] = 0;
@ -91,23 +92,63 @@ public:
this_fname_ += "." + std::to_string(thisrank);
#endif
// Now that we modified the filename for MPI support, propagate the change to the super-class
// as it does operations on the file behind our back.
fname_ = this_fname_;
// delete output file if it exists
unlink(this_fname_.c_str());
// create output HDF5 file
HDFCreateFile(this_fname_);
// Write UNITS header
// Write UNITS header using the physical constants assumed internally by SWIFT
HDFCreateGroup(this_fname_, "Units");
HDFWriteGroupAttribute(this_fname_, "Units", "Unit mass in cgs (U_M)", 1.98848e43); // 10^10 Msun in grams
HDFWriteGroupAttribute(this_fname_, "Units", "Unit length in cgs (U_L)", 3.08567758e24); // 1 Mpc in cm
HDFWriteGroupAttribute(this_fname_, "Units", "Unit time in cgs (U_t)", 3.08567758e19); // so that unit vel is 1 km/s
HDFWriteGroupAttribute(this_fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere
HDFWriteGroupAttribute(this_fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin
HDFWriteGroupAttribute(this_fname_, "Units", "Unit mass in cgs (U_M)", 1.98841e43); // 10^10 Msun in grams
HDFWriteGroupAttribute(this_fname_, "Units", "Unit length in cgs (U_L)", 3.08567758149e24); // 1 Mpc in cm
HDFWriteGroupAttribute(this_fname_, "Units", "Unit time in cgs (U_t)", 3.08567758149e19); // so that unit vel is 1 km/s
HDFWriteGroupAttribute(this_fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere
HDFWriteGroupAttribute(this_fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin
// TODO: Write MUSIC configuration header
HDFCreateGroup(fname_, "ICs_parameters");
// ...
// Write MUSIC configuration header
int order = cf_.get_value<int>("setup", "LPTorder");
std::string load = cf_.get_value<std::string>("setup", "ParticleLoad");
std::string tf = cf_.get_value<std::string>("cosmology", "transfer");
std::string cosmo_set = cf_.get_value<std::string>("cosmology", "ParameterSet");
std::string rng = cf_.get_value<std::string>("random", "generator");
int do_fixing = cf_.get_value<bool>("setup", "DoFixing");
int do_invert = cf_.get_value<bool>("setup", "DoInversion");
int do_baryons = cf_.get_value<bool>("setup", "DoBaryons");
int do_baryonsVrel = cf_.get_value<bool>("setup", "DoBaryonVrel");
int L = cf_.get_value<int>("setup", "GridRes");
HDFCreateGroup(this_fname_, "ICs_parameters");
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Code", std::string("MUSIC2 - monofonIC"));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Revision", std::string(GIT_REV));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Tag", std::string(GIT_TAG));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Branch", std::string(GIT_BRANCH));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Precision", std::string(CMAKE_PRECISION_STR));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Convolutions", std::string(CMAKE_CONVOLVER_STR));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "PLT", std::string(CMAKE_PLT_STR));
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "LPT Order", order);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Particle Load", load);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Transfer Function", tf);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Cosmology Parameter Set", cosmo_set);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Random Generator", rng);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Mode Fixing", do_fixing);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Mode inversion", do_invert);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Baryons", do_baryons);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Baryons Relative Velocity", do_baryonsVrel);
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Grid Resolution", L);
if (tf == "CLASS") {
double ztarget = cf_.get_value<double>("cosmology", "ztarget");
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Target Redshift", ztarget);
}
if (rng == "PANPHASIA") {
std::string desc = cf_.get_value<std::string>("random", "descriptor");
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Descriptor", desc);
}
}
// use destructor to write header post factum
@ -117,19 +158,20 @@ public:
{
HDFCreateGroup(this_fname_, "Header");
HDFWriteGroupAttribute(fname_, "Header", "Dimension", 3);
HDFWriteGroupAttribute(fname_, "Header", "BoxSize", lunit_); // in Mpc, not Mpc/h
HDFWriteGroupAttribute(fname_, "Header", "BoxSize", boxsize_ / hubble_param_); // in Mpc, not Mpc/h
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_6array<unsigned>(npartTotal_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_6array<unsigned>(npartTotalHighWord_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_6array<unsigned>(npart_));
HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_6array<double>(mass_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_7array<unsigned>(npartTotal_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_7array<unsigned>(npartTotalHighWord_));
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_7array<unsigned>(npart_));
HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_7array<double>(mass_));
HDFWriteGroupAttribute(this_fname_, "Header", "Time", from_value<double>(time_));
HDFWriteGroupAttribute(this_fname_, "Header", "Redshift", from_value<double>(zstart_));
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Entropy_ICs", from_value<int>(0));
HDFWriteGroupAttribute(this_fname_, "Header", "NumFilesPerSnapshot", from_value<int>(num_files_));
// write GAS internal energy if baryons are enabled
// write GAS internal energy and smoothing length if baryons are enabled
if( bdobaryons_ )
{
const double gamma = cf_.get_value_safe<double>("cosmology", "gamma", 5.0 / 3.0);
@ -147,12 +189,17 @@ public:
const double mu = (Tini > 1.e4) ? 4.0 / (8. - 5. * YHe) : 4.0 / (1. + 3. * (1. - YHe));
const double ceint = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / unitv / unitv;
music::ilog.Print("Swift : set initial gas temperature to %.2f K/mu", Tini / mu);
music::ilog.Print("Swift : Calculated initial gas temperature: %.2f K/mu", Tini / mu);
music::ilog.Print("Swift : set initial internal energy to %.2e km^2/s^2", ceint);
std::vector<write_real_t> data( npart_[0], ceint );
HDFWriteDataset(this_fname_, "PartType0/InternalEnergy", data);
data.assign( npart_[0], boxsize_ / cf_.get_value<double>("setup","GridRes") );
const double h = boxsize_ / hubble_param_ / cf_.get_value<double>("setup","GridRes");
music::ilog.Print("Swift : set initial smoothing length to mean inter-part separation: %.2f Mpc", h);
data.assign( npart_[0], h);
HDFWriteDataset(this_fname_, "PartType0/SmoothingLength", data);
}
@ -192,7 +239,7 @@ public:
case cosmo_species::baryon:
return 0;
case cosmo_species::neutrino:
return 3;
return 6;
}
return -1;
}

11
src/plugins/random_panphasia.cc
View file

@ -171,22 +171,27 @@ protected:
{
// if ngrid is not a multiple of i_base, then we need to enlarge and then sample down
ngrid_ = pcf_->get_value<size_t>("setup", "GridRes");
int ngridminsize_panphasia = pcf_->get_value_safe<int>("random", "PanphasiaMinRootResolution",512);
grid_p_ = pdescriptor_->i_base;
lextra_ = (log10((double)ngrid_ / (double)grid_p_) + 0.001) / log10(2.0);
// lmin
ngrid_panphasia_ = (1 << lextra_) * grid_p_;
if( ngrid_panphasia_ < ngrid_ ){
while( ngrid_panphasia_ < ngridminsize_panphasia ){
lextra_++;
ngrid_panphasia_*=2;
}
assert( ngrid_panphasia_ >= ngrid_ );
assert( ngrid_panphasia_ >= ngridminsize_panphasia);
clear_panphasia_thread_states();
music::ilog.Print("PANPHASIA: using grid size %lld (level=%d)",ngrid_panphasia_, lextra_);
if (ngridminsize_panphasia<512)
music::ilog.Print("PANPHASIA WARNING: PanphasiaMinRootResolution = %d below minimum recommended of 512",ngridminsize_panphasia);
music::ilog.Print("PANPHASIA: running with %d threads", num_threads_, ngrid_panphasia_ );
coordinate_system_shift_[0] = -pcf_->get_value_safe<int>("setup", "shift_x", 0);