mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
fix genericio plugin
This commit is contained in:
parent
0c4cb7e5ae
commit
c28dc782e3
1 changed files with 37 additions and 18 deletions
|
@ -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
|
Loading…
Reference in a new issue