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 'class_update' of bitbucket.org:ohahn/monofonic into generic_cosmo

Browse source
This commit is contained in:
Oliver Hahn 2023-04-04 18:32:11 +02:00
commit 20bd8697ac
26 changed files with 1182 additions and 218 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
CMakeLists.txt
View file

@ -263,7 +263,7 @@ endif(ENABLE_PANPHASIA)
if(ENABLE_PLT)
target_compile_definitions(${PRGNAME} PRIVATE "ENABLE_PLT")
endif(ENABLE_PLT)
target_include_directories(${PRGNAME} PRIVATE ${GSL_INCLUDE_DIR})
target_link_libraries(${PRGNAME} PRIVATE GSL::gsl)

6
example.conf
View file

@ -20,7 +20,11 @@ DoFixing = no # do mode fixing à la Angulo&Pontzen (https://arxiv.o
DoInversion = no # invert phases (for paired simulations)
ParticleLoad = sc # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x)
# (increases number of particles by given factor!), or 'glass'
# (increases number of particles by given factor!),
# or 'glass' or 'masked'
## if `ParticleLoad = masked' then you can specify here how masking should take place
# ParticleMaskType = 3 # bit mask for particle mask (0=center,1=center+edges,2=center+faces,3=center+edges+faces)
## if `ParticleLoad = glass' then specify here where to load the glass distribution from
# GlassFileName = glass128.hdf5

18
external/panphasia_ho/uniform_rand_threefry4x64.c vendored
View file

@ -100,8 +100,8 @@ void threefry4x64_test_(int verbose)
if ((rand.v[0] != result.v[0]) || (rand.v[1] != result.v[1]) || (rand.v[2] != result.v[2]) || (rand.v[3] != result.v[3]))
{
printf("Serious error occured !!!!!!!!!! Random generator is not working correctly \n");
printf("Random generated: %lu %lu %lu %lu\n", rand.v[0], rand.v[1], rand.v[2], rand.v[3]);
printf("Random expected: %lu %lu %lu %lu\n", result.v[0], result.v[1], result.v[2], result.v[3]);
printf("Random generated: %llu %llu %llu %llu\n", rand.v[0], rand.v[1], rand.v[2], rand.v[3]);
printf("Random expected: %llu %llu %llu %llu\n", result.v[0], result.v[1], result.v[2], result.v[3]);
//abort();
}
else
@ -121,8 +121,8 @@ void threefry4x64_test_(int verbose)
if ((rand.v[0] != result.v[0]) || (rand.v[1] != result.v[1]) || (rand.v[2] != result.v[2]) || (rand.v[3] != result.v[3]))
{
printf("Serious error occured !!!!!!!!!! Random generator is not working correctly \n");
printf("Random generated: %lu %lu %lu %lu\n", rand.v[0], rand.v[1], rand.v[2], rand.v[3]);
printf("Random expected: %lu %lu %lu %lu\n", result.v[0], result.v[1], result.v[2], result.v[3]);
printf("Random generated: %llu %llu %llu %llu\n", rand.v[0], rand.v[1], rand.v[2], rand.v[3]);
printf("Random expected: %llu %llu %llu %llu\n", result.v[0], result.v[1], result.v[2], result.v[3]);
//abort();
}
else
@ -143,8 +143,8 @@ void threefry4x64_test_(int verbose)
if ((rand.v[0] != result.v[0]) || (rand.v[1] != result.v[1]) || (rand.v[2] != result.v[2]) || (rand.v[3] != result.v[3]))
{
printf("Serious error occured !!!!!!!!!! Random generator is not working correctly \n");
printf("Random generated: %lu %lu %lu %lu\n", rand.v[0], rand.v[1], rand.v[2], rand.v[3]);
printf("Random expected: %lu %lu %lu %lu\n", result.v[0], result.v[1], result.v[2], result.v[3]);
printf("Random generated: %llu %llu %llu %llu\n", rand.v[0], rand.v[1], rand.v[2], rand.v[3]);
printf("Random expected: %llu %llu %llu %llu\n", result.v[0], result.v[1], result.v[2], result.v[3]);
//abort();
}
else
@ -176,7 +176,7 @@ void set_panphasia_key_(int verbose)
verbose = 0; //ARJ
if (verbose)
printf("Setting the threefry4x64 key to\n(%0lu %0lu %0lu %0lu)\n\n",
printf("Setting the threefry4x64 key to\n(%0llu %0llu %0llu %0llu)\n\n",
panphasia_key.v[0], panphasia_key.v[1], panphasia_key.v[2], panphasia_key.v[3]);
panphasia_key_initialised = 999;
@ -237,10 +237,10 @@ void check_panphasia_key_(int verbose)
if (panphasia_check_key.v[0] != panphasia_key.v[0] || panphasia_check_key.v[1] != panphasia_key.v[1] || panphasia_check_key.v[2] != panphasia_key.v[2] || panphasia_check_key.v[2] != panphasia_key.v[2])
{
printf("A serious error has happened - the threefry4x64 key has become corrupted!\n");
printf("Should be: (%0lu %0lu %0lu %0lu)\n", panphasia_check_key.v[0],
printf("Should be: (%0llu %0llu %0llu %0llu)\n", panphasia_check_key.v[0],
panphasia_check_key.v[1], panphasia_check_key.v[2], panphasia_check_key.v[3]);
printf("But now is: (%0lu %0lu %0lu %0lu)\n", panphasia_key.v[0],
printf("But now is: (%0llu %0llu %0llu %0llu)\n", panphasia_key.v[0],
panphasia_key.v[1], panphasia_key.v[2], panphasia_key.v[3]);
printf("The fact that it has changed suggests the key has been overwritten in memory.\n");
abort();

2
include/HDF_IO.hh
View file

@ -226,7 +226,7 @@ template<typename T >
inline void HDFReadSelect( const std::string Filename, const std::string ObjName, const std::vector<unsigned>& ii, std::vector<T> &Data ){
hid_t HDF_Type, HDF_FileID, HDF_DatasetID, HDF_DataspaceID, HDF_MemspaceID;
hsize_t HDF_StorageSize;
// hsize_t HDF_StorageSize;
HDF_Type = GetDataType<T>();

165
include/config_file.hh
View file

@ -1,17 +1,17 @@
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn
//
//
// monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//
// monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
@ -25,7 +25,7 @@
#include <stdexcept>
#include <string>
#include <typeinfo>
#include <vector>
#include <logger.hh>
/*!
@ -36,7 +36,8 @@
* configuration is stored in hash-pairs and can be queried and
* validated by the responsible class/routine
*/
class config_file {
class config_file
{
//! current line number
unsigned iline_;
@ -52,7 +53,8 @@ public:
* @return trimmed string
*/
std::string trim(std::string const &source,
char const *delims = " \t\r\n") const {
char const *delims = " \t\r\n") const
{
std::string result(source);
//... skip initial whitespace ...
std::string::size_type index = result.find_last_not_of(delims);
@ -76,27 +78,58 @@ public:
* @param oval the interpreted/converted value
*/
template <class in_value, class out_value>
void convert(const in_value &ival, out_value &oval) const {
void convert(const in_value &ival, out_value &oval) const
{
std::stringstream ss;
ss << ival; //.. insert value into stream
ss >> oval; //.. retrieve value from stream
if (!ss.eof()) {
if (!ss.eof())
{
//.. conversion error
music::elog << "Error: conversion of \'" << ival << "\' failed."
<< std::endl;
<< std::endl;
throw except_invalid_conversion(std::string("invalid conversion to ") +
typeid(out_value).name() + '.');
typeid(out_value).name() + '.');
}
}
template <class in_value, typename T>
void convert(const in_value &ival, std::vector<T> &oval) const
{
std::stringstream ss;
ss << ival; //.. insert value into stream
while (ss.good())
{
std::stringstream ss2;
std::string substr;
std::getline(ss, substr, ',');
T oval_single;
ss2 << substr;
ss2 >> oval_single;
oval.push_back(oval_single);
}
if (!ss.eof())
{
//.. conversion error
music::elog << "Error: conversion of \'" << ival << "\' failed."
<< std::endl;
throw except_invalid_conversion(std::string("invalid conversion to ") +
typeid(T).name() + '.');
}
}
//! constructor of class config_file
/*! @param filename the path/name of the configuration file to be parsed
*/
explicit config_file(std::string const &filename) : iline_(0), items_() {
explicit config_file(std::string const &filename) : iline_(0), items_()
{
std::ifstream file(filename.c_str());
if (!file.is_open()){
if (!file.is_open())
{
music::elog << "Could not open config file \'" << filename << "\'." << std::endl;
throw std::runtime_error(
std::string("Error: Could not open config file \'") + filename +
@ -110,7 +143,8 @@ public:
int pos_equal;
iline_ = 0;
//.. walk through all lines ..
while (std::getline(file, line)) {
while (std::getline(file, line))
{
++iline_;
//.. encounterd EOL ?
if (!line.length())
@ -122,7 +156,8 @@ public:
line.erase(idx);
//.. encountered section tag ?
if (line[0] == '[') {
if (line[0] == '[')
{
in_section = trim(line.substr(1, line.find(']') - 1));
continue;
}
@ -133,21 +168,24 @@ public:
value = trim(line.substr(pos_equal + 1));
if ((size_t)pos_equal == std::string::npos &&
(name.size() != 0 || value.size() != 0)) {
(name.size() != 0 || value.size() != 0))
{
music::wlog << "Ignoring non-assignment in " << filename << ":"
<< iline_ << std::endl;
<< iline_ << std::endl;
continue;
}
if (name.length() == 0 && value.size() != 0) {
if (name.length() == 0 && value.size() != 0)
{
music::wlog << "Ignoring assignment missing entry name in "
<< filename << ":" << iline_ << std::endl;
<< filename << ":" << iline_ << std::endl;
continue;
}
if (value.length() == 0 && name.size() != 0) {
if (value.length() == 0 && name.size() != 0)
{
music::wlog << "Empty entry will be ignored in " << filename << ":"
<< iline_ << std::endl;
<< iline_ << std::endl;
continue;
}
@ -155,9 +193,10 @@ public:
continue;
//.. add key/value pair to hash table ..
if (items_.find(in_section + '/' + name) != items_.end()) {
if (items_.find(in_section + '/' + name) != items_.end())
{
music::wlog << "Redeclaration overwrites previous value in "
<< filename << ":" << iline_ << std::endl;
<< filename << ":" << iline_ << std::endl;
}
items_[in_section + '/' + name] = value;
@ -168,7 +207,8 @@ public:
/*! @param key the key value, usually "section/key"
* @param value the value of the key, also a string
*/
void insert_value(std::string const &key, std::string const &value) {
void insert_value(std::string const &key, std::string const &value)
{
items_[key] = value;
}
@ -178,7 +218,8 @@ public:
* @param value the value of the key, also a string
*/
void insert_value(std::string const &section, std::string const &key,
std::string const &value) {
std::string const &value)
{
items_[section + '/' + key] = value;
}
@ -187,7 +228,8 @@ public:
* @param key the key name to be checked
* @return true if the key is present, false otherwise
*/
bool contains_key(std::string const &section, std::string const &key) {
bool contains_key(std::string const &section, std::string const &key)
{
std::map<std::string, std::string>::const_iterator i =
items_.find(section + '/' + key);
if (i == items_.end())
@ -199,7 +241,8 @@ public:
/*! @param key the key name to be checked
* @return true if the key is present, false otherwise
*/
bool contains_key(std::string const &key) {
bool contains_key(std::string const &key)
{
std::map<std::string, std::string>::const_iterator i = items_.find(key);
if (i == items_.end())
return false;
@ -213,7 +256,9 @@ public:
* @return the value of the key
* @sa except_item_not_found
*/
template <class T> T get_value(std::string const &key) const {
template <class T>
T get_value(std::string const &key) const
{
return get_value<T>("", key);
}
@ -226,13 +271,15 @@ public:
* @sa except_item_not_found
*/
template <class T>
T get_value_basic(std::string const &section, std::string const &key) const {
T get_value_basic(std::string const &section, std::string const &key) const
{
T r;
std::map<std::string, std::string>::const_iterator i =
items_.find(section + '/' + key);
if (i == items_.end()){
if (i == items_.end())
{
throw except_item_not_found('\'' + section + '/' + key +
std::string("\' not found."));
std::string("\' not found."));
}
convert(i->second, r);
@ -247,7 +294,7 @@ public:
{
r = get_value_basic<T>(section, key);
}
catch (except_item_not_found& e)
catch (except_item_not_found &e)
{
music::elog << e.what() << std::endl;
throw;
@ -265,11 +312,15 @@ public:
*/
template <class T>
T get_value_safe(std::string const &section, std::string const &key,
T default_value) const {
T default_value) const
{
T r;
try {
try
{
r = get_value_basic<T>(section, key);
} catch (except_item_not_found&) {
}
catch (except_item_not_found &)
{
r = default_value;
music::dlog << "Item \'" << section << "/" << key << " not found in config. Default = \'" << default_value << "\'" << std::endl;
}
@ -284,14 +335,17 @@ public:
* @return the key value (if key found) otherwise default_value
*/
template <class T>
T get_value_safe(std::string const &key, T default_value) const {
T get_value_safe(std::string const &key, T default_value) const
{
return get_value_safe("", key, default_value);
}
//! dumps all key-value pairs to a std::ostream
void dump(std::ostream &out) {
void dump(std::ostream &out)
{
std::map<std::string, std::string>::const_iterator i = items_.begin();
while (i != items_.end()) {
while (i != items_.end())
{
if (i->second.length() > 0)
out << std::setw(24) << std::left << i->first << " = " << i->second
<< std::endl;
@ -299,13 +353,15 @@ public:
}
}
void dump_to_log(void) {
void dump_to_log(void)
{
music::ilog << "List of all configuration options:" << std::endl;
std::map<std::string, std::string>::const_iterator i = items_.begin();
while (i != items_.end()) {
while (i != items_.end())
{
if (i->second.length() > 0)
music::ilog << std::setw(28) << i->first << " = " << i->second
<< std::endl;
<< std::endl;
++i;
}
}
@ -313,20 +369,23 @@ public:
//--- EXCEPTIONS ---
//! runtime error that is thrown if key is not found in getValue
class except_item_not_found : public std::runtime_error {
class except_item_not_found : public std::runtime_error
{
public:
except_item_not_found(std::string itemname)
: std::runtime_error(itemname.c_str()) {}
};
//! runtime error that is thrown if type conversion fails
class except_invalid_conversion : public std::runtime_error {
class except_invalid_conversion : public std::runtime_error
{
public:
except_invalid_conversion(std::string errmsg) : std::runtime_error(errmsg) {}
};
//! runtime error that is thrown if identifier is not found in keys
class ErrIllegalIdentifier : public std::runtime_error {
class ErrIllegalIdentifier : public std::runtime_error
{
public:
ErrIllegalIdentifier(std::string errmsg) : std::runtime_error(errmsg) {}
};
@ -342,7 +401,8 @@ public:
//... converts the string to type bool, returns type bool ...
template <>
inline bool config_file::get_value<bool>(std::string const &strSection,
std::string const &strEntry) const {
std::string const &strEntry) const
{
std::string r1 = get_value<std::string>(strSection, strEntry);
if (r1 == "true" || r1 == "yes" || r1 == "on" || r1 == "1")
return true;
@ -357,17 +417,21 @@ inline bool config_file::get_value<bool>(std::string const &strSection,
template <>
inline bool config_file::get_value_safe<bool>(std::string const &strSection,
std::string const &strEntry,
bool defaultValue) const {
std::string const &strEntry,
bool defaultValue) const
{
std::string r1;
try {
try
{
r1 = get_value_basic<std::string>(strSection, strEntry);
std::transform(r1.begin(), r1.end(),r1.begin(), ::toupper);
std::transform(r1.begin(), r1.end(), r1.begin(), ::toupper);
if (r1 == "YAY" || r1 == "TRUE" || r1 == "YES" || r1 == "ON" || r1 == "1")
return true;
if (r1 == "NAY" || r1 == "FALSE" || r1 == "NO" || r1 == "OFF" || r1 == "0")
return false;
} catch (except_item_not_found&) {
}
catch (except_item_not_found &)
{
return defaultValue;
}
return defaultValue;
@ -376,6 +440,7 @@ inline bool config_file::get_value_safe<bool>(std::string const &strSection,
template <>
inline void
config_file::convert<std::string, std::string>(const std::string &ival,
std::string &oval) const {
std::string &oval) const
{
oval = ival;
}

27
include/cosmology_calculator.hh
View file

@ -255,7 +255,7 @@ public:
music::ilog << "Linear growth factors: D+_target = " << Dplus_target_ << ", D+_start = " << Dplus_start_ << std::endl;
// set up transfer functions and compute normalisation
transfer_function_ = std::move(select_TransferFunction_plugin(cf, cosmo_param_));
transfer_function_ = select_TransferFunction_plugin(cf, cosmo_param_);
transfer_function_->intialise();
if( !transfer_function_->tf_isnormalised_ ){
cosmo_param_.set("pnorm", this->compute_pnorm_from_sigma8() );
@ -268,8 +268,13 @@ public:
music::ilog << std::setw(32) << std::left << "TF supports distinct CDM+baryons"
<< " : " << (transfer_function_->tf_is_distinct() ? "yes" : "no") << std::endl;
music::ilog << std::setw(32) << std::left << "TF minimum wave number"
<< " : " << transfer_function_->get_kmin() << " h/Mpc" << std::endl;
music::ilog << std::setw(32) << std::left << "TF maximum wave number"
<< " : " << transfer_function_->get_kmax() << " h/Mpc" << std::endl;
if( std::sqrt(3.0)* 2.0*M_PI / cf.get_value<double>("setup","BoxLength") * cf.get_value<double>("setup","GridRes")/2 >= transfer_function_->get_kmax() ){
music::elog << "Simulation nyquist mode kny = " << std::sqrt(3.0)* 2.0*M_PI / cf.get_value<double>("setup","BoxLength") * cf.get_value<double>("setup","GridRes")/2 << " h/Mpc is beyond valid range of transfer function!" << std::endl;
}
m_n_s_ = cosmo_param_["n_s"];
m_sqrtpnorm_ = cosmo_param_["sqrtpnorm"];
@ -305,19 +310,19 @@ public:
for (double k = kmin; k < transfer_function_->get_kmax(); k *= 1.01)
{
ofs << std::setw(20) << std::setprecision(10) << k
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm)*Dplus_start_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon)*Dplus_start_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter0)* Dplus_start_ / Dplus_target_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm0)* Dplus_start_ / Dplus_target_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon0)* Dplus_start_ / Dplus_target_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm0)* Dplus_start_ / Dplus_target_, 2.0)
// << std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon)*Dplus_start_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_matter0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_cdm0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, delta_baryon0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_cdm0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::setw(20) << std::setprecision(10) << std::pow(this->get_amplitude(k, theta_baryon0)* Dplus_start_ / Dplus_target_, 2.0)
<< std::endl;
}
}

18
include/cosmology_parameters.hh
View file

@ -107,14 +107,19 @@ namespace cosmology
pmap_["sigma_8"] = cf.get_value_safe<double>("cosmology", "sigma_8", -1.0);
// baryon and non-relativistic matter content
pmap_["Omega_b"] = cf.get_value_safe<double>("cosmology", "Omega_b", defaultp["Omega_b"]);
pmap_["Omega_m"] = cf.get_value_safe<double>("cosmology", "Omega_m", defaultp["Omega_m"]);
pmap_["Omega_b"] = cf.get_value_safe<double>("cosmology", "Omega_b", defaultp["Omega_b"]);
pmap_["Omega_m"] = cf.get_value_safe<double>("cosmology", "Omega_m", defaultp["Omega_m"]);
pmap_["Omega_ncdm"] = cf.get_value_safe<double>("cosmology", "Omega_ncdm", defaultp["Omega_ncdm"]); // TOMA
// pmap_["m_ncdm"] = cf.get_value_safe<double>("cosmology", "m_ncdm", defaultp["m_ncdm"]); //TOMA
// pmap_["N_ncdm"] = this->get("Omega_ncdm")>1e-9? 1 : 0;
// pmap_["T_ncdm"] = cf.get_value_safe<double>("cosmology", "T_ncdm", defaultp["T_ncdm"]); //TOMA
// const int N_ncdm = this->get("N_ncdm");
// massive neutrino species
pmap_["m_nu1"] = cf.get_value_safe<double>("cosmology", "m_nu1", defaultp["m_nu1"]);
pmap_["m_nu2"] = cf.get_value_safe<double>("cosmology", "m_nu2", defaultp["m_nu2"]);
pmap_["m_nu3"] = cf.get_value_safe<double>("cosmology", "m_nu3", defaultp["m_nu3"]);
int N_nu_massive = int(this->get("m_nu1") > 1e-9) + int(this->get("m_nu2") > 1e-9) + int(this->get("m_nu3") > 1e-9);;
const int N_nu_massive = int(this->get("m_nu1") > 1e-9) + int(this->get("m_nu2") > 1e-9) + int(this->get("m_nu3") > 1e-9);
// number ultrarelativistic neutrinos
pmap_["N_ur"] = cf.get_value_safe<double>("cosmology", "N_ur", 3.046 - N_nu_massive);
@ -168,8 +173,8 @@ namespace cosmology
// total relativistic
pmap_["Omega_r"] = this->get("Omega_gamma") + this->get("Omega_nu_massless");
// compute amount of cold dark matter as the rest
pmap_["Omega_c"] = this->get("Omega_m") - this->get("Omega_b") - this->get("Omega_nu_massive");
// compute amount of non-relativistic matter
pmap_["Omega_c"] = std::max(this->get("Omega_m") - this->get("Omega_b") - this->get("Omega_nu_massive") - this->get("Omega_ncdm"),0.0);
if (cf.get_value_safe<bool>("cosmology", "ZeroRadiation", false))
{
@ -177,7 +182,7 @@ namespace cosmology
}
pmap_["f_b"] = this->get("Omega_b") / this->get("Omega_m");
pmap_["f_c"] = 1.0 - this->get("f_b"); // this means we add massive neutrinos to CDM here
pmap_["f_c"] = 1.0 - this->get("f_b"); // this means we add massive neutrinos and ncdm to CDM here
#if 1
// assume zero curvature, take difference from dark energy
@ -202,6 +207,7 @@ namespace cosmology
music::ilog << "sigma_8 = " << std::setw(16) << this->get("sigma_8");
music::ilog << "n_s = " << std::setw(16) << this->get("n_s") << std::endl;
music::ilog << " Omega_c = " << std::setw(16) << this->get("Omega_c") << "Omega_b = " << std::setw(16) << this->get("Omega_b") << "Omega_m = " << std::setw(16) << this->get("Omega_m") << std::endl;
// music::ilog << " Omega_ncdm = " << std::setw(16) << this->get("Omega_ncdm") << "m_ncdm = " << std::setw(16) << this->get("m_ncdm") << this->get("T_ncdm") << "T_ncdm = " << std::setw(16) << this->get("T_ncdm") << std::endl;
music::ilog << " Omega_r = " << std::setw(16) << this->get("Omega_r") << "Omega_nu = " << std::setw(16) << this->get("Omega_nu_massive") << "∑m_nu = " << sum_m_nu << "eV" << std::endl;
music::ilog << " Omega_DE = " << std::setw(16) << this->get("Omega_DE") << "w_0 = " << std::setw(16) << this->get("w_0") << "w_a = " << std::setw(16) << this->get("w_a") << std::endl;
//music::ilog << " Omega_k = " << 1.0 - this->get("Omega_m") - this->get("Omega_r") - this->get("Omega_DE") << std::endl;

3
include/general.hh
View file

@ -20,6 +20,7 @@
#include <complex>
#include <map>
#include <memory>
#if defined(USE_MPI)
#include <mpi.h>
@ -184,6 +185,8 @@ inline void multitask_sync_barrier(void)
#endif
}
extern size_t global_mem_high_mark, local_mem_high_mark;
namespace CONFIG
{
extern int MPI_thread_support;

209
include/grid_ghosts.hh Normal file
View file

@ -0,0 +1,209 @@
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2022 by Oliver Hahn
//
// monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#pragma once
#include <array>
#include <vector>
#include <numeric>
#include <general.hh>
#include <math/vec3.hh>
template <int numghosts, bool haveleft, bool haveright, typename grid_t>
struct grid_with_ghosts
{
using data_t = typename grid_t::data_t;
using vec3 = std::array<real_t, 3>;
static constexpr bool is_distributed_trait = grid_t::is_distributed_trait;
static constexpr int num_ghosts = numghosts;
static constexpr bool have_left = haveleft, have_right = haveright;
std::vector<data_t> boundary_left_, boundary_right_;
std::vector<int> local0starts_;
const grid_t &gridref;
size_t nx_, ny_, nz_, nzp_;
//... determine communication offsets
std::vector<ptrdiff_t> offsets_, sizes_;
int get_task(ptrdiff_t index) const
{
int itask = 0;
while (itask < MPI::get_size() - 1 && offsets_[itask + 1] <= index)
++itask;
return itask;
}
explicit grid_with_ghosts(const grid_t &g)
: gridref(g), nx_(g.n_[0]), ny_(g.n_[1]), nz_(g.n_[2]), nzp_(g.n_[2]+2)
{
if (is_distributed_trait)
{
int ntasks(MPI::get_size());
offsets_.assign(ntasks+1, 0);
sizes_.assign(ntasks, 0);
MPI_Allgather(&g.local_0_size_, 1, MPI_LONG_LONG, &sizes_[0], 1,
MPI_LONG_LONG, MPI_COMM_WORLD);
MPI_Allgather(&g.local_0_start_, 1, MPI_LONG_LONG, &offsets_[0], 1,
MPI_LONG_LONG, MPI_COMM_WORLD);
for( int i=0; i< CONFIG::MPI_task_size; i++ ){
if( offsets_[i+1] < offsets_[i] + sizes_[i] ) offsets_[i+1] = offsets_[i] + sizes_[i];
}
update_ghosts_allow_multiple( g );
}
}
void update_ghosts_allow_multiple( const grid_t &g )
{
#if defined(USE_MPI)
//... exchange boundary
if( have_left ) boundary_left_.assign(num_ghosts * ny_ * nzp_, data_t{0.0});
if( have_right ) boundary_right_.assign(num_ghosts * ny_ * nzp_, data_t{0.0});
size_t slicesz = ny_ * nzp_;
MPI_Status status;
std::vector<MPI_Request> req;
MPI_Request temp_req;
if( have_right ){
for( int itask=0; itask<CONFIG::MPI_task_size; ++itask ){
for( size_t i=0; i<num_ghosts; ++i ){
ptrdiff_t iglobal_request = (offsets_[itask] + sizes_[itask] + i) % g.n_[0];
if( iglobal_request >= g.local_0_start_ && iglobal_request < g.local_0_start_ + g.local_0_size_ ){
size_t ii = iglobal_request - g.local_0_start_;
MPI_Isend( &g.relem(ii*slicesz), slicesz, MPI::get_datatype<data_t>(), itask, iglobal_request, MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
}
}
}
//--- receive data ------------------------------------------------------------
#pragma omp parallel if(CONFIG::MPI_threads_ok)
{
MPI_Status status;
#pragma omp for
for( size_t i=0; i<num_ghosts; ++i ){
ptrdiff_t iglobal_request = (g.local_0_start_ + g.local_0_size_ + i) % g.n_[0];
int recvfrom = get_task(iglobal_request);
//#pragma omp critical // need critical region here if we do "MPI_THREAD_FUNNELED",
{
// receive data slice and check for MPI errors when in debug mode
status.MPI_ERROR = MPI_SUCCESS;
MPI_Recv(&boundary_right_[i*slicesz], (int)slicesz, MPI::get_datatype<data_t>(), recvfrom, (int)iglobal_request, MPI_COMM_WORLD, &status);
assert(status.MPI_ERROR == MPI_SUCCESS);
}
}
}
}
MPI_Barrier( MPI_COMM_WORLD );
if( have_left ){
for( int itask=0; itask<CONFIG::MPI_task_size; ++itask ){
for( size_t i=0; i<num_ghosts; ++i ){
ptrdiff_t iglobal_request = (offsets_[itask] + g.n_[0] - num_ghosts + i) % g.n_[0];
if( iglobal_request >= g.local_0_start_ && iglobal_request < g.local_0_start_ + g.local_0_size_ ){
size_t ii = iglobal_request - g.local_0_start_;
MPI_Isend( &g.relem(ii*slicesz), slicesz, MPI::get_datatype<data_t>(), itask, iglobal_request, MPI_COMM_WORLD, &temp_req);
req.push_back(temp_req);
}
}
}
//--- receive data ------------------------------------------------------------
#pragma omp parallel if(CONFIG::MPI_threads_ok)
{
MPI_Status status;
#pragma omp for
for( size_t i=0; i<num_ghosts; ++i ){
ptrdiff_t iglobal_request = (g.local_0_start_ + g.n_[0] - num_ghosts + i) % g.n_[0];
int recvfrom = get_task(iglobal_request);
//#pragma omp critical // need critical region here if we do "MPI_THREAD_FUNNELED",
{
// receive data slice and check for MPI errors when in debug mode
status.MPI_ERROR = MPI_SUCCESS;
MPI_Recv(&boundary_left_[i*slicesz], (int)slicesz, MPI::get_datatype<data_t>(), recvfrom, (int)iglobal_request, MPI_COMM_WORLD, &status);
assert(status.MPI_ERROR == MPI_SUCCESS);
}
}
}
}
MPI_Barrier( MPI_COMM_WORLD );
for (size_t i = 0; i < req.size(); ++i)
{
// need to set status as wait does not necessarily modify it
// c.f. http://www.open-mpi.org/community/lists/devel/2007/04/1402.php
status.MPI_ERROR = MPI_SUCCESS;
// std::cout << "task " << CONFIG::MPI_task_rank << " : checking request No" << i << std::endl;
int flag(1);
MPI_Test(&req[i], &flag, &status);
if( !flag ){
std::cout << "task " << CONFIG::MPI_task_rank << " : request No" << i << " unsuccessful" << std::endl;
}
MPI_Wait(&req[i], &status);
// std::cout << "---> ok!" << std::endl;
assert(status.MPI_ERROR == MPI_SUCCESS);
}
MPI_Barrier(MPI_COMM_WORLD);
#endif
}
data_t relem(const ptrdiff_t& i, const ptrdiff_t& j, const ptrdiff_t&k ) const noexcept
{
return this->relem({i,j,k});
}
data_t relem(const std::array<ptrdiff_t, 3> &pos) const noexcept
{
const ptrdiff_t ix = pos[0];
const ptrdiff_t iy = (pos[1]+gridref.n_[1])%gridref.n_[1];
const ptrdiff_t iz = (pos[2]+gridref.n_[2])%gridref.n_[2];
if( is_distributed_trait ){
const ptrdiff_t localix = ix;
if( localix < 0 ){
return boundary_left_[((localix+num_ghosts)*ny_+iy)*nzp_+iz];
}else if( localix >= gridref.local_0_size_ ){
return boundary_right_[((localix-gridref.local_0_size_)*ny_+iy)*nzp_+iz];
}else{
return gridref.relem(localix, iy, iz);
}
}
return gridref.relem((ix+gridref.n_[0])%gridref.n_[0], iy, iz);
}
};

5
include/grid_interpolate.hh
View file

@ -175,13 +175,14 @@ struct grid_interpolate
MPI_Alltoall(&sendcounts[0], 1, MPI_INT, &recvcounts[0], 1, MPI_INT, MPI_COMM_WORLD);
size_t tot_receive = recvcounts[0], tot_send = sendcounts[0];
size_t tot_receive = recvcounts[0];
// size_t tot_send = sendcounts[0];
for (int i = 1; i < MPI::get_size(); ++i)
{
sendoffsets[i] = sendcounts[i - 1] + sendoffsets[i - 1];
recvoffsets[i] = recvcounts[i - 1] + recvoffsets[i - 1];
tot_receive += recvcounts[i];
tot_send += sendcounts[i];
// tot_send += sendcounts[i];
}
std::vector<vec3> recvbuf(tot_receive/3,{0.,0.,0.});

17
include/math/vec3.hh
View file

@ -30,7 +30,7 @@ public:
//! empty constructor
vec3_t()
: x(data_[0]),y(data_[1]),z(data_[2]){}
: data_{{T(0),T(0),T(0)}},x(data_[0]),y(data_[1]),z(data_[2]){}
//! copy constructor
vec3_t( const vec3_t<T> &v)
@ -45,12 +45,15 @@ public:
: data_(std::move(v.data_)), x(data_[0]), y(data_[1]), z(data_[2]){}
//! construct vec3_t from initializer list
template<typename ...E>
vec3_t(E&&...e)
: data_{{std::forward<E>(e)...}}, x{data_[0]}, y{data_[1]}, z{data_[2]}
{}
// vec3_t( T a, T b, T c )
// : data_{{a,b,c}}, x(data_[0]), y(data_[1]), z(data_[2]){}
// template<typename ...E>
// vec3_t(E&&...e)
// : data_{{std::forward<E>(e)...}}, x{data_[0]}, y{data_[1]}, z{data_[2]}
// {}
vec3_t( T a, T b, T c )
: data_{{a,b,c}}, x(data_[0]), y(data_[1]), z(data_[2]){}
explicit vec3_t( T a)
: data_{{a,a,a}}, x(data_[0]), y(data_[1]), z(data_[2]){}
//! bracket index access to vector components
T &operator[](size_t i) noexcept{ return data_[i];}

125
include/memory_stat.hh Normal file
View file

@ -0,0 +1,125 @@
/*
* Author: David Robert Nadeau
* Site: http://NadeauSoftware.com/
* License: Creative Commons Attribution 3.0 Unported License
* http://creativecommons.org/licenses/by/3.0/deed.en_US
*/
#pragma once
namespace memory
{
#if defined(_WIN32)
#include <windows.h>
#include <psapi.h>
#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
#include <unistd.h>
#include <sys/resource.h>
#if defined(__APPLE__) && defined(__MACH__)
#include <mach/mach.h>
#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
#include <fcntl.h>
#include <procfs.h>
#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
#include <stdio.h>
#endif
#else
#error "Cannot define getPeakRSS( ) or getCurrentRSS( ) for an unknown OS."
#endif
/**
* Returns the peak (maximum so far) resident set size (physical
* memory use) measured in bytes, or zero if the value cannot be
* determined on this OS.
*/
inline size_t getPeakRSS( )
{
#if defined(_WIN32)
/* Windows -------------------------------------------------- */
PROCESS_MEMORY_COUNTERS info;
GetProcessMemoryInfo( GetCurrentProcess( ), &info, sizeof(info) );
return (size_t)info.PeakWorkingSetSize;
#elif (defined(_AIX) || defined(__TOS__AIX__)) || (defined(__sun__) || defined(__sun) || defined(sun) && (defined(__SVR4) || defined(__svr4__)))
/* AIX and Solaris ------------------------------------------ */
struct psinfo psinfo;
int fd = -1;
if ( (fd = open( "/proc/self/psinfo", O_RDONLY )) == -1 )
return (size_t)0L; /* Can't open? */
if ( read( fd, &psinfo, sizeof(psinfo) ) != sizeof(psinfo) )
{
close( fd );
return (size_t)0L; /* Can't read? */
}
close( fd );
return (size_t)(psinfo.pr_rssize * 1024L);
#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
/* BSD, Linux, and OSX -------------------------------------- */
struct rusage rusage;
getrusage( RUSAGE_SELF, &rusage );
#if defined(__APPLE__) && defined(__MACH__)
return (size_t)rusage.ru_maxrss;
#else
return (size_t)(rusage.ru_maxrss * 1024L);
#endif
#else
/* Unknown OS ----------------------------------------------- */
return (size_t)0L; /* Unsupported. */
#endif
}
/**
* Returns the current resident set size (physical memory use) measured
* in bytes, or zero if the value cannot be determined on this OS.
*/
inline size_t getCurrentRSS( )
{
#if defined(_WIN32)
/* Windows -------------------------------------------------- */
PROCESS_MEMORY_COUNTERS info;
GetProcessMemoryInfo( GetCurrentProcess( ), &info, sizeof(info) );
return (size_t)info.WorkingSetSize;
#elif defined(__APPLE__) && defined(__MACH__)
/* OSX ------------------------------------------------------ */
struct mach_task_basic_info info;
mach_msg_type_number_t infoCount = MACH_TASK_BASIC_INFO_COUNT;
if ( task_info( mach_task_self( ), MACH_TASK_BASIC_INFO,
(task_info_t)&info, &infoCount ) != KERN_SUCCESS )
return (size_t)0L; /* Can't access? */
return (size_t)info.resident_size;
#elif defined(__linux__) || defined(__linux) || defined(linux) || defined(__gnu_linux__)
/* Linux ---------------------------------------------------- */
long rss = 0L;
FILE* fp = NULL;
if ( (fp = fopen( "/proc/self/statm", "r" )) == NULL )
return (size_t)0L; /* Can't open? */
if ( fscanf( fp, "%*s%ld", &rss ) != 1 )
{
fclose( fp );
return (size_t)0L; /* Can't read? */
}
fclose( fp );
return (size_t)rss * (size_t)sysconf( _SC_PAGESIZE);
#else
/* AIX, BSD, Solaris, and Unknown OS ------------------------ */
return (size_t)0L; /* Unsupported. */
#endif
}
};

334
include/particle_generator.hh
View file

@ -17,6 +17,7 @@
#pragma once
#include <math/vec3.hh>
#include <grid_ghosts.hh>
#include <grid_interpolate.hh>
#if defined(USE_HDF5)
@ -29,13 +30,20 @@ namespace particle
enum lattice
{
lattice_glass = -1,
lattice_masked = -2, // masked lattices are SC lattices with a specifiable mask leaving out particles
lattice_glass = -1, // glass: needs a glass file
lattice_sc = 0, // SC : simple cubic
lattice_bcc = 1, // BCC: body-centered cubic
lattice_fcc = 2, // FCC: face-centered cubic
lattice_rsc = 3, // RSC: refined simple cubic
};
// const std::vector<std::vector<bool>> lattice_masks =
// {
// // mask from Richings et al. https://arxiv.org/pdf/2005.14495.pdf
// {true,true,true,true,true,true,true,false},
// };
const std::vector<std::vector<vec3_t<real_t>>> lattice_shifts =
{
// first shift must always be zero! (otherwise set_positions and set_velocities break)
@ -149,20 +157,75 @@ namespace particle
private:
particle::container particles_;
size_t global_num_particles_;
static constexpr int masksize_ = 2;
std::array<int, masksize_*masksize_*masksize_> particle_type_mask_;
inline int get_mask_value( const vec3_t<size_t>& global_idx_3d ) const
{
int sig = ((global_idx_3d[0]%masksize_)*masksize_
+global_idx_3d[1]%masksize_)*masksize_
+global_idx_3d[2]%masksize_;
return particle_type_mask_[sig];
}
template< typename ggrid_t >
inline real_t get_mean_mask_value( const ggrid_t& gg_field, const vec3_t<size_t>& global_idx_3d, size_t i, size_t j, size_t k, int lattice_index ) const
{
ptrdiff_t ox = (ptrdiff_t)i-(ptrdiff_t)(global_idx_3d[0]%masksize_);
ptrdiff_t oy = (ptrdiff_t)j-(ptrdiff_t)(global_idx_3d[1]%masksize_);
ptrdiff_t oz = (ptrdiff_t)k-(ptrdiff_t)(global_idx_3d[2]%masksize_);
size_t count_full{0}, count_masked{0};
real_t mean_full{0.0}, mean_masked{0.0};
for( size_t i=0; i<masksize_; ++i ){
for( size_t j=0; j<masksize_; ++j ){
for( size_t k=0; k<masksize_; ++k ){
mean_full += gg_field.relem( ox+i, oy+j, oz+k );
++count_full;
if( particle_type_mask_[4*i+2*j+k] == lattice_index ){
mean_masked += gg_field.relem( ox+i, oy+j, oz+k );
++count_masked;
}
}
}
}
mean_full /= count_full;
mean_masked /= count_masked;
return mean_masked - mean_full;
}
public:
lattice_generator(lattice lattice_type, const bool b64reals, const bool b64ids, const bool bwithmasses, size_t IDoffset, const field_t &field, config_file &cf)
{
if (lattice_type != lattice_glass)
{
music::wlog << "Glass ICs will currently be incorrect due to disabled ghost zone updates! ";
/**
* @brief Construct a new lattice generator object
*
* @param lattice_type Type of the lattice (currently: 0=SC, 1=BCC, 2=FCC, 3=double SC, -1=glass, -2=masked SC)
* @param lattice_index Index of the 'atom type', i.e. whether CDM, baryons, ... (currently only supports 0 or 1)
* @param b64reals Boolean whether 64bit floating point shall be used to store positions/velocities/masses
* @param b64ids Boolean whether 64bit integers should be used for IDS
* @param bwithmasses Boolean whether distinct masses for each particle shall be stored
* @param IDoffset The global ID offset to be applied to this generation of particles
* @param field Reference to the field from which particles shall be generated (used only to set dimensions)
* @param cf Reference to the config_file object
*/
lattice_generator(lattice lattice_type, int lattice_index, const bool b64reals, const bool b64ids, const bool bwithmasses, size_t IDoffset, const field_t &field, config_file &cf)
: global_num_particles_(0)
{
// initialise the particle mask with zeros (only used if lattice_type==lattice_masked)
for( auto& m : particle_type_mask_) m = 0;
if (lattice_type >= 0) // These are the Bravais lattices
{
// number of modes present in the field
const size_t num_p_in_load = field.local_size();
// unless SC lattice is used, particle number is a multiple of the number of modes (=num_p_in_load):
const size_t overload = 1ull << std::max<int>(0, lattice_type); // 1 for sc, 2 for bcc, 4 for fcc, 8 for rsc
// allocate memory for all local particles
particles_.allocate(overload * num_p_in_load, b64reals, b64ids, bwithmasses);
// set the global number of particles for this lattice_type and lattice_index
global_num_particles_ = field.global_size() * overload;
// set particle IDs to the Lagrangian coordinate (1D encoded) with additionally the field shift encoded as well
IDoffset = IDoffset * overload * field.global_size();
@ -188,8 +251,99 @@ namespace particle
}
}
}
else
else if( lattice_type == lattice_masked )
{
if( field.global_size()%8 != 0 ){
music::elog << "For masked lattice type, linear field resolution must be a multiple of two." << std::endl;
abort();
}
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For masked lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
// set the particle mask
if( cf.get_value_safe("setup","DoBaryons",false) )
{
unsigned mask_type = cf.get_value_safe("setup","ParticleMaskType",3);
// mask everywehere 0, except the last element
for( auto& m : particle_type_mask_) m = -1;
particle_type_mask_[0] = 0; // CDM at corner of unit cube
if( mask_type & 1<<0 ) {
// edge centers
particle_type_mask_[1] = 0; // CDM
particle_type_mask_[2] = 0; // CDM
particle_type_mask_[4] = 0; // CDM
}
if( mask_type & 1<<1 ){
// face centers
particle_type_mask_[3] = 0; // CDM
particle_type_mask_[5] = 0; // CDM
particle_type_mask_[6] = 0; // CDM
}
particle_type_mask_[7] = 1; // baryon at cell center (aka opposite corner)
}else{
// mask everywhere 0, all particle locations are occupied by CDM
for( auto& m : particle_type_mask_) m = 0;
}
// count number of particles taking into account masking
size_t ipcount = 0;
for (size_t i = 0; i < field.rsize(0); ++i)
{
for (size_t j = 0; j < field.rsize(1); ++j)
{
for (size_t k = 0; k < field.rsize(2); ++k)
{
if( this->get_mask_value(field.get_cell_idx_3d(i,j,k)) != lattice_index ) continue;
++ipcount;
}
}
}
// set global number of particles
#if defined(USE_MPI)
MPI_Allreduce( &ipcount, &global_num_particles_, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD );
#else
global_num_particles_ = ipcount;
#endif
// number of modes present in the field
const size_t num_p_in_load = ipcount;
// allocate memory for all local particles
particles_.allocate(num_p_in_load, b64reals, b64ids, bwithmasses);
// set particle IDs to the Lagrangian coordinate (1D encoded) with additionally the field shift encoded as well
IDoffset = IDoffset * field.global_size();
for (size_t i = 0, ipcount = 0; i < field.rsize(0); ++i)
{
for (size_t j = 0; j < field.rsize(1); ++j)
{
for (size_t k = 0; k < field.rsize(2); ++k)
{
if( this->get_mask_value(field.get_cell_idx_3d(i,j,k)) != lattice_index ) continue;
if (b64ids)
{
particles_.set_id64(ipcount, IDoffset + field.get_cell_idx_1d(i, j, k));
}
else
{
particles_.set_id32(ipcount, IDoffset + field.get_cell_idx_1d(i, j, k));
}
++ipcount;
}
}
}
}
else if( lattice_type == lattice_glass )
{
music::wlog << "Glass ICs will currently be incorrect due to disabled ghost zone updates! ";
glass_ptr_ = std::make_unique<glass>( cf, field );
particles_.allocate(glass_ptr_->size(), b64reals, b64ids, false);
@ -206,14 +360,24 @@ namespace particle
}
}
}
music::ilog << "Created Particles [" << lattice_index << "] : " << global_num_particles_ << std::endl;
}
// invalidates field, phase shifted to unspecified position after return
void set_masses(const lattice lattice_type, bool is_second_lattice, const real_t munit, const bool b64reals, field_t &field, config_file &cf)
void set_masses(const lattice lattice_type, int lattice_index, const real_t munit, const bool b64reals, field_t &field, config_file &cf)
{
// works only for Bravais types
if (lattice_type >= 0)
// works only for Bravais types and masked type
assert( lattice_type>=0 || lattice_type==lattice_masked );
if (lattice_type >= 0) // Bravais lattices
{
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For Bravais lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
const size_t overload = 1ull << std::max<int>(0, lattice_type); // 1 for sc, 2 for bcc, 4 for fcc, 8 for rsc
const size_t num_p_in_load = field.local_size();
const real_t pmeanmass = munit / real_t(field.global_size()* overload);
@ -225,7 +389,7 @@ namespace particle
for (int ishift = 0; ishift < (1 << lattice_type); ++ishift)
{
// if we are dealing with the secondary lattice, apply a global shift
if (ishift == 0 && is_second_lattice)
if (ishift == 0 && lattice_index > 0)
{
field.shift_field(second_lattice_shift[lattice_type]);
}
@ -261,7 +425,69 @@ namespace particle
music::ilog << "Particle Mass : mean/munit = " << mean_pm/munit << " ; fractional RMS = " << std_pm / mean_pm * 100.0 << "%" << std::endl;
if(std_pm / mean_pm > 0.1 ) music::wlog << "Particle mass perturbation larger than 10%, consider decreasing \n\t the starting redshift or disabling baryon decaying modes." << std::endl;
if(bmass_negative) music::elog << "Negative particle mass produced! Decrease the starting \n\t redshift or disable baryon decaying modes!" << std::endl;
} else if( lattice_type == lattice_masked ) {
if( field.global_size()%8 != 0 ){
music::elog << "For masked lattice type, linear field resolution must be a multiple of two." << std::endl;
abort();
}
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For masked lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
grid_with_ghosts<1, true, true, field_t> gg_field(field);
const real_t pmeanmass = munit / global_num_particles_;
bool bmass_negative = false;
real_t mean_pm = 0.0;//field.mean() * pmeanmass;
real_t std_pm = 0.0; //field.std() * pmeanmass;
// read out values from phase shifted field and set assoc. particle's value
for (size_t i = 0, ipcount = 0; i < field.size(0); ++i)
{
for (size_t j = 0; j < field.size(1); ++j)
{
for (size_t k = 0; k < field.size(2); ++k)
{
const auto idx3 = field.get_cell_idx_3d(i,j,k);
if( this->get_mask_value(idx3) != lattice_index ) continue;
const auto mean_mask = this->get_mean_mask_value( gg_field, idx3, i, j, k, lattice_index );
// get
const auto pmass = pmeanmass * (field.relem(i, j, k) - mean_mask);
// check for negative mass
bmass_negative |= pmass<0.0;
// set
if (b64reals) particles_.set_mass64(ipcount++, pmass);
else particles_.set_mass32(ipcount++, pmass);
// statistics
mean_pm += pmass;
std_pm += pmass*pmass;
}
}
}
#if defined(USE_MPI)
{
double local_mean_pm = mean_pm, local_std_pm = std_pm;
MPI_Allreduce( &local_mean_pm, &mean_pm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
MPI_Allreduce( &local_std_pm, &std_pm, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD );
}
#endif
mean_pm /= global_num_particles_;
std_pm /= global_num_particles_;
std_pm -= mean_pm*mean_pm;
// diagnostics
music::ilog << "Particle Mass : mean/munit = " << mean_pm/munit << " ; fractional RMS = " << std_pm / mean_pm * 100.0 << "%" << std::endl;
if(std_pm / mean_pm > 0.1 ) music::wlog << "Particle mass perturbation larger than 10%, consider decreasing \n\t the starting redshift or disabling baryon decaying modes." << std::endl;
if(bmass_negative) music::elog << "Negative particle mass produced! Decrease the starting \n\t redshift or disable baryon decaying modes!" << std::endl;
}else{
// should not happen
music::elog << "Cannot have individual particle masses for glasses!" << std::endl;
@ -270,15 +496,21 @@ namespace particle
}
// invalidates field, phase shifted to unspecified position after return
void set_positions(const lattice lattice_type, bool is_second_lattice, int idim, real_t lunit, const bool b64reals, field_t &field, config_file &cf)
void set_positions(const lattice lattice_type, int lattice_index, int idim, real_t lunit, const bool b64reals, field_t &field, config_file &cf)
{
if (lattice_type >= 0)
if (lattice_type >= 0) // Bravais lattice
{
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For Bravais lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
const size_t num_p_in_load = field.local_size();
for (int ishift = 0; ishift < (1 << lattice_type); ++ishift)
{
// if we are dealing with the secondary lattice, apply a global shift
if (ishift == 0 && is_second_lattice)
if (ishift == 0 && lattice_index==1)
{
field.shift_field(second_lattice_shift[lattice_type]);
}
@ -296,7 +528,7 @@ namespace particle
{
for (size_t k = 0; k < field.size(2); ++k)
{
auto pos = field.template get_unit_r_shifted<real_t>(i, j, k, lattice_shifts[lattice_type][ishift] + (is_second_lattice ? second_lattice_shift[lattice_type] : vec3_t<real_t>{real_t(0.), real_t(0.), real_t(0.)}));
auto pos = field.template get_unit_r_shifted<real_t>(i, j, k, lattice_shifts[lattice_type][ishift] + (lattice_index==1 ? second_lattice_shift[lattice_type] : vec3_t<real_t>{real_t(0.), real_t(0.), real_t(0.)}));
if (b64reals)
{
particles_.set_pos64(ipcount++, idim, pos[idim] * lunit + field.relem(i, j, k));
@ -310,6 +542,38 @@ namespace particle
}
}
}
else if( lattice_type == lattice_masked )
{
if( field.global_size()%8 != 0 ){
music::elog << "For masked lattice type, linear field resolution must be a multiple of two." << std::endl;
abort();
}
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For masked lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
for (size_t i = 0, ipcount = 0; i < field.size(0); ++i)
{
for (size_t j = 0; j < field.size(1); ++j)
{
for (size_t k = 0; k < field.size(2); ++k)
{
if( this->get_mask_value(field.get_cell_idx_3d(i,j,k)) != lattice_index ) continue;
// get position (in box units) of the current cell of 3d array 'field'
auto pos = field.template get_unit_r<real_t>(i, j, k);
// add the displacement to get the particle position
if (b64reals){
particles_.set_pos64(ipcount++, idim, pos[idim] * lunit + field.relem(i, j, k));
}else{
particles_.set_pos32(ipcount++, idim, pos[idim] * lunit + field.relem(i, j, k));
}
}
}
}
}
else
{
glass_ptr_->update_ghosts( field );
@ -330,15 +594,20 @@ namespace particle
}
}
void set_velocities(lattice lattice_type, bool is_second_lattice, int idim, const bool b64reals, field_t &field, config_file &cf)
void set_velocities(lattice lattice_type, int lattice_index, int idim, const bool b64reals, field_t &field, config_file &cf)
{
if (lattice_type >= 0)
if (lattice_type >= 0) // Bravais lattice
{
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For Bravais lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
const size_t num_p_in_load = field.local_size();
for (int ishift = 0; ishift < (1 << lattice_type); ++ishift)
{
// if we are dealing with the secondary lattice, apply a global shift
if (ishift == 0 && is_second_lattice)
if (ishift == 0 && lattice_index==1)
{
field.shift_field(second_lattice_shift[lattice_type]);
}
@ -368,6 +637,35 @@ namespace particle
}
}
}
else if( lattice_type == lattice_masked )
{
if( field.global_size()%8 != 0 ){
music::elog << "For masked lattice type, linear field resolution must be a multiple of two." << std::endl;
abort();
}
if( lattice_index > 1 || lattice_index < 0 ){
music::elog << "For masked lattice type, lattice index must be 0 or 1." << std::endl;
abort();
}
for (size_t i = 0, ipcount = 0; i < field.size(0); ++i)
{
for (size_t j = 0; j < field.size(1); ++j)
{
for (size_t k = 0; k < field.size(2); ++k)
{
if( this->get_mask_value(field.get_cell_idx_3d(i,j,k)) != lattice_index ) continue;
if (b64reals){
particles_.set_vel64(ipcount++, idim, field.relem(i, j, k));
}else{
particles_.set_vel32(ipcount++, idim, field.relem(i, j, k));
}
}
}
}
}
else
{
glass_ptr_->update_ghosts( field );

4
include/particle_plt.hh
View file

@ -421,7 +421,7 @@ private:
if( !is_in(i,j,k,mat_invrecip) ){
auto average_lv = [&]( const auto& t1, const auto& t2, const auto& t3, vec3_t<real_t>& v, vec3_t<real_t>& l ) {
v = real_t(0.0); l = real_t(0.0);
v = vec3_t<real_t>(0.0); l = vec3_t<real_t>(0.0);
int count(0);
auto add_lv = [&]( auto it ) -> void {
@ -586,4 +586,4 @@ public:
};
}
}

4
src/cosmology_parameters.cc
View file

@ -37,6 +37,7 @@ parameters::defaultmmap_t parameters::default_pmaps_
{"Omega_m", 0.3158},
{"Omega_b", 0.04938898},
{"Omega_DE", 0.6842},
{"Omega_ncdm", 0.0},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96605},
@ -55,6 +56,7 @@ parameters::defaultmmap_t parameters::default_pmaps_
{"Omega_m", 0.3106},
{"Omega_b", 0.04897284},
{"Omega_DE", 0.6894},
{"Omega_ncdm", 0.0},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96824},
@ -73,6 +75,7 @@ parameters::defaultmmap_t parameters::default_pmaps_
{"Omega_m", 0.3134},
{"Omega_b", 0.04919537},
{"Omega_DE", 0.6866},
{"Omega_ncdm", 0.0},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96654},
@ -91,6 +94,7 @@ parameters::defaultmmap_t parameters::default_pmaps_
{"Omega_m", 0.3099},
{"Omega_b", 0.048891054},
{"Omega_DE", 0.6901},
{"Omega_ncdm", 0.0},
{"w_0", -1.0},
{"w_a", 0.0},
{"n_s", 0.96822},

19
src/grid_fft.cc
View file

@ -19,6 +19,24 @@
#include <grid_fft.hh>
#include <thread>
#include "memory_stat.hh"
void memory_report(void)
{
//... report memory usage
size_t curr_mem_high_mark = 0;
local_mem_high_mark = memory::getCurrentRSS();
#if defined(USE_MPI)
MPI_Allreduce(&local_mem_high_mark, &curr_mem_high_mark, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, MPI_COMM_WORLD);
#else
curr_mem_high_mark = local_mem_high_mark;
#endif
if( curr_mem_high_mark > 1.1*global_mem_high_mark ){
music::ilog << std::setw(57) << std::setfill(' ') << std::right << "mem high: " << std::setw(8) << curr_mem_high_mark/(1ull<<20) << " MBytes / task" << std::endl;
global_mem_high_mark = curr_mem_high_mark;
}
}
template <typename data_t, bool bdistributed>
void Grid_FFT<data_t, bdistributed>::allocate(void)
{
@ -175,6 +193,7 @@ void Grid_FFT<data_t, bdistributed>::allocate(void)
#endif //// of #ifdef #else USE_MPI ////////////////////////////////////////////////////////////////////////////////////
}
ballocated_ = true;
memory_report();
}
template <typename data_t, bool bdistributed>

47
src/ic_generator.cc
View file

@ -62,9 +62,9 @@ std::unique_ptr<cosmology::calculator> the_cosmo_calc;
*/
int initialise( config_file& the_config )
{
the_random_number_generator = std::move(select_RNG_plugin(the_config));
the_random_number_generator = select_RNG_plugin(the_config);
the_cosmo_calc = std::make_unique<cosmology::calculator>(the_config);
the_output_plugin = std::move(select_output_plugin(the_config, the_cosmo_calc));
the_output_plugin = select_output_plugin(the_config, the_cosmo_calc);
return 0;
}
@ -116,7 +116,8 @@ int run( config_file& the_config )
: ((lattice_str=="fcc")? particle::lattice_fcc
: ((lattice_str=="rsc")? particle::lattice_rsc
: ((lattice_str=="glass")? particle::lattice_glass
: particle::lattice_sc))));
: ((lattice_str=="masked")? particle::lattice_masked
: particle::lattice_sc)))));
//--------------------------------------------------------------------------------------------------------
//! apply fixing of the complex mode amplitude following Angulo & Pontzen (2016) [https://arxiv.org/abs/1603.05253]
@ -149,12 +150,12 @@ int run( config_file& the_config )
//--------------------------------------------------------------------------------------------------------
//! add beyond box tidal field modes following Schmidt et al. (2018) [https://arxiv.org/abs/1803.03274]
bool bAddExternalTides = the_config.contains_key("cosmology", "LSS_aniso_lx")
& the_config.contains_key("cosmology", "LSS_aniso_ly")
& the_config.contains_key("cosmology", "LSS_aniso_lz");
&& the_config.contains_key("cosmology", "LSS_aniso_ly")
&& the_config.contains_key("cosmology", "LSS_aniso_lz");
if( bAddExternalTides && !( the_config.contains_key("cosmology", "LSS_aniso_lx")
| the_config.contains_key("cosmology", "LSS_aniso_ly")
| the_config.contains_key("cosmology", "LSS_aniso_lz") ))
|| the_config.contains_key("cosmology", "LSS_aniso_ly")
|| the_config.contains_key("cosmology", "LSS_aniso_lz") ))
{
music::elog << "Not all dimensions of LSS_aniso_l{x,y,z} specified! Will ignore external tidal field!" << std::endl;
bAddExternalTides = false;
@ -367,10 +368,10 @@ int run( config_file& the_config )
// phi = - delta / k^2
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
music::ilog << "Generating LPT fields...." << std::endl;
music::ilog << "\n>>> Generating LPT fields.... <<<\n" << std::endl;
double wtime = get_wtime();
music::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(1) term" << std::flush;
music::ilog << std::setw(79) << std::setfill('.') << std::left << ">> Computing phi(1) term" << std::endl;
phi.FourierTransformForward(false);
phi.assign_function_of_grids_kdep([&](auto k, auto wn) {
@ -381,7 +382,7 @@ int run( config_file& the_config )
phi.zero_DC_mode();
music::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime() - wtime << "s" << std::endl;
music::ilog << std::setw(70) << std::setfill(' ') << std::right << "took : " << std::setw(8) << get_wtime() - wtime << "s" << std::endl;
//======================================================================
//... compute 2LPT displacement potential ....
@ -392,7 +393,7 @@ int run( config_file& the_config )
phi2.FourierTransformForward(false);
wtime = get_wtime();
music::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(2) term" << std::flush;
music::ilog << std::setw(79) << std::setfill('.') << std::left << ">> Computing phi(2) term" << std::endl;
Conv.convolve_SumOfHessians(phi, {0, 0}, phi, {1, 1}, {2, 2}, op::assign_to(phi2));
Conv.convolve_Hessians(phi, {1, 1}, phi, {2, 2}, op::add_to(phi2));
Conv.convolve_Hessians(phi, {0, 1}, phi, {0, 1}, op::subtract_from(phi2));
@ -411,7 +412,7 @@ int run( config_file& the_config )
}
phi2.apply_InverseLaplacian();
music::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime() - wtime << "s" << std::endl;
music::ilog << std::setw(70) << std::setfill(' ') << std::right << "took : " << std::setw(8) << get_wtime() - wtime << "s" << std::endl;
if (bAddExternalTides)
{
@ -432,19 +433,18 @@ int run( config_file& the_config )
//... phi3 = phi3a - 10/7 phi3b
//... 3a term ...
wtime = get_wtime();
music::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3a) term" << std::flush;
music::ilog << std::setw(79) << std::setfill('.') << std::left << ">> Computing phi(3a) term" << std::endl;
Conv.convolve_Hessians(phi, {0, 0}, phi, {1, 1}, phi, {2, 2}, op::assign_to(phi3a));
Conv.convolve_Hessians(phi, {0, 1}, phi, {0, 2}, phi, {1, 2}, op::multiply_add_to(phi3a,2.0));
Conv.convolve_Hessians(phi, {1, 2}, phi, {1, 2}, phi, {0, 0}, op::subtract_from(phi3a));
Conv.convolve_Hessians(phi, {0, 2}, phi, {0, 2}, phi, {1, 1}, op::subtract_from(phi3a));
Conv.convolve_Hessians(phi, {0, 1}, phi, {0, 1}, phi, {2, 2}, op::subtract_from(phi3a));
phi3a.apply_InverseLaplacian();
music::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime() - wtime << "s" << std::endl;
music::ilog << std::setw(70) << std::setfill(' ') << std::right << "took : " << std::setw(8) << get_wtime() - wtime << "s" << std::endl;
//... 3b term ...
wtime = get_wtime();
music::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing phi(3b) term" << std::flush;
music::ilog << std::setw(71) << std::setfill('.') << std::left << ">> Computing phi(3b) term" << std::endl;
phi3b.allocate();
phi3b.FourierTransformForward(false);
Conv.convolve_SumOfHessians(phi, {0, 0}, phi2, {1, 1}, {2, 2}, op::assign_to(phi3b));
@ -456,11 +456,11 @@ int run( config_file& the_config )
Conv.convolve_Hessians(phi, {1, 2}, phi2, {1, 2}, op::multiply_add_to(phi3b,+10.0/7.0));
phi3b.apply_InverseLaplacian();
//phi3b *= 0.5; // factor 1/2 from definition of phi(3b)!
music::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime() - wtime << "s" << std::endl;
music::ilog << std::setw(70) << std::setfill(' ') << std::right << "took : " << std::setw(8) << get_wtime() - wtime << "s" << std::endl;
//... transversal term ...
wtime = get_wtime();
music::ilog << std::setw(40) << std::setfill('.') << std::left << "Computing A(3) term" << std::flush;
music::ilog << std::setw(71) << std::setfill('.') << std::left << ">> Computing A(3) term" << std::endl;
for (int idim = 0; idim < 3; ++idim)
{
// cyclic rotations of indices
@ -473,7 +473,7 @@ int run( config_file& the_config )
Conv.convolve_DifferenceOfHessians(phi2, {idimp, idimpp}, phi, {idimp, idimp}, {idimpp, idimpp}, op::subtract_from(*A3[idim]));
A3[idim]->apply_InverseLaplacian();
}
music::ilog << std::setw(20) << std::setfill(' ') << std::right << "took " << get_wtime() - wtime << "s" << std::endl;
music::ilog << std::setw(70) << std::setfill(' ') << std::right << "took : " << std::setw(8) << get_wtime() - wtime << "s" << std::endl;
}
///... scale all potentials with respective growth factors
@ -557,8 +557,11 @@ int run( config_file& the_config )
size_t IDoffset = (this_species == cosmo_species::baryon)? ((the_output_plugin->has_64bit_ids())? 1 : 1): 0 ;
// allocate particle structure and generate particle IDs
bool secondary_lattice = (this_species == cosmo_species::baryon &&
the_output_plugin->write_species_as(this_species) == output_type::particles) ? true : false;
particle_lattice_generator_ptr =
std::make_unique<particle::lattice_generator<Grid_FFT<real_t>>>( lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(),
std::make_unique<particle::lattice_generator<Grid_FFT<real_t>>>( lattice_type, secondary_lattice, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(),
bDoBaryons, IDoffset, tmp, the_config );
}
@ -566,7 +569,7 @@ int run( config_file& the_config )
if( bDoBaryons && (the_output_plugin->write_species_as( this_species ) == output_type::particles
|| the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian) )
{
bool shifted_lattice = (this_species == cosmo_species::baryon &&
bool secondary_lattice = (this_species == cosmo_species::baryon &&
the_output_plugin->write_species_as(this_species) == output_type::particles) ? true : false;
const real_t munit = the_output_plugin->mass_unit();
@ -589,7 +592,7 @@ int run( config_file& the_config )
});
if( the_output_plugin->write_species_as( this_species ) == output_type::particles ){
particle_lattice_generator_ptr->set_masses( lattice_type, shifted_lattice, 1.0, the_output_plugin->has_64bit_reals(), rho, the_config );
particle_lattice_generator_ptr->set_masses( lattice_type, secondary_lattice, 1.0, the_output_plugin->has_64bit_reals(), rho, the_config );
}else if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian ){
the_output_plugin->write_grid_data( rho, this_species, fluid_component::mass );
}

21
src/main.cc
View file

@ -43,8 +43,11 @@ bool FFTW_threads_ok = false;
int num_threads = 1;
}
size_t global_mem_high_mark, local_mem_high_mark;
#include "system_stat.hh"
#include "memory_stat.hh"
#include <exception>
#include <stdexcept>
@ -76,6 +79,8 @@ int main( int argc, char** argv )
music::logger::set_level(music::log_level::debug);
#endif
global_mem_high_mark = local_mem_high_mark = 0;
//------------------------------------------------------------------------------
// initialise MPI
//------------------------------------------------------------------------------
@ -259,13 +264,25 @@ int main( int argc, char** argv )
ic_generator::reset();
///////////////////////////////////////////////////////////////////////
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
size_t peak_mem = memory::getPeakRSS();
#if defined(USE_MPI)
size_t peak_mem_max{0};
MPI_Allreduce(&peak_mem, &peak_mem_max, 1, MPI_UNSIGNED_LONG_LONG, MPI_MAX, MPI_COMM_WORLD);
peak_mem = peak_mem_max;
#endif
if( peak_mem > (1ull<<30) )
music::ilog << "Peak memory usage was " << peak_mem /(1ull<<30) << " GBytes / task" << std::endl;
else
music::ilog << "Peak memory usage was " << peak_mem /(1ull<<20) << " MBytes / task" << std::endl;
#if defined(USE_MPI)
MPI_Barrier(MPI_COMM_WORLD);
MPI_Finalize();
#endif
music::ilog << "-------------------------------------------------------------------------------" << std::endl;
music::ilog << "Done. Have a nice day!\n" << std::endl;
return 0;

2
src/output_plugin.cc
View file

@ -74,7 +74,7 @@ std::unique_ptr<output_plugin> select_output_plugin( config_file& cf, std::uniqu
music::ilog << std::setw(32) << std::left << "Output plugin" << " : " << formatname << std::endl;
}
return std::move(the_output_plugin_creator->create( cf, pcc ));
return the_output_plugin_creator->create( cf, pcc );
}

20
src/plugins/output_gadget_hdf5.cc
View file

@ -2,17 +2,17 @@
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn
//
//
// monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//
// monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
#ifdef USE_HDF5
@ -122,8 +122,8 @@ public:
// use destructor to write header post factum
~gadget_hdf5_output_plugin()
{
if (!std::uncaught_exception())
{
if (!std::uncaught_exception())
{
HDFCreateGroup(this_fname_, "Header");
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_6array<unsigned>(header_.npart));
HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_6array<double>(header_.mass));
@ -144,6 +144,14 @@ public:
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Entropy_ICs", from_value<int>(header_.flag_entropy_instead_u));
music::ilog << "Wrote Gadget-HDF5 file(s) to " << this_fname_ << std::endl;
music::ilog << "You can use the following values in param.txt:" << std::endl;
music::ilog << "Omega0 " << header_.Omega0 << std::endl;
music::ilog << "OmegaLambda " << header_.OmegaLambda << std::endl;
music::ilog << "OmegaBaryon " << pcc_->cosmo_param_["Omega_b"] << std::endl;
music::ilog << "HubbleParam " << header_.HubbleParam << std::endl;
music::ilog << "Hubble 100.0" << std::endl;
music::ilog << "BoxSize " << header_.BoxSize << std::endl;
}
}
@ -239,4 +247,4 @@ output_plugin_creator_concrete<gadget_hdf5_output_plugin<double>> creator3("gadg
#endif
} // namespace
#endif
#endif

3
src/plugins/output_genericio.cc
View file

@ -74,8 +74,7 @@ public:
void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species)
{
double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
double boxmass = Omega_species * rhoc * std::pow(cf_.get_value<double>("setup", "BoxLength"), 3.);
double boxmass = Omega_species * munit_;
double particle_mass = boxmass / pc.get_global_num_particles();
size_t npart = pc.get_local_num_particles();

2
src/plugins/output_grafic2.cc
View file

@ -212,7 +212,7 @@ void grafic2_output_plugin::write_grid_data(const Grid_FFT<real_t> &g, const cos
unlink(file_name.c_str());
}
std::ofstream *pofs;
std::ofstream *pofs = nullptr;
// write header or seek to end of file
if (CONFIG::MPI_task_rank == 0)

4
src/plugins/random_music.cc
View file

@ -180,7 +180,7 @@ void RNG_music::parse_random_parameters(void)
//... generate some dummy seed which only depends on the level, negative so we know it's not
//... an actual seed and thus should not be used as a constraint for coarse levels
// rngseeds_.push_back( -abs((unsigned)(ltemp-i)%123+(unsigned)(ltemp+827342523521*i)%123456789) );
rngseeds_.push_back(-abs((long)(ltemp - i) % 123 + (long)(ltemp + 7342523521 * i) % 123456789));
rngseeds_.push_back(-std::abs((long)(ltemp - i) % 123 + (long)(ltemp + 7342523521 * i) % 123456789));
else
{
if (ltemp <= 0)
@ -357,4 +357,4 @@ void RNG_music::compute_random_numbers(void)
namespace
{
RNG_plugin_creator_concrete<RNG_music> creator("MUSIC1");
}
}

339
src/plugins/transfer_CLASS.cc
View file

@ -1,17 +1,17 @@
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn
//
//
// monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
//
// monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
@ -35,13 +35,16 @@
class transfer_CLASS_plugin : public TransferFunction_plugin
{
private:
using TransferFunction_plugin::cosmo_params_;
interpolated_function_1d<true, true, false> delta_c_, delta_b_, delta_n_, delta_m_, theta_c_, theta_b_, theta_n_, theta_m_;
interpolated_function_1d<true, true, false> delta_c0_, delta_b0_, delta_n0_, delta_m0_, theta_c0_, theta_b0_, theta_n0_, theta_m0_;
interpolated_function_1d<true, false, false> delta_c_, delta_b_, delta_n_, delta_t_, theta_c_, theta_b_, theta_n_, theta_t_;
interpolated_function_1d<true, false, false> delta_c0_, delta_b0_, delta_n0_, delta_t0_, theta_c0_, theta_b0_, theta_n0_, theta_t0_;
double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_;
double omega_ncdm_total_, omega_nu_total_, omega_c_;
std::vector<double> omega_ncdm_, m_ncdm_, T_ncdm_, omega_nu_;
size_t N_ncdm_, N_nu_massive_;
double zstart_, ztarget_, astart_, atarget_, kmax_, kmin_, h_, tnorm_, N_nu_massive;
ClassParams pars_;
std::unique_ptr<ClassEngine> the_ClassEngine_;
@ -58,10 +61,10 @@ private:
void init_ClassEngine(void)
{
//--- general parameters ------------------------------------------
add_class_parameter("z_max_pk", std::max(std::max(zstart_, ztarget_),199.0)); // use 1.2 as safety
add_class_parameter("P_k_max_h/Mpc", std::max(2.0,kmax_));
add_class_parameter("z_max_pk", std::max(std::max(zstart_, ztarget_), 199.0)); // use 1.2 as safety
add_class_parameter("P_k_max_h/Mpc", std::max(2.0, kmax_));
add_class_parameter("output", "dTk,vTk");
add_class_parameter("extra metric transfer functions","yes");
add_class_parameter("extra metric transfer functions", "yes");
// add_class_parameter("lensing", "no");
//--- choose gauge ------------------------------------------------
@ -72,12 +75,10 @@ private:
add_class_parameter("h", cosmo_params_.get("h"));
add_class_parameter("Omega_b", cosmo_params_.get("Omega_b"));
add_class_parameter("Omega_cdm", cosmo_params_.get("Omega_c"));
add_class_parameter("Omega_k", cosmo_params_.get("Omega_k"));
add_class_parameter("Omega_fld", 0.0);
add_class_parameter("Omega_scf", 0.0);
// add_class_parameter("fluid_equation_of_state","CLP");
// add_class_parameter("w0_fld", -1 );
// add_class_parameter("wa_fld", 0. );
@ -91,29 +92,97 @@ private:
add_class_parameter("N_ncdm", 0);
#else
add_class_parameter("N_ur", cosmo_params_.get("N_ur"));
add_class_parameter("N_ncdm", cosmo_params_.get("N_nu_massive"));
if( cosmo_params_.get("N_nu_massive") > 0 ){
std::stringstream sstr;
if( cosmo_params_.get("m_nu1") > 1e-9 ) sstr << cosmo_params_.get("m_nu1");
if( cosmo_params_.get("m_nu2") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu2");
if( cosmo_params_.get("m_nu3") > 1e-9 ) sstr << ", " << cosmo_params_.get("m_nu3");
add_class_parameter("m_ncdm", sstr.str().c_str());
// treating of massive trans-relativistic species
std::stringstream sstr_m, sstr_omega, sstr_T;
const double sum_m_nu = cosmo_params_.get("m_nu1") + cosmo_params_.get("m_nu2") + cosmo_params_.get("m_nu3");
if (N_ncdm_ > 0)
{
for (size_t i = 0; i < omega_ncdm_.size(); ++i)
{
sstr_m << m_ncdm_[i];
sstr_omega << omega_ncdm_[i];
sstr_T << T_ncdm_[i];
omega_ncdm_total_ += omega_ncdm_[i];
if( i<omega_ncdm_.size()-1 )
{
sstr_m << ", ";
sstr_omega << ", ";
sstr_T << ", ";
}
music::ilog.Print("Found extra nCDM species: m=%g eV, T=%g (T_gamma), Omega=%g", m_ncdm_[i], T_ncdm_[i], omega_ncdm_[i]);
}
music::ilog.Print("Energy density of all nCDM species: Omega = %g", omega_ncdm_total_);
omega_c_ -= omega_ncdm_total_;
}
// change above to enable
//add_class_parameter("omega_ncdm", 0.0006451439);
//add_class_parameter("m_ncdm", "0.4");
//add_class_parameter("T_ncdm", 0.71611);
N_nu_massive_ = 0;
if (cosmo_params_.get("N_nu_massive") > 0)
{
// if already nCDM species were written, need to add a comma
if (N_ncdm_ > 0)
{
sstr_m << ", ";
sstr_omega << ", ";
sstr_T << ", ";
}
if (cosmo_params_.get("m_nu1") > 1e-9)
{
const double omeganu1 = cosmo_params_.get("m_nu1") / sum_m_nu * cosmo_params_.get("Omega_nu_massive");
sstr_m << cosmo_params_.get("m_nu1");
sstr_omega << omeganu1;
sstr_T << 0.71611; // CLASS recommended value for active massive neutrinos, yields m/omega = 93.14 eV
omega_nu_.push_back(omeganu1);
omega_nu_total_ += omeganu1;
N_nu_massive_ += 1;
if (cosmo_params_.get("m_nu2") > 1e-9)
{
const double omeganu2 = cosmo_params_.get("m_nu2") / sum_m_nu * cosmo_params_.get("Omega_nu_massive");
sstr_m << ", " << cosmo_params_.get("m_nu2");
sstr_omega << ", " << omeganu2;
sstr_T << ", " << 0.71611; // CLASS recommended value for active massive neutrinos, yields m/omega = 93.14 eV
omega_nu_.push_back(omeganu2);
omega_nu_total_ += omeganu2;
N_nu_massive_ += 1;
if (cosmo_params_.get("m_nu3") > 1e-9)
{
const double omeganu3 = cosmo_params_.get("m_nu3") / sum_m_nu * cosmo_params_.get("Omega_nu_massive");
sstr_m << ", " << cosmo_params_.get("m_nu3");
sstr_omega << ", " << omeganu3;
sstr_T << ", " << 0.71611; // CLASS recommended value for active massive neutrinos, yields m/omega = 93.14 eV
omega_nu_.push_back(omeganu3);
omega_nu_total_ += omeganu3;
N_nu_massive_ += 1;
}
}
}
}
add_class_parameter("Omega_cdm", std::max(omega_c_, 1e-9));
add_class_parameter("Omega_ncdm", sstr_omega.str().c_str());
add_class_parameter("m_ncdm", sstr_m.str().c_str());
add_class_parameter("T_ncdm", sstr_T.str().c_str());
add_class_parameter("N_ncdm", cosmo_params_.get("N_nu_massive") + N_ncdm_);
#endif
//--- cosmological parameters, primordial -------------------------
add_class_parameter("P_k_ini type", "analytic_Pk");
if( cosmo_params_.get("A_s") > 0.0 ){
if (cosmo_params_.get("A_s") > 0.0)
{
add_class_parameter("A_s", cosmo_params_.get("A_s"));
}else{
}
else
{
add_class_parameter("sigma8", cosmo_params_.get("sigma_8"));
}
add_class_parameter("n_s", cosmo_params_.get("n_s"));
@ -128,7 +197,7 @@ private:
add_class_parameter("k_per_decade_for_pk", 100);
add_class_parameter("k_per_decade_for_bao", 100);
add_class_parameter("compute damping scale", "yes");
add_class_parameter("tol_perturb_integration", 1.e-8);
add_class_parameter("tol_perturbations_integration", 1.e-8);
add_class_parameter("tol_background_integration", 1e-9);
// high precision options from cl_permille.pre:
@ -147,15 +216,15 @@ private:
add_class_parameter("perturbations_verbose", class_verbosity);
add_class_parameter("transfer_verbose", class_verbosity);
add_class_parameter("primordial_verbose", class_verbosity);
add_class_parameter("spectra_verbose", class_verbosity);
add_class_parameter("nonlinear_verbose", class_verbosity);
add_class_parameter("harmonic_verbose", class_verbosity);
add_class_parameter("fourier_verbose", class_verbosity);
add_class_parameter("lensing_verbose", class_verbosity);
add_class_parameter("output_verbose", class_verbosity);
// output parameters, only needed for the control CLASS .ini file that we output
std::stringstream zlist;
if (ztarget_ == zstart_)
zlist << ztarget_ << ((ztarget_!=0.0)? ", 0.0" : "");
zlist << ztarget_ << ((ztarget_ != 0.0) ? ", 0.0" : "");
else
zlist << std::max(ztarget_, zstart_) << ", " << std::min(ztarget_, zstart_) << ", 0.0";
add_class_parameter("z_pk", zlist.str());
@ -163,23 +232,81 @@ private:
music::ilog << "Computing transfer function via ClassEngine..." << std::endl;
double wtime = get_wtime();
the_ClassEngine_ = std::move(std::make_unique<ClassEngine>(pars_, false));
the_ClassEngine_ = std::make_unique<ClassEngine>(pars_, false);
wtime = get_wtime() - wtime;
music::ilog << "CLASS took " << wtime << " s." << std::endl;
}
//! run ClassEngine with parameters set up
void run_ClassEngine(double z, std::vector<double> &k, std::vector<double> &dc, std::vector<double> &tc, std::vector<double> &db, std::vector<double> &tb,
std::vector<double> &dn, std::vector<double> &tn, std::vector<double> &dm, std::vector<double> &tm)
void run_ClassEngine(double z, int gauge, std::vector<double> &k, std::vector<double> &dc, std::vector<double> &tc, std::vector<double> &db, std::vector<double> &tb,
std::vector<double> &dn, std::vector<double> &tn, std::vector<double> &dm, std::vector<double> &tm, std::vector<double> &dt, std::vector<double> &tt,
std::vector<double> &phi_or_h_prime, std::vector<double> &psi_or_eta_prime)
{
k.clear();
dc.clear(); db.clear(); dn.clear(); dm.clear();
tc.clear(); tb.clear(); tn.clear(); tm.clear();
the_ClassEngine_->getTk(z, k, dc, db, dn, dm, tc, tb, tn, tm);
k.clear();
dc.clear();
db.clear();
dn.clear();
dm.clear();
dt.clear();
tc.clear();
tb.clear();
tn.clear();
tm.clear();
tt.clear();
const double h = cosmo_params_.get("h");
std::vector<std::vector<double>> dncdm, tncdm;
// if gauge < 0 : run class in synchronous gauge, and convert *only* velocities to conformal Newtonian gauge below, phi_or_h_prime is h_prime, psi_or_eta_prime is eta_prime
// if gauge == 0 : run class in synchronous gauge, do not convert velocities, phi_or_h_prime is h_prime, psi_or_eta_prime is eta_prime
// if gauge == 1 : run class in Newtonian gauge, do not convert velocities, phi_or_h is phi, psi_or_eta_prime is psi
// ...
const int class_gauge = (gauge < 0) ? 0 : gauge;
double fHa{0.0};
the_ClassEngine_->getTk(z, class_gauge, k, dc, db, dncdm, dm, dt, tc, tb, tncdm, tm, tt, phi_or_h_prime, psi_or_eta_prime, fHa);
// allocate memory for neutrinos
dn.assign(k.size(), 0.0);
tn.assign(k.size(), 0.0);
for (size_t index_k = 0; index_k < k.size(); ++index_k)
{
// add all neutrino species into one single massive neutrino TF
for (size_t i = N_ncdm_, j = 0; i < N_ncdm_ + N_nu_massive_; ++i, ++j)
{
dn[index_k] += omega_nu_[j] / omega_nu_total_ * dncdm[i][index_k];
tn[index_k] += omega_nu_[j] / omega_nu_total_ * tncdm[i][index_k];
}
// add all dark matter species into one single mixed-dark-matter TF (we do not support multiple DM species yet)
dc[index_k] *= omega_c_ / (omega_c_ + omega_ncdm_total_);
tc[index_k] *= omega_c_ / (omega_c_ + omega_ncdm_total_);
for (size_t i = 0; i < N_ncdm_; ++i)
{
dc[index_k] += omega_ncdm_[i] / (omega_c_ + omega_ncdm_total_) * dncdm[i][index_k];
tc[index_k] += omega_ncdm_[i] / (omega_c_ + omega_ncdm_total_) * tncdm[i][index_k];
}
}
// convert velocities to conformal Newtonian gauge if gauge < 0
if (gauge < 0)
{
for (size_t index_k = 0; index_k < k.size(); ++index_k)
{
// gauge transformation to conformal Newtonian velocity, Ma & Bertschinger eq.(27b)
double alphak2 = (phi_or_h_prime[index_k] + 6 * psi_or_eta_prime[index_k]) / 2;
tc[index_k] = (-alphak2) / fHa;
tb[index_k] = (-alphak2 + tb[index_k]) / fHa;
tn[index_k] = (-alphak2 + tn[index_k]) / fHa;
tm[index_k] = (-alphak2 + tm[index_k]) / fHa;
tt[index_k] = (-alphak2 + tt[index_k]) / fHa;
}
}
else
{
_unused(fHa);
}
const double h = cosmo_params_.get("h");
for (size_t i = 0; i < k.size(); ++i)
{
@ -190,38 +317,46 @@ private:
db[i] = -db[i] * ik2;
dn[i] = -dn[i] * ik2;
dm[i] = -dm[i] * ik2;
dt[i] = -dt[i] * ik2;
tc[i] = -tc[i] * ik2;
tb[i] = -tb[i] * ik2;
tn[i] = -tn[i] * ik2;
tm[i] = -tm[i] * ik2;
tt[i] = -tt[i] * ik2;
phi_or_h_prime[i] = -phi_or_h_prime[i] * ik2;
psi_or_eta_prime[i] = -psi_or_eta_prime[i] * ik2;
}
}
public:
explicit transfer_CLASS_plugin(config_file &cf, const cosmology::parameters& cosmo_params)
: TransferFunction_plugin(cf,cosmo_params)
explicit transfer_CLASS_plugin(config_file &cf, const cosmology::parameters &cosmo_params)
: TransferFunction_plugin(cf, cosmo_params)
{
this->tf_isnormalised_ = true;
ofs_class_input_.open("input_class_parameters.ini", std::ios::trunc);
// all cosmological parameters need to be passed through the_cosmo_calc
ztarget_ = pcf_->get_value_safe<double>("cosmology", "ztarget", 0.0);
atarget_ = 1.0 / (1.0 + ztarget_);
zstart_ = pcf_->get_value<double>("setup", "zstart");
astart_ = 1.0 / (1.0 + zstart_);
h_ = cosmo_params_["h"];
if (cosmo_params_["A_s"] > 0.0) {
if (cosmo_params_["A_s"] > 0.0)
{
music::ilog << "CLASS: Using A_s=" << cosmo_params_["A_s"] << " to normalise the transfer function." << std::endl;
}else{
}
else
{
double sigma8 = cosmo_params_["sigma_8"];
if( sigma8 < 0 ){
if (sigma8 < 0)
{
throw std::runtime_error("Need to specify either A_s or sigma_8 for CLASS plugin...");
}
music::ilog << "CLASS: Using sigma8_ =" << sigma8<< " to normalise the transfer function." << std::endl;
music::ilog << "CLASS: Using sigma8_ =" << sigma8 << " to normalise the transfer function." << std::endl;
}
// determine highest k we will need for the resolution selected
@ -229,17 +364,62 @@ public:
int nres = pcf_->get_value<double>("setup", "GridRes");
kmax_ = std::max(20.0, 2.0 * M_PI / lbox * nres / 2 * sqrt(3) * 2.0); // 120% of spatial diagonal, or k=10h Mpc-1
// initialize parameters for alternative dark matter models
omega_c_ = cosmo_params.get("Omega_c");
omega_ncdm_total_ = 0.0;
omega_nu_total_ = 0.0;
if (pcf_->contains_key("cosmology", "m_ncdm"))
{
m_ncdm_ = pcf_->get_value<std::vector<double>>("cosmology", "m_ncdm");
if (m_ncdm_.size() > 1 && !pcf_->contains_key("cosmology", "Omega_ncdm"))
{
throw std::runtime_error("Need to specify \'cosmology/Omega_ncdm\' if multiple m_ncdm are specified");
}
else
{
if (m_ncdm_.size() == 1)
{
// if m_ncdm given for one non-CDM species, default to Omega_ncdm = Omega_c if not given explicitly
omega_ncdm_.push_back(pcf_->get_value_safe<double>("cosmology", "Omega_ncdm", omega_c_));
T_ncdm_.push_back(pcf_->get_value_safe<double>("cosmology", "T_ncdm", 0.7137658555));
}
else
{
omega_ncdm_ = pcf_->get_value<std::vector<double>>("cosmology", "Omega_ncdm");
if (!pcf_->contains_key("cosmology", "T_ncdm"))
{
T_ncdm_.assign(m_ncdm_.size(), 0.7137658555); // (4/11)**(1/3)
}
else
{
T_ncdm_ = pcf_->get_value<std::vector<double>>("cosmology", "T_ncdm");
}
}
}
if (m_ncdm_.size() != omega_ncdm_.size())
{
throw std::runtime_error("Need to specify as many \'cosmology/Omega_ncdm\' as \'cosmology/m_ncdm\' (if using more than one).");
}
if (m_ncdm_.size() != T_ncdm_.size())
{
throw std::runtime_error("Need to specify as many \'cosmology/T_ncdm\' as \'cosmology/m_ncdm\' (if specifying at least one).");
}
}
N_ncdm_ = m_ncdm_.size();
// initialise CLASS and get the normalisation
this->init_ClassEngine();
double A_s_ = the_ClassEngine_->get_A_s(); // this either the input one, or the one computed from sigma8
// compute the normalisation to interface with MUSIC
double k_p = cosmo_params["k_p"] / cosmo_params["h"];
tnorm_ = std::sqrt(2.0 * M_PI * M_PI * A_s_ * std::pow(1.0 / k_p, cosmo_params["n_s"] - 1) / std::pow(2.0 * M_PI, 3.0));
// compute the transfer function at z=0 using CLASS engine
std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm;
this->run_ClassEngine(0.0, k, dc, tc, db, tb, dn, tn, dm, tm);
constexpr int gauge{-1}; // always use synchronous gauge
std::vector<double> k, dc, tc, db, tb, dn, tn, dm, tm, dt, tt, phi_or_h, psi_or_eta;
this->run_ClassEngine(0.0, gauge, k, dc, tc, db, tb, dn, tn, dm, tm, dt, tt, phi_or_h, psi_or_eta);
delta_c0_.set_data(k, dc);
theta_c0_.set_data(k, tc);
@ -247,19 +427,20 @@ public:
theta_b0_.set_data(k, tb);
delta_n0_.set_data(k, dn);
theta_n0_.set_data(k, tn);
delta_m0_.set_data(k, dm);
theta_m0_.set_data(k, tm);
delta_t0_.set_data(k, dm);
theta_t0_.set_data(k, tm);
// compute the transfer function at z=z_target using CLASS engine
this->run_ClassEngine(ztarget_, gauge, k, dc, tc, db, tb, dn, tn, dm, tm, dt, tt, phi_or_h, psi_or_eta);
// compute the transfer function at z=z_target using CLASS engine
this->run_ClassEngine(ztarget_, k, dc, tc, db, tb, dn, tn, dm, tm);
delta_c_.set_data(k, dc);
theta_c_.set_data(k, tc);
delta_b_.set_data(k, db);
theta_b_.set_data(k, tb);
delta_n_.set_data(k, dn);
theta_n_.set_data(k, tn);
delta_m_.set_data(k, dm);
theta_m_.set_data(k, tm);
delta_t_.set_data(k, dm);
theta_t_.set_data(k, tm);
kmin_ = k[0];
kmax_ = k.back();
@ -289,35 +470,49 @@ public:
{
// values at ztarget:
case delta_matter:
val = delta_m_(k); break;
val = delta_t_(k);
break;
case delta_cdm:
val = delta_c_(k); break;
val = delta_c_(k);
break;
case delta_baryon:
val = delta_b_(k); break;
val = delta_b_(k);
break;
case theta_matter:
val = theta_m_(k); break;
val = theta_t_(k);
break;
case theta_cdm:
val = theta_c_(k); break;
val = theta_c_(k);
break;
case theta_baryon:
val = theta_b_(k); break;
val = theta_b_(k);
break;
case delta_bc:
val = delta_b_(k)-delta_c_(k); break;
val = delta_b_(k) - delta_c_(k);
break;
case theta_bc:
val = theta_b_(k)-theta_c_(k); break;
val = theta_b_(k) - theta_c_(k);
break;
// values at zstart:
case delta_matter0:
val = delta_m0_(k); break;
val = delta_t0_(k);
break;
case delta_cdm0:
val = delta_c0_(k); break;
val = delta_c0_(k);
break;
case delta_baryon0:
val = delta_b0_(k); break;
val = delta_b0_(k);
break;
case theta_matter0:
val = theta_m0_(k); break;
val = theta_t0_(k);
break;
case theta_cdm0:
val = theta_c0_(k); break;
val = theta_c0_(k);
break;
case theta_baryon0:
val = theta_b0_(k); break;
val = theta_b0_(k);
break;
default:
throw std::runtime_error("Invalid type requested in transfer function evaluation");
}
@ -330,7 +525,7 @@ public:
namespace
{
TransferFunction_plugin_creator_concrete<transfer_CLASS_plugin> creator("CLASS");
TransferFunction_plugin_creator_concrete<transfer_CLASS_plugin> creator("CLASS");
}
#endif // USE_CLASS

2
src/random_plugin.cc
View file

@ -77,5 +77,5 @@ std::unique_ptr<RNG_plugin> select_RNG_plugin(config_file &cf)
music::ilog << std::setw(32) << std::left << "Random number generator plugin" << " : " << rngname << std::endl;
}
return std::move(the_RNG_plugin_creator->Create(cf));
return the_RNG_plugin_creator->Create(cf);
}

2
src/transfer_function_plugin.cc
View file

@ -75,5 +75,5 @@ std::unique_ptr<TransferFunction_plugin> select_TransferFunction_plugin(config_f
music::ilog << std::setw(32) << std::left << "Transfer function plugin" << " : " << tfname << std::endl;
}
return std::move(the_TransferFunction_plugin_creator->create(cf, cosmo_param));
return the_TransferFunction_plugin_creator->create(cf, cosmo_param);
}