mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
Merged in update_swift_plugin_mpi (pull request #21)
Update MPI-flavour of swift plugin Approved-by: Oliver Hahn
This commit is contained in:
commit
3fe061b10a
4 changed files with 369 additions and 119 deletions
|
@ -17,7 +17,7 @@ Currently supported features (the list is growing, so check back):
|
||||||
|
|
||||||
- Multiple Einstein-Boltzmann modules: direct interface with [CLASS](https://lesgourg.github.io/class_public/class.html), file input from CAMB, and fitting formulae (Eisenstein&Hu).
|
- Multiple Einstein-Boltzmann modules: direct interface with [CLASS](https://lesgourg.github.io/class_public/class.html), file input from CAMB, and fitting formulae (Eisenstein&Hu).
|
||||||
|
|
||||||
- Multiple output modules for RAMSES, Arepo, Gadget-2/3, and HACC (courtesy M.Buehlmann) via plugins (Swift and Nyx are next).
|
- Multiple output modules for RAMSES, Arepo, Gadget-2/3, SWIFT, and HACC (courtesy M.Buehlmann) via plugins (Nyx is next).
|
||||||
|
|
||||||
- Multiple random number modules (MUSIC1,NGenIC,Panphasia,...) (A new MUSIC2 module is in development)
|
- Multiple random number modules (MUSIC1,NGenIC,Panphasia,...) (A new MUSIC2 module is in development)
|
||||||
|
|
||||||
|
|
|
@ -144,7 +144,8 @@ NumThreads = 8
|
||||||
# format = genericio
|
# format = genericio
|
||||||
# filename = ics_hacc
|
# filename = ics_hacc
|
||||||
|
|
||||||
##> SWIFT compatible HDF5 format
|
##> SWIFT compatible HDF5 format. Format broadly similar to gadget_hdf5 but in a single
|
||||||
|
##> file even when using MPI. No h-factors and no sqrt(a)-factor for the velocities.
|
||||||
# format = SWIFT
|
# format = SWIFT
|
||||||
# filename = ics_swift.hdf5
|
# filename = ics_swift.hdf5
|
||||||
|
|
||||||
|
|
|
@ -886,6 +886,178 @@ inline void HDFWriteDatasetVector( const std::string Filename, const std::string
|
||||||
H5Fclose( HDF_FileID );
|
H5Fclose( HDF_FileID );
|
||||||
}
|
}
|
||||||
|
|
||||||
|
template< typename T >
|
||||||
|
inline void HDFCreateEmptyDataset( const std::string Filename, const std::string ObjName, const size_t num_particles, const bool filter = false)
|
||||||
|
{
|
||||||
|
|
||||||
|
hid_t
|
||||||
|
HDF_FileID,
|
||||||
|
HDF_DatasetID,
|
||||||
|
HDF_DataspaceID,
|
||||||
|
HDF_Type,
|
||||||
|
HDF_Prop;
|
||||||
|
|
||||||
|
hsize_t HDF_Dims;
|
||||||
|
|
||||||
|
HDF_Prop = H5Pcreate(H5P_DATASET_CREATE);
|
||||||
|
|
||||||
|
if (filter)
|
||||||
|
{
|
||||||
|
// 1MB chunking
|
||||||
|
hsize_t HDF_Dims[1] = {1024 * 1024 / sizeof(T)};
|
||||||
|
H5Pset_chunk(HDF_Prop, 1, HDF_Dims);
|
||||||
|
|
||||||
|
// md5 checksum
|
||||||
|
H5Pset_fletcher32(HDF_Prop);
|
||||||
|
}
|
||||||
|
|
||||||
|
HDF_FileID = H5Fopen( Filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
|
||||||
|
|
||||||
|
HDF_Type = GetDataType<T>();
|
||||||
|
|
||||||
|
HDF_Dims = (hsize_t) (num_particles);
|
||||||
|
HDF_DataspaceID = H5Screate_simple(1, &HDF_Dims, NULL);
|
||||||
|
HDF_DatasetID = H5Dcreate( HDF_FileID, ObjName.c_str(), HDF_Type,
|
||||||
|
HDF_DataspaceID, HDF_Prop );
|
||||||
|
H5Dclose( HDF_DatasetID );
|
||||||
|
H5Sclose( HDF_DataspaceID );
|
||||||
|
H5Pclose( HDF_Prop );
|
||||||
|
|
||||||
|
H5Fclose( HDF_FileID );
|
||||||
|
}
|
||||||
|
|
||||||
|
template< typename T >
|
||||||
|
inline void HDFCreateEmptyDatasetVector( const std::string Filename, const std::string ObjName, const size_t num_particles, const bool filter = false)
|
||||||
|
{
|
||||||
|
|
||||||
|
hid_t
|
||||||
|
HDF_FileID,
|
||||||
|
HDF_DatasetID,
|
||||||
|
HDF_DataspaceID,
|
||||||
|
HDF_Type,
|
||||||
|
HDF_Prop;
|
||||||
|
|
||||||
|
hsize_t HDF_Dims[2];
|
||||||
|
|
||||||
|
HDF_Prop = H5Pcreate(H5P_DATASET_CREATE);
|
||||||
|
|
||||||
|
if (filter)
|
||||||
|
{
|
||||||
|
// ~1MB chunking
|
||||||
|
hsize_t HDF_Dims[2] = {256 * 1024 / sizeof(T), 3};
|
||||||
|
H5Pset_chunk(HDF_Prop, 2, HDF_Dims);
|
||||||
|
|
||||||
|
// md5 checksum
|
||||||
|
H5Pset_fletcher32(HDF_Prop);
|
||||||
|
}
|
||||||
|
|
||||||
|
HDF_FileID = H5Fopen( Filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
|
||||||
|
|
||||||
|
HDF_Type = GetDataType<T>();
|
||||||
|
|
||||||
|
HDF_Dims[0] = (hsize_t) (num_particles);
|
||||||
|
HDF_Dims[1] = (hsize_t) 3;
|
||||||
|
|
||||||
|
HDF_DataspaceID = H5Screate_simple(2, HDF_Dims, NULL);
|
||||||
|
HDF_DatasetID = H5Dcreate( HDF_FileID, ObjName.c_str(), HDF_Type,
|
||||||
|
HDF_DataspaceID, HDF_Prop );
|
||||||
|
H5Dclose( HDF_DatasetID );
|
||||||
|
H5Sclose( HDF_DataspaceID );
|
||||||
|
H5Pclose( HDF_Prop );
|
||||||
|
|
||||||
|
H5Fclose( HDF_FileID );
|
||||||
|
}
|
||||||
|
|
||||||
|
template< typename T >
|
||||||
|
inline void HDFWriteDatasetChunk( const std::string Filename, const std::string ObjName, const std::vector<T> &Data, const size_t offset )
|
||||||
|
{
|
||||||
|
|
||||||
|
hid_t
|
||||||
|
HDF_FileID,
|
||||||
|
HDF_DatasetID,
|
||||||
|
HDF_MemDataspaceID,
|
||||||
|
HDF_FileDataspaceID,
|
||||||
|
HDF_Type;
|
||||||
|
|
||||||
|
hsize_t HDF_Dims,
|
||||||
|
HDF_Shape,
|
||||||
|
HDF_Offset;
|
||||||
|
|
||||||
|
HDF_FileID = H5Fopen( Filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
|
||||||
|
|
||||||
|
HDF_Type = GetDataType<T>();
|
||||||
|
|
||||||
|
HDF_Dims = (hsize_t)Data.size();
|
||||||
|
HDF_MemDataspaceID = H5Screate_simple(1, &HDF_Dims, NULL);
|
||||||
|
|
||||||
|
HDF_Shape = (hsize_t)Data.size();
|
||||||
|
HDF_Offset = (hsize_t)offset;
|
||||||
|
|
||||||
|
HDF_DatasetID = H5Dopen( HDF_FileID, ObjName.c_str() );
|
||||||
|
|
||||||
|
HDF_FileDataspaceID = H5Dget_space(HDF_DatasetID);
|
||||||
|
|
||||||
|
H5Sselect_hyperslab(HDF_FileDataspaceID, H5S_SELECT_SET, &HDF_Offset, NULL, &HDF_Shape, NULL);
|
||||||
|
|
||||||
|
H5Dwrite( HDF_DatasetID, HDF_Type, HDF_MemDataspaceID, HDF_FileDataspaceID, H5P_DEFAULT, &Data[0] );
|
||||||
|
|
||||||
|
H5Dclose( HDF_DatasetID );
|
||||||
|
H5Sclose( HDF_MemDataspaceID );
|
||||||
|
H5Sclose( HDF_FileDataspaceID );
|
||||||
|
H5Fclose( HDF_FileID );
|
||||||
|
}
|
||||||
|
|
||||||
|
template< typename T >
|
||||||
|
inline void HDFWriteDatasetVectorChunk( const std::string Filename, const std::string ObjName, const std::vector<T> &Data, const size_t offset )
|
||||||
|
{
|
||||||
|
|
||||||
|
hid_t
|
||||||
|
HDF_FileID,
|
||||||
|
HDF_DatasetID,
|
||||||
|
HDF_MemDataspaceID,
|
||||||
|
HDF_FileDataspaceID,
|
||||||
|
HDF_Type;
|
||||||
|
|
||||||
|
hsize_t HDF_Dims[2],
|
||||||
|
HDF_Shape[2],
|
||||||
|
HDF_Offset[2];
|
||||||
|
|
||||||
|
HDF_FileID = H5Fopen( Filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT );
|
||||||
|
|
||||||
|
HDF_Type = GetDataType<T>();
|
||||||
|
|
||||||
|
HDF_Dims[0] = (hsize_t)(Data.size()/3);
|
||||||
|
HDF_Dims[1] = 3;
|
||||||
|
|
||||||
|
HDF_MemDataspaceID = H5Screate_simple(2, HDF_Dims, NULL);
|
||||||
|
|
||||||
|
if( Data.size() % 3 != 0 ){
|
||||||
|
std::cerr << " - Warning: Trying to write vector data in HDFWriteDatasetVector\n"
|
||||||
|
<< " but array length not divisible by 3!\n\n";
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
HDF_Shape[0] = (hsize_t)(Data.size()/3);
|
||||||
|
HDF_Shape[1] = (hsize_t)3;
|
||||||
|
HDF_Offset[0] = (hsize_t)offset;
|
||||||
|
HDF_Offset[1] = (hsize_t)0;
|
||||||
|
|
||||||
|
HDF_DatasetID = H5Dopen( HDF_FileID, ObjName.c_str() );
|
||||||
|
|
||||||
|
HDF_FileDataspaceID = H5Dget_space(HDF_DatasetID);
|
||||||
|
|
||||||
|
H5Sselect_hyperslab(HDF_FileDataspaceID, H5S_SELECT_SET, HDF_Offset, NULL, HDF_Shape, NULL);
|
||||||
|
|
||||||
|
H5Dwrite( HDF_DatasetID, HDF_Type, HDF_MemDataspaceID, HDF_FileDataspaceID, H5P_DEFAULT, &Data[0] );
|
||||||
|
|
||||||
|
H5Dclose( HDF_DatasetID );
|
||||||
|
H5Sclose( HDF_MemDataspaceID );
|
||||||
|
H5Sclose( HDF_FileDataspaceID );
|
||||||
|
H5Fclose( HDF_FileID );
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
inline void HDFCreateGroup( const std::string Filename, const std::string GroupName )
|
inline void HDFCreateGroup( const std::string Filename, const std::string GroupName )
|
||||||
{
|
{
|
||||||
hid_t HDF_FileID, HDF_GroupID;
|
hid_t HDF_FileID, HDF_GroupID;
|
||||||
|
|
|
@ -38,29 +38,32 @@ class swift_output_plugin : public output_plugin
|
||||||
{
|
{
|
||||||
|
|
||||||
protected:
|
protected:
|
||||||
int num_files_, num_simultaneous_writers_;
|
int num_ranks_, this_rank_;
|
||||||
|
int num_files_;
|
||||||
|
|
||||||
real_t lunit_, vunit_, munit_;
|
real_t lunit_, vunit_, munit_;
|
||||||
real_t boxsize_, hubble_param_, astart_, zstart_;
|
real_t boxsize_, hubble_param_, astart_, zstart_;
|
||||||
bool blongids_, bdobaryons_;
|
bool blongids_, bdobaryons_;
|
||||||
std::string this_fname_;
|
|
||||||
|
|
||||||
std::array<uint32_t,7> npart_, npartTotal_, npartTotalHighWord_;
|
std::array<uint64_t,7> npart_;
|
||||||
|
std::array<uint32_t,7> npartTotal_, npartTotalHighWord_;
|
||||||
std::array<double,7> mass_;
|
std::array<double,7> mass_;
|
||||||
double time_;
|
double time_;
|
||||||
|
double ceint_, h_;
|
||||||
|
|
||||||
public:
|
public:
|
||||||
//! constructor
|
//! constructor
|
||||||
explicit swift_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
|
explicit swift_output_plugin(config_file &cf, std::unique_ptr<cosmology::calculator> &pcc)
|
||||||
: output_plugin(cf, pcc, "SWIFT")
|
: output_plugin(cf, pcc, "SWIFT")
|
||||||
{
|
{
|
||||||
|
|
||||||
|
// SWIFT uses a single file as IC */
|
||||||
num_files_ = 1;
|
num_files_ = 1;
|
||||||
#ifdef USE_MPI
|
this_rank_ = 0;
|
||||||
// use as many output files as we have MPI tasks
|
num_ranks_ = 1;
|
||||||
MPI_Comm_size(MPI_COMM_WORLD, &num_files_);
|
|
||||||
#endif
|
|
||||||
real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
|
real_t astart = 1.0 / (1.0 + cf_.get_value<double>("setup", "zstart"));
|
||||||
const double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3
|
const double rhoc = 27.7536609198; // in h^2 1e10 M_sol / Mpc^3 assuming SWIFT's internal constants
|
||||||
|
|
||||||
hubble_param_ = pcc->cosmo_param_["h"];
|
hubble_param_ = pcc->cosmo_param_["h"];
|
||||||
zstart_ = cf_.get_value<double>("setup","zstart");
|
zstart_ = cf_.get_value<double>("setup","zstart");
|
||||||
|
@ -84,31 +87,49 @@ public:
|
||||||
|
|
||||||
time_ = astart;
|
time_ = astart;
|
||||||
|
|
||||||
this_fname_ = fname_;
|
|
||||||
#ifdef USE_MPI
|
#ifdef USE_MPI
|
||||||
int thisrank = 0;
|
MPI_Comm_rank(MPI_COMM_WORLD, &this_rank_);
|
||||||
MPI_Comm_rank(MPI_COMM_WORLD, &thisrank);
|
MPI_Comm_size(MPI_COMM_WORLD, &num_ranks_);
|
||||||
if (num_files_ > 1)
|
|
||||||
this_fname_ += "." + std::to_string(thisrank);
|
|
||||||
#endif
|
#endif
|
||||||
|
|
||||||
// Now that we modified the filename for MPI support, propagate the change to the super-class
|
if (bdobaryons_) {
|
||||||
// as it does operations on the file behind our back.
|
|
||||||
fname_ = this_fname_;
|
const double gamma = cf_.get_value_safe<double>("cosmology", "gamma", 5.0 / 3.0);
|
||||||
|
const double YHe = pcc_->cosmo_param_["YHe"];
|
||||||
|
const double omegab = pcc_->cosmo_param_["Omega_b"];
|
||||||
|
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;
|
||||||
|
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));
|
||||||
|
ceint_ = 1.3806e-16 / 1.6726e-24 * Tini * npol / mu / unitv / unitv;
|
||||||
|
|
||||||
|
music::ilog.Print("Swift : Calculated initial gas temperature: %.2f K/mu", Tini / mu);
|
||||||
|
music::ilog.Print("Swift : set initial internal energy to %.2e km^2/s^2", ceint_);
|
||||||
|
|
||||||
|
h_ = boxsize_ / hubble_param_ / cf_.get_value<double>("setup","GridRes");
|
||||||
|
music::ilog.Print("Swift : set initial smoothing length to mean inter-part separation: %.2f Mpc", h_);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Only ranks 0 writes the header
|
||||||
|
if (this_rank_ != 0) return;
|
||||||
|
|
||||||
// delete output file if it exists
|
// delete output file if it exists
|
||||||
unlink(this_fname_.c_str());
|
unlink(fname_.c_str());
|
||||||
|
|
||||||
// create output HDF5 file
|
// create output HDF5 file
|
||||||
HDFCreateFile(this_fname_);
|
HDFCreateFile(fname_);
|
||||||
|
|
||||||
// Write UNITS header using the physical constants assumed internally by SWIFT
|
// Write UNITS header using the physical constants assumed internally by SWIFT
|
||||||
HDFCreateGroup(this_fname_, "Units");
|
HDFCreateGroup(fname_, "Units");
|
||||||
HDFWriteGroupAttribute(this_fname_, "Units", "Unit mass in cgs (U_M)", 1.98841e43); // 10^10 Msun in grams
|
HDFWriteGroupAttribute(fname_, "Units", "Unit mass in cgs (U_M)", 1.98841e43); // 10^10 Msun in grams
|
||||||
HDFWriteGroupAttribute(this_fname_, "Units", "Unit length in cgs (U_L)", 3.08567758149e24); // 1 Mpc in cm
|
HDFWriteGroupAttribute(fname_, "Units", "Unit length in cgs (U_L)", 3.08567758149e24); // 1 Mpc in cm
|
||||||
HDFWriteGroupAttribute(this_fname_, "Units", "Unit time in cgs (U_t)", 3.08567758149e19); // so that unit vel is 1 km/s
|
HDFWriteGroupAttribute(fname_, "Units", "Unit time in cgs (U_t)", 3.08567758149e19); // so that unit vel is 1 km/s
|
||||||
HDFWriteGroupAttribute(this_fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere
|
HDFWriteGroupAttribute(fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere
|
||||||
HDFWriteGroupAttribute(this_fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin
|
HDFWriteGroupAttribute(fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin
|
||||||
|
|
||||||
// Write MUSIC configuration header
|
// Write MUSIC configuration header
|
||||||
int order = cf_.get_value<int>("setup", "LPTorder");
|
int order = cf_.get_value<int>("setup", "LPTorder");
|
||||||
|
@ -122,32 +143,32 @@ public:
|
||||||
int do_baryonsVrel = cf_.get_value<bool>("setup", "DoBaryonVrel");
|
int do_baryonsVrel = cf_.get_value<bool>("setup", "DoBaryonVrel");
|
||||||
int L = cf_.get_value<int>("setup", "GridRes");
|
int L = cf_.get_value<int>("setup", "GridRes");
|
||||||
|
|
||||||
HDFCreateGroup(this_fname_, "ICs_parameters");
|
HDFCreateGroup(fname_, "ICs_parameters");
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Code", std::string("MUSIC2 - monofonIC"));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Code", std::string("MUSIC2 - monofonIC"));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Revision", std::string(GIT_REV));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Revision", std::string(GIT_REV));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Tag", std::string(GIT_TAG));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Tag", std::string(GIT_TAG));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Branch", std::string(GIT_BRANCH));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Branch", std::string(GIT_BRANCH));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Precision", std::string(CMAKE_PRECISION_STR));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Precision", std::string(CMAKE_PRECISION_STR));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Convolutions", std::string(CMAKE_CONVOLVER_STR));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Convolutions", std::string(CMAKE_CONVOLVER_STR));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "PLT", std::string(CMAKE_PLT_STR));
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "PLT", std::string(CMAKE_PLT_STR));
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "LPT Order", order);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "LPT Order", order);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Particle Load", load);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Particle Load", load);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Transfer Function", tf);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Transfer Function", tf);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Cosmology Parameter Set", cosmo_set);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Cosmology Parameter Set", cosmo_set);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Random Generator", rng);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Random Generator", rng);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Mode Fixing", do_fixing);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Mode Fixing", do_fixing);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Mode inversion", do_invert);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Mode inversion", do_invert);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Baryons", do_baryons);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Baryons", do_baryons);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Baryons Relative Velocity", do_baryonsVrel);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Baryons Relative Velocity", do_baryonsVrel);
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Grid Resolution", L);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Grid Resolution", L);
|
||||||
|
|
||||||
if (tf == "CLASS") {
|
if (tf == "CLASS") {
|
||||||
double ztarget = cf_.get_value<double>("cosmology", "ztarget");
|
double ztarget = cf_.get_value<double>("cosmology", "ztarget");
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Target Redshift", ztarget);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Target Redshift", ztarget);
|
||||||
}
|
}
|
||||||
if (rng == "PANPHASIA") {
|
if (rng == "PANPHASIA") {
|
||||||
std::string desc = cf_.get_value<std::string>("random", "descriptor");
|
std::string desc = cf_.get_value<std::string>("random", "descriptor");
|
||||||
HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Descriptor", desc);
|
HDFWriteGroupAttribute(fname_, "ICs_parameters", "Descriptor", desc);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -156,55 +177,28 @@ public:
|
||||||
{
|
{
|
||||||
if (!std::uncaught_exception())
|
if (!std::uncaught_exception())
|
||||||
{
|
{
|
||||||
HDFCreateGroup(this_fname_, "Header");
|
if (this_rank_ == 0) {
|
||||||
HDFWriteGroupAttribute(fname_, "Header", "Dimension", 3);
|
|
||||||
HDFWriteGroupAttribute(fname_, "Header", "BoxSize", boxsize_ / hubble_param_); // in Mpc, not Mpc/h
|
|
||||||
|
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total", from_7array<unsigned>(npartTotal_));
|
// Write Standard Gadget / SWIFT hdf5 header
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_7array<unsigned>(npartTotalHighWord_));
|
HDFCreateGroup(fname_, "Header");
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_7array<unsigned>(npart_));
|
HDFWriteGroupAttribute(fname_, "Header", "Dimension", 3);
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_7array<double>(mass_));
|
HDFWriteGroupAttribute(fname_, "Header", "BoxSize", boxsize_ / hubble_param_); // in Mpc, not Mpc/h
|
||||||
|
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "Time", from_value<double>(time_));
|
HDFWriteGroupAttribute(fname_, "Header", "NumPart_Total", from_7array<unsigned>(npartTotal_));
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "Redshift", from_value<double>(zstart_));
|
HDFWriteGroupAttribute(fname_, "Header", "NumPart_Total_HighWord", from_7array<unsigned>(npartTotalHighWord_));
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Entropy_ICs", from_value<int>(0));
|
HDFWriteGroupAttribute(fname_, "Header", "NumPart_ThisFile", from_7array<uint64_t>(npart_));
|
||||||
|
HDFWriteGroupAttribute(fname_, "Header", "MassTable", from_7array<double>(mass_));
|
||||||
|
|
||||||
HDFWriteGroupAttribute(this_fname_, "Header", "NumFilesPerSnapshot", from_value<int>(num_files_));
|
HDFWriteGroupAttribute(fname_, "Header", "Time", from_value<double>(time_));
|
||||||
|
HDFWriteGroupAttribute(fname_, "Header", "Redshift", from_value<double>(zstart_));
|
||||||
|
HDFWriteGroupAttribute(fname_, "Header", "Flag_Entropy_ICs", from_value<int>(0));
|
||||||
|
|
||||||
// write GAS internal energy and smoothing length if baryons are enabled
|
HDFWriteGroupAttribute(fname_, "Header", "NumFilesPerSnapshot", from_value<int>(num_files_));
|
||||||
if( bdobaryons_ )
|
|
||||||
{
|
|
||||||
const double gamma = cf_.get_value_safe<double>("cosmology", "gamma", 5.0 / 3.0);
|
|
||||||
const double YHe = pcc_->cosmo_param_["YHe"];
|
|
||||||
//const double Omega0 = cf_.get_value<double>("cosmology", "Omega_m");
|
|
||||||
const double omegab = pcc_->cosmo_param_["Omega_b"];
|
|
||||||
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;
|
|
||||||
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;
|
|
||||||
|
|
||||||
music::ilog.Print("Swift : Calculated initial gas temperature: %.2f K/mu", Tini / mu);
|
|
||||||
music::ilog.Print("Swift : set initial internal energy to %.2e km^2/s^2", ceint);
|
|
||||||
|
|
||||||
std::vector<write_real_t> data( npart_[0], ceint );
|
|
||||||
HDFWriteDataset(this_fname_, "PartType0/InternalEnergy", data);
|
|
||||||
|
|
||||||
const double h = boxsize_ / hubble_param_ / cf_.get_value<double>("setup","GridRes");
|
|
||||||
|
|
||||||
music::ilog.Print("Swift : set initial smoothing length to mean inter-part separation: %.2f Mpc", h);
|
|
||||||
|
|
||||||
data.assign( npart_[0], h);
|
|
||||||
HDFWriteDataset(this_fname_, "PartType0/SmoothingLength", data);
|
|
||||||
|
|
||||||
|
music::ilog << "Done writing SWIFT IC file to " << fname_ << std::endl;
|
||||||
|
music::ilog << "Note that the IC file does not contain any h-factors nor any extra sqrt(a)-factors for the velocities" << std::endl;
|
||||||
|
music::ilog << "The SWIFT parameters 'InitialConditions:cleanup_h_factors' and 'InitialConditions:cleanup_velocity_factors' should hence *NOT* be set!" << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
music::ilog << "Wrote SWIFT IC file(s) to " << this_fname_ << std::endl;
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
@ -247,46 +241,129 @@ public:
|
||||||
void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species)
|
void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species)
|
||||||
{
|
{
|
||||||
int sid = get_species_idx(s);
|
int sid = get_species_idx(s);
|
||||||
|
|
||||||
assert(sid != -1);
|
assert(sid != -1);
|
||||||
|
|
||||||
npart_[sid] = (pc.get_local_num_particles());
|
npart_[sid] = (pc.get_global_num_particles());
|
||||||
npartTotal_[sid] = (uint32_t)(pc.get_global_num_particles());
|
npartTotal_[sid] = (uint32_t)(pc.get_global_num_particles());
|
||||||
npartTotalHighWord_[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32);
|
npartTotalHighWord_[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32);
|
||||||
|
|
||||||
|
std::array<double,7> particle_masses;
|
||||||
if( pc.bhas_individual_masses_ )
|
if( pc.bhas_individual_masses_ )
|
||||||
mass_[sid] = 0.0;
|
particle_masses[sid] = 0.0;
|
||||||
else
|
else
|
||||||
mass_[sid] = Omega_species * munit_ / pc.get_global_num_particles();
|
particle_masses[sid] = Omega_species * munit_ / pc.get_global_num_particles();
|
||||||
|
|
||||||
HDFCreateGroup(this_fname_, std::string("PartType") + std::to_string(sid));
|
const size_t global_num_particles = pc.get_global_num_particles();
|
||||||
|
|
||||||
//... write positions and velocities.....
|
// start by creating the full empty datasets in the file
|
||||||
if (this->has_64bit_reals())
|
if (this_rank_ == 0) {
|
||||||
{
|
|
||||||
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions64_);
|
HDFCreateGroup(fname_, std::string("PartType") + std::to_string(sid));
|
||||||
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities64_);
|
|
||||||
}
|
if (this->has_64bit_reals())
|
||||||
else
|
{
|
||||||
{
|
HDFCreateEmptyDatasetVector<double>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), global_num_particles, true);
|
||||||
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions32_);
|
HDFCreateEmptyDatasetVector<double>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), global_num_particles, true);
|
||||||
HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities32_);
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
HDFCreateEmptyDatasetVector<float>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), global_num_particles, true);
|
||||||
|
HDFCreateEmptyDatasetVector<float>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), global_num_particles, true);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (this->has_64bit_ids())
|
||||||
|
HDFCreateEmptyDataset<uint64_t>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), global_num_particles, true);
|
||||||
|
else
|
||||||
|
HDFCreateEmptyDataset<uint32_t>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), global_num_particles, true);
|
||||||
|
|
||||||
|
if (this->has_64bit_reals())
|
||||||
|
HDFCreateEmptyDataset<double>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), global_num_particles, true);
|
||||||
|
else
|
||||||
|
HDFCreateEmptyDataset<float>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), global_num_particles, true);
|
||||||
|
|
||||||
|
if( bdobaryons_ && s == cosmo_species::baryon) {
|
||||||
|
|
||||||
|
// note: despite this being a constant array we still need to handle it in a distributed way
|
||||||
|
HDFCreateEmptyDataset<write_real_t>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/InternalEnergy"), global_num_particles, true);
|
||||||
|
HDFCreateEmptyDataset<write_real_t>(fname_, std::string("PartType") + std::to_string(sid) + std::string("/SmoothingLength"), global_num_particles, true);
|
||||||
|
}
|
||||||
|
|
||||||
|
music::ilog << "Created empty arrays for PartType" << std::to_string(sid) << " into file " << fname_ << "." << std::endl;
|
||||||
}
|
}
|
||||||
|
|
||||||
//... write ids.....
|
// compute each rank's offset in the global array
|
||||||
if (this->has_64bit_ids())
|
const size_t n_local = pc.get_local_num_particles();
|
||||||
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids64_);
|
size_t offset = 0;
|
||||||
else
|
#ifdef USE_MPI
|
||||||
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids32_);
|
MPI_Exscan(&n_local, &offset, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, MPI_COMM_WORLD);
|
||||||
|
#endif
|
||||||
|
|
||||||
//... write masses.....
|
// now each node writes its own chunk in a round-robin fashion, appending at the end of the currently existing data
|
||||||
if( pc.bhas_individual_masses_ ){
|
for (int rank = 0; rank < num_ranks_; ++rank) {
|
||||||
if (this->has_64bit_reals()){
|
|
||||||
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass64_);
|
#ifdef USE_MPI
|
||||||
}else{
|
// wait until the initialisation or the previous rank in the loop is done
|
||||||
HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_);
|
MPI_Barrier(MPI_COMM_WORLD);
|
||||||
|
#endif
|
||||||
|
|
||||||
|
if (rank == this_rank_) {
|
||||||
|
|
||||||
|
//... write positions and velocities.....
|
||||||
|
if (this->has_64bit_reals())
|
||||||
|
{
|
||||||
|
HDFWriteDatasetVectorChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions64_, offset);
|
||||||
|
HDFWriteDatasetVectorChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities64_, offset);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
HDFWriteDatasetVectorChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions32_, offset);
|
||||||
|
HDFWriteDatasetVectorChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities32_, offset);
|
||||||
|
}
|
||||||
|
|
||||||
|
//... write ids.....
|
||||||
|
if (this->has_64bit_ids())
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids64_, offset);
|
||||||
|
else
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids32_, offset);
|
||||||
|
|
||||||
|
//... write masses.....
|
||||||
|
if( pc.bhas_individual_masses_ )
|
||||||
|
{
|
||||||
|
if (this->has_64bit_reals())
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass64_, offset);
|
||||||
|
else
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_, offset);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
std::vector<write_real_t> data( pc.get_local_num_particles(), particle_masses[sid]);
|
||||||
|
|
||||||
|
if (this->has_64bit_reals())
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), data, offset);
|
||||||
|
else
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), data, offset);
|
||||||
|
}
|
||||||
|
|
||||||
|
// write GAS internal energy and smoothing length if baryons are enabled
|
||||||
|
if(bdobaryons_ && s == cosmo_species::baryon) {
|
||||||
|
|
||||||
|
std::vector<write_real_t> data( pc.get_local_num_particles(), ceint_ );
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/InternalEnergy"), data, offset);
|
||||||
|
|
||||||
|
data.assign( pc.get_local_num_particles(), h_);
|
||||||
|
HDFWriteDatasetChunk(fname_, std::string("PartType") + std::to_string(sid) + std::string("/SmoothingLength"), data, offset);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
if (this_rank_ == 0) {
|
||||||
|
music::ilog << "Rank " << rank << " wrote its PartType" << std::to_string(sid) << " data to the IC file." << std::endl;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
#ifdef USE_MPI
|
||||||
|
// end with a barrier to make sure everyone is done before the destructor does its job
|
||||||
|
MPI_Barrier(MPI_COMM_WORLD);
|
||||||
|
#endif
|
||||||
}
|
}
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
Loading…
Reference in a new issue