diff --git a/src/plugins/output_swift.cc b/src/plugins/output_swift.cc index 636a936..aa70ed1 100644 --- a/src/plugins/output_swift.cc +++ b/src/plugins/output_swift.cc @@ -38,14 +38,15 @@ class swift_output_plugin : public output_plugin { protected: - int num_files_, num_simultaneous_writers_; + int num_ranks_, this_rank_; + int num_files_; real_t lunit_, vunit_, munit_; real_t boxsize_, hubble_param_, astart_, zstart_; bool blongids_, bdobaryons_; - std::string this_fname_; - std::array npart_, npartTotal_, npartTotalHighWord_; + std::array npart_; + std::array npartTotal_, npartTotalHighWord_; std::array mass_; double time_; @@ -54,11 +55,12 @@ public: explicit swift_output_plugin(config_file &cf, std::unique_ptr &pcc) : output_plugin(cf, pcc, "SWIFT") { + + // SWIFT uses a single file as IC */ num_files_ = 1; -#ifdef USE_MPI - // use as many output files as we have MPI tasks - MPI_Comm_size(MPI_COMM_WORLD, &num_files_); -#endif + this_rank_ = 0; + num_files_ = 1; + real_t astart = 1.0 / (1.0 + cf_.get_value("setup", "zstart")); const double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3 @@ -83,32 +85,28 @@ public: } time_ = astart; - - this_fname_ = fname_; + #ifdef USE_MPI - int thisrank = 0; - MPI_Comm_rank(MPI_COMM_WORLD, &thisrank); - if (num_files_ > 1) - this_fname_ += "." + std::to_string(thisrank); + MPI_Comm_rank(MPI_COMM_WORLD, &this_rank_); + MPI_Comm_size(MPI_COMM_WORLD, &num_ranks_); #endif - // Now that we modified the filename for MPI support, propagate the change to the super-class - // as it does operations on the file behind our back. - fname_ = this_fname_; + // Only ranks 0 writes the header + if (this_rank_ != 0) return; // delete output file if it exists - unlink(this_fname_.c_str()); + unlink(fname_.c_str()); // create output HDF5 file - HDFCreateFile(this_fname_); + HDFCreateFile(fname_); // Write UNITS header using the physical constants assumed internally by SWIFT - HDFCreateGroup(this_fname_, "Units"); - HDFWriteGroupAttribute(this_fname_, "Units", "Unit mass in cgs (U_M)", 1.98841e43); // 10^10 Msun in grams - HDFWriteGroupAttribute(this_fname_, "Units", "Unit length in cgs (U_L)", 3.08567758149e24); // 1 Mpc in cm - HDFWriteGroupAttribute(this_fname_, "Units", "Unit time in cgs (U_t)", 3.08567758149e19); // so that unit vel is 1 km/s - HDFWriteGroupAttribute(this_fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere - HDFWriteGroupAttribute(this_fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin + HDFCreateGroup(fname_, "Units"); + HDFWriteGroupAttribute(fname_, "Units", "Unit mass in cgs (U_M)", 1.98841e43); // 10^10 Msun in grams + HDFWriteGroupAttribute(fname_, "Units", "Unit length in cgs (U_L)", 3.08567758149e24); // 1 Mpc in cm + HDFWriteGroupAttribute(fname_, "Units", "Unit time in cgs (U_t)", 3.08567758149e19); // so that unit vel is 1 km/s + HDFWriteGroupAttribute(fname_, "Units", "Unit current in cgs (U_I)", 1.0); // 1 Ampere + HDFWriteGroupAttribute(fname_, "Units", "Unit temperature in cgs (U_T)", 1.0); // 1 Kelvin // Write MUSIC configuration header int order = cf_.get_value("setup", "LPTorder"); @@ -122,32 +120,32 @@ public: int do_baryonsVrel = cf_.get_value("setup", "DoBaryonVrel"); int L = cf_.get_value("setup", "GridRes"); - HDFCreateGroup(this_fname_, "ICs_parameters"); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Code", std::string("MUSIC2 - monofonIC")); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Revision", std::string(GIT_REV)); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Tag", std::string(GIT_TAG)); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Git Branch", std::string(GIT_BRANCH)); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Precision", std::string(CMAKE_PRECISION_STR)); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Convolutions", std::string(CMAKE_CONVOLVER_STR)); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "PLT", std::string(CMAKE_PLT_STR)); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "LPT Order", order); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Particle Load", load); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Transfer Function", tf); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Cosmology Parameter Set", cosmo_set); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Random Generator", rng); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Mode Fixing", do_fixing); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Mode inversion", do_invert); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Baryons", do_baryons); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Baryons Relative Velocity", do_baryonsVrel); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Grid Resolution", L); + HDFCreateGroup(fname_, "ICs_parameters"); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Code", std::string("MUSIC2 - monofonIC")); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Revision", std::string(GIT_REV)); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Tag", std::string(GIT_TAG)); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Git Branch", std::string(GIT_BRANCH)); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Precision", std::string(CMAKE_PRECISION_STR)); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Convolutions", std::string(CMAKE_CONVOLVER_STR)); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "PLT", std::string(CMAKE_PLT_STR)); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "LPT Order", order); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Particle Load", load); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Transfer Function", tf); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Cosmology Parameter Set", cosmo_set); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Random Generator", rng); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Mode Fixing", do_fixing); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Mode inversion", do_invert); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Baryons", do_baryons); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Baryons Relative Velocity", do_baryonsVrel); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Grid Resolution", L); if (tf == "CLASS") { double ztarget = cf_.get_value("cosmology", "ztarget"); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Target Redshift", ztarget); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Target Redshift", ztarget); } if (rng == "PANPHASIA") { std::string desc = cf_.get_value("random", "descriptor"); - HDFWriteGroupAttribute(this_fname_, "ICs_parameters", "Descriptor", desc); + HDFWriteGroupAttribute(fname_, "ICs_parameters", "Descriptor", desc); } } @@ -155,56 +153,27 @@ public: ~swift_output_plugin() { if (!std::uncaught_exception()) - { - HDFCreateGroup(this_fname_, "Header"); - 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(npartTotal_)); - HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_Total_HighWord", from_7array(npartTotalHighWord_)); - HDFWriteGroupAttribute(this_fname_, "Header", "NumPart_ThisFile", from_7array(npart_)); - HDFWriteGroupAttribute(this_fname_, "Header", "MassTable", from_7array(mass_)); - - HDFWriteGroupAttribute(this_fname_, "Header", "Time", from_value(time_)); - HDFWriteGroupAttribute(this_fname_, "Header", "Redshift", from_value(zstart_)); - HDFWriteGroupAttribute(this_fname_, "Header", "Flag_Entropy_ICs", from_value(0)); + { + if (this_rank_ == 0) { + + // Write Standard Gadget / SWIFT hdf5 header + HDFCreateGroup(fname_, "Header"); + HDFWriteGroupAttribute(fname_, "Header", "Dimension", 3); + HDFWriteGroupAttribute(fname_, "Header", "BoxSize", boxsize_ / hubble_param_); // in Mpc, not Mpc/h + + HDFWriteGroupAttribute(fname_, "Header", "NumPart_Total", from_7array(npartTotal_)); + HDFWriteGroupAttribute(fname_, "Header", "NumPart_Total_HighWord", from_7array(npartTotalHighWord_)); + HDFWriteGroupAttribute(fname_, "Header", "NumPart_ThisFile", from_7array(npart_)); + HDFWriteGroupAttribute(fname_, "Header", "MassTable", from_7array(mass_)); + + HDFWriteGroupAttribute(fname_, "Header", "Time", from_value(time_)); + HDFWriteGroupAttribute(fname_, "Header", "Redshift", from_value(zstart_)); + HDFWriteGroupAttribute(fname_, "Header", "Flag_Entropy_ICs", from_value(0)); + + HDFWriteGroupAttribute(fname_, "Header", "NumFilesPerSnapshot", from_value(num_files_)); - HDFWriteGroupAttribute(this_fname_, "Header", "NumFilesPerSnapshot", from_value(num_files_)); - - // write GAS internal energy and smoothing length if baryons are enabled - if( bdobaryons_ ) - { - const double gamma = cf_.get_value_safe("cosmology", "gamma", 5.0 / 3.0); - const double YHe = pcc_->cosmo_param_["YHe"]; - //const double Omega0 = cf_.get_value("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 data( npart_[0], ceint ); - HDFWriteDataset(this_fname_, "PartType0/InternalEnergy", data); - - const double h = boxsize_ / hubble_param_ / cf_.get_value("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 << "Wrote SWIFT IC file(s) to " << fname_ << std::endl; } - - music::ilog << "Wrote SWIFT IC file(s) to " << this_fname_ << std::endl; } } @@ -247,46 +216,77 @@ public: void write_particle_data(const particle::container &pc, const cosmo_species &s, double Omega_species) { int sid = get_species_idx(s); - 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()); npartTotalHighWord_[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32); - + if( pc.bhas_individual_masses_ ) mass_[sid] = 0.0; else mass_[sid] = Omega_species * munit_ / pc.get_global_num_particles(); - - HDFCreateGroup(this_fname_, std::string("PartType") + std::to_string(sid)); + + HDFCreateGroup(fname_, std::string("PartType") + std::to_string(sid)); //... write positions and velocities..... if (this->has_64bit_reals()) { - HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions64_); - HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities64_); + HDFWriteDatasetVector(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions64_); + HDFWriteDatasetVector(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities64_); } else { - HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions32_); - HDFWriteDatasetVector(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities32_); + HDFWriteDatasetVector(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Coordinates"), pc.positions32_); + HDFWriteDatasetVector(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Velocities"), pc.velocities32_); } //... write ids..... if (this->has_64bit_ids()) - HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids64_); + HDFWriteDataset(fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids64_); else - HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids32_); + HDFWriteDataset(fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids32_); //... write masses..... if( pc.bhas_individual_masses_ ){ if (this->has_64bit_reals()){ - HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass64_); + HDFWriteDataset(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass64_); }else{ - HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_); + HDFWriteDataset(fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_); } } + + // write GAS internal energy and smoothing length if baryons are enabled + if( bdobaryons_ && s == cosmo_species::baryon) { + + const double gamma = cf_.get_value_safe("cosmology", "gamma", 5.0 / 3.0); + const double YHe = pcc_->cosmo_param_["YHe"]; + //const double Omega0 = cf_.get_value("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 data( npart_[0], ceint ); + HDFWriteDataset(fname_, std::string("PartType") + std::to_string(sid) + std::string("/InternalEnergy"), data); + + const double h = boxsize_ / hubble_param_ / cf_.get_value("setup","GridRes"); + + music::ilog.Print("Swift : set initial smoothing length to mean inter-part separation: %.2f Mpc", h); + + data.assign( npart_[0], h); + HDFWriteDataset(fname_, std::string("PartType") + std::to_string(sid) + std::string("/SmoothingLength"), data); + + } } };