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

Merged in master (pull request #19)

Browse source 
fix genericio plugin

Approved-by: Oliver Hahn
This commit is contained in:
Michael Bühlmann 2020-12-29 09:09:12 +00:00 committed by Oliver Hahn
commit 1ef2b417e3
1 changed files with 37 additions and 18 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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