From a8e135c5e1000ce032d84c5bfb7defc615a96d9a Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Thu, 2 Jul 2020 01:27:53 +0200 Subject: [PATCH] added first version of perturbed masses for multi-fluid --- example.conf | 2 ++ include/output_plugin.hh | 3 ++ include/particle_container.hh | 35 +++++++++++++++--- include/particle_generator.hh | 57 ++++++++++++++++++++++++++--- src/ic_generator.cc | 59 ++++++++++++++++++++++++++++--- src/plugins/output_arepo.cc | 30 ++++++++++++---- src/plugins/output_gadget2.cc | 6 +++- src/plugins/output_gadget_hdf5.cc | 24 ++++++++++--- src/plugins/output_generic.cc | 2 ++ src/plugins/output_grafic2.cc | 6 +++- 10 files changed, 198 insertions(+), 26 deletions(-) diff --git a/example.conf b/example.conf index 073b887..e2fa656 100644 --- a/example.conf +++ b/example.conf @@ -13,6 +13,8 @@ DoBaryons = no DoFixing = yes # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x) (increases number of particles by factor!) ParticleLoad = sc +# perturb masses instead of displacements in two-component simulations? +PerturbedMasses = yes # Add a possible constraint field here: #ConstraintFieldFile = initial_conditions.h5 #ConstraintFieldName = ic_white_noise diff --git a/include/output_plugin.hh b/include/output_plugin.hh index fff657c..7f75f3a 100644 --- a/include/output_plugin.hh +++ b/include/output_plugin.hh @@ -70,6 +70,9 @@ public: //! routine to return a multiplicative factor that contains the desired velocity units for the output virtual real_t velocity_unit() const = 0; + + //! routine to return a multiplicative factor that contains critical density * box volume in desired mass units for output + virtual real_t mass_unit() const = 0; }; /*! diff --git a/include/particle_container.hh b/include/particle_container.hh index 92b683c..09f09e6 100644 --- a/include/particle_container.hh +++ b/include/particle_container.hh @@ -20,29 +20,40 @@ namespace particle{ class container { public: - std::vector positions32_, velocities32_; - std::vector positions64_, velocities64_; + std::vector positions32_, velocities32_, mass32_; + std::vector positions64_, velocities64_, mass64_; std::vector ids32_; std::vector ids64_; + bool bhas_individual_masses_; - container(){ } + container() : bhas_individual_masses_(false) { } container(const container &) = delete; - void allocate(size_t nump, bool b64reals, bool b64ids) + void allocate(size_t nump, bool b64reals, bool b64ids, bool bindividualmasses) { + bhas_individual_masses_ = bindividualmasses; + if( b64reals ){ positions64_.resize(3 * nump); velocities64_.resize(3 * nump); positions32_.clear(); velocities32_.clear(); + if( bindividualmasses ){ + mass64_.resize(nump); + mass32_.clear(); + } }else{ positions32_.resize(3 * nump); velocities32_.resize(3 * nump); positions64_.clear(); velocities64_.clear(); + if( bindividualmasses ){ + mass32_.resize(nump); + mass64_.clear(); + } } if( b64ids ){ @@ -102,6 +113,22 @@ public: ids64_[ipart] = id; } + const void* get_mass32_ptr() const{ + return reinterpret_cast( &mass32_[0] ); + } + + void set_mass32(size_t ipart, float m){ + mass32_[ipart] = m; + } + + const void* get_mass64_ptr() const{ + return reinterpret_cast( &mass64_[0] ); + } + + void set_mass64(size_t ipart, double m){ + mass64_[ipart] = m; + } + size_t get_local_num_particles(void) const { return std::max(ids32_.size(),ids64_.size()); diff --git a/include/particle_generator.hh b/include/particle_generator.hh index dc56c62..7550f57 100644 --- a/include/particle_generator.hh +++ b/include/particle_generator.hh @@ -142,7 +142,7 @@ namespace particle particle::container particles_; public: - lattice_generator(lattice lattice_type, const bool b64reals, const bool b64ids, size_t IDoffset, const field_t &field, config_file &cf) + 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) { @@ -151,7 +151,7 @@ namespace particle // 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(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); + particles_.allocate(overload * 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 * overload * field.global_size(); @@ -180,7 +180,7 @@ namespace particle else { glass_ptr_ = std::make_unique( cf, field ); - particles_.allocate(glass_ptr_->size(), b64reals, b64ids); + particles_.allocate(glass_ptr_->size(), b64reals, b64ids, false); #pragma omp parallel for for (size_t i = 0; i < glass_ptr_->size(); ++i) @@ -198,9 +198,57 @@ 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_masses(const lattice lattice_type, bool is_second_lattice, const real_t munit, const bool b64reals, field_t &field, config_file &cf) { // works only for Bravais types + if (lattice_type >= 0) + { + const size_t num_p_in_load = field.local_size(); + const real_t pmeanmass = munit / real_t(num_p_in_load); + + 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) + { + field.shift_field(second_lattice_shift[lattice_type]); + } + + // can omit first shift since zero by convention, unless shifted already above, otherwise apply relative phase shift + if (ishift > 0) + { + field.shift_field(lattice_shifts[lattice_type][ishift] - lattice_shifts[lattice_type][ishift - 1]); + } + // read out values from phase shifted field and set assoc. particle's value + const auto ipcount0 = ishift * num_p_in_load; + for (size_t i = 0, ipcount = ipcount0; 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 (b64reals) + { + particles_.set_mass64(ipcount++, pmeanmass * field.relem(i, j, k)); + } + else + { + particles_.set_mass32(ipcount++, pmeanmass * field.relem(i, j, k)); + } + } + } + } + } + }else{ + // should not happen + music::elog << "Cannot have individual particle masses for glasses!" << std::endl; + throw std::runtime_error("cannot have individual particle masses for glasses"); + } + } + + // 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) + { if (lattice_type >= 0) { const size_t num_p_in_load = field.local_size(); @@ -261,7 +309,6 @@ namespace particle void set_velocities(lattice lattice_type, bool is_second_lattice, int idim, const bool b64reals, field_t &field, config_file &cf) { - // works only for Bravais types if (lattice_type >= 0) { const size_t num_p_in_load = field.local_size(); diff --git a/src/ic_generator.cc b/src/ic_generator.cc index 11525ed..e52e1c3 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -53,7 +53,12 @@ int Run( config_file& the_config ) //-------------------------------------------------------------------------------------------------------- //! order of the LPT approximation - int LPTorder = the_config.get_value_safe("setup","LPTorder",100); + const int LPTorder = the_config.get_value_safe("setup","LPTorder",100); + + //-------------------------------------------------------------------------------------------------------- + //! use mass perturbations instead of displacements in multi-fluid ICs + const bool bPerturbedMasses = the_config.get_value_safe("setup","PerturbedMasses",true); + //-------------------------------------------------------------------------------------------------------- //! initialice particles on a bcc or fcc lattice instead of a standard sc lattice (doubles and quadruples the number of particles) @@ -401,7 +406,8 @@ int Run( config_file& the_config ) //====================================================================== //... for two-fluid ICs we need to evolve the offset field //====================================================================== - if( bDoBaryons && LPTorder > 1 ) + + if( bDoBaryons && !bPerturbedMasses && LPTorder > 1 ) { // compute phi_bc tmp.FourierTransformForward(false); @@ -457,6 +463,20 @@ int Run( config_file& the_config ) } } + // count how many species should be written out as Lagrangian particles + int num_particle_species = 0; + for( const auto& this_species : species_list ) + { + if( the_output_plugin->write_species_as( this_species ) == output_type::particles ) + ++num_particle_species; + } + + // use pertubed masses if switched on and using more than one species as particles + const bool bUsePerturbedMasses = bPerturbedMasses && (num_particle_species > 1); + + //==============================================================// + // main output loop, loop over all species that are enabled + //==============================================================// for( const auto& this_species : species_list ) { music::ilog << std::endl @@ -464,6 +484,7 @@ int Run( config_file& the_config ) const real_t C_species = (this_species == cosmo_species::baryon)? (1.0-the_cosmo_calc->cosmo_param_.f_b) : -the_cosmo_calc->cosmo_param_.f_b; + // main loop block { std::unique_ptr>> particle_lattice_generator_ptr; @@ -475,9 +496,37 @@ int Run( config_file& the_config ) // allocate particle structure and generate particle IDs particle_lattice_generator_ptr = - std::make_unique>>( lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), IDoffset, tmp, the_config ); + std::make_unique>>( lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), + bUsePerturbedMasses, IDoffset, tmp, the_config ); } + // if we use perturbed particle masses, set them + if( bUsePerturbedMasses && the_output_plugin->write_species_as( this_species ) == output_type::particles ) + { + bool shifted_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(); + + //====================================================================== + // initialise rho + //====================================================================== + Grid_FFT rho({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen}); + + wnoise.FourierTransformForward(); + rho.FourierTransformForward(false); + rho.assign_function_of_grids_kdep( [&]( auto k, auto wn ){ + return wn * the_cosmo_calc->get_amplitude_rhobc(k.norm());; + }, wnoise ); + rho.zero_DC_mode(); + rho.FourierTransformBackward(); + + rho.apply_function_r( [&]( auto prho ){ + return 1.0 + C_species * prho; + }); + + particle_lattice_generator_ptr->set_masses( lattice_type, shifted_lattice, Omega[this_species] * munit, the_output_plugin->has_64bit_reals(), rho, the_config ); + } //if( the_output_plugin->write_species_as( cosmo_species::dm ) == output_type::field_eulerian ){ if( the_output_plugin->write_species_as(this_species) == output_type::field_eulerian ) @@ -639,7 +688,7 @@ int Run( config_file& the_config ) tmp.kelem(idx) = lg.gradient(idim,tmp.get_k3(i,j,k)) * phitot + lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx); - if( bDoBaryons ){ + if( bDoBaryons && !bUsePerturbedMasses ){ vec3_t kvec = phi.get_k(i,j,k); real_t phi_bc = the_cosmo_calc->get_amplitude_phibc(kvec.norm()); tmp.kelem(idx) -= lg.gradient(idim, tmp.get_k3(i,j,k)) * wnoise.kelem(idx) * C_species * phi_bc; @@ -693,7 +742,7 @@ int Run( config_file& the_config ) tmp.kelem(idx) = lg.gradient(idim,tmp.get_k3(i,j,k)) * phitot_v + vfac3 * (lg.gradient(idimp,tmp.get_k3(i,j,k)) * A3[idimpp]->kelem(idx) - lg.gradient(idimpp,tmp.get_k3(i,j,k)) * A3[idimp]->kelem(idx)); - if( bDoBaryons && LPTorder > 1 ){ + if( bDoBaryons && !bUsePerturbedMasses && LPTorder > 1 ){ tmp.kelem(idx) += C_species * vfac * dbc3[idim]->kelem(idx); } diff --git a/src/plugins/output_arepo.cc b/src/plugins/output_arepo.cc index f2a51ed..65ef087 100644 --- a/src/plugins/output_arepo.cc +++ b/src/plugins/output_arepo.cc @@ -44,7 +44,7 @@ class gadget_hdf5_output_plugin : public output_plugin protected: int num_files_, num_simultaneous_writers_; header_t header_; - real_t lunit_, vunit_; + real_t lunit_, vunit_, munit_; bool blongids_; std::string this_fname_; double Tini_; @@ -65,8 +65,12 @@ public: MPI_Comm_size(MPI_COMM_WORLD, &num_files_); #endif 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 + lunit_ = cf_.get_value("setup", "BoxLength"); vunit_ = lunit_ / std::sqrt(astart); + munit_ = rhoc * std::pow(cf_.get_value("setup", "BoxLength"), 3); // in 1e10 h^-1 M_sol + blongids_ = cf_.get_value_safe("output", "UseLongids", false); num_simultaneous_writers_ = cf_.get_value_safe("output", "NumSimWriters", num_files_); @@ -165,6 +169,8 @@ public: real_t velocity_unit() const { return vunit_; } + real_t mass_unit() const { return munit_; } + bool has_64bit_reals() const { if (typeid(write_real_t) == typeid(double)) @@ -203,9 +209,11 @@ public: header_.npartTotal[sid] = (uint32_t)(pc.get_global_num_particles()); header_.npartTotalHighWord[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32); - double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3 - double boxmass = Omega_species * rhoc * std::pow(header_.BoxSize, 3); - header_.mass[sid] = boxmass / pc.get_global_num_particles(); + if( pc.bhas_individual_masses_ ){ + header_.mass[sid] = 0.0; + }else{ + header_.mass[sid] = Omega_species * munit_ / pc.get_global_num_particles(); + } HDFCreateGroup(this_fname_, std::string("PartType") + std::to_string(sid)); @@ -222,10 +230,20 @@ public: } //... write ids..... - if (this->has_64bit_ids()) + if (this->has_64bit_ids()){ HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/ParticleIDs"), pc.ids64_); - else + }else{ HDFWriteDataset(this_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_); + }else{ + HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_); + } + } // std::cout << ">>>A> " << header_.npart[sid] << std::endl; } diff --git a/src/plugins/output_gadget2.cc b/src/plugins/output_gadget2.cc index 0a3afbb..1b685e9 100644 --- a/src/plugins/output_gadget2.cc +++ b/src/plugins/output_gadget2.cc @@ -33,7 +33,7 @@ public: protected: int num_files_; header this_header_; - real_t lunit_, vunit_; + real_t lunit_, vunit_, munit_; bool blongids_; public: @@ -46,9 +46,11 @@ public: // use as many output files as we have MPI tasks MPI_Comm_size(MPI_COMM_WORLD, &num_files_); #endif + const double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3 real_t astart = 1.0 / (1.0 + cf_.get_value("setup", "zstart")); lunit_ = cf_.get_value("setup", "BoxLength"); vunit_ = lunit_ / std::sqrt(astart); + munit_ = rhoc * std::pow(cf_.get_value("setup", "BoxLength"), 3);; blongids_ = cf_.get_value_safe("output", "UseLongids", false); } @@ -58,6 +60,8 @@ public: real_t velocity_unit() const { return vunit_; } + real_t mass_unit() const { return munit_; } + bool has_64bit_reals() const { if (typeid(write_real_t) == typeid(double)) diff --git a/src/plugins/output_gadget_hdf5.cc b/src/plugins/output_gadget_hdf5.cc index c3b8a52..12e0f92 100644 --- a/src/plugins/output_gadget_hdf5.cc +++ b/src/plugins/output_gadget_hdf5.cc @@ -44,7 +44,7 @@ class gadget_hdf5_output_plugin : public output_plugin protected: int num_files_, num_simultaneous_writers_; header_t header_; - real_t lunit_, vunit_; + real_t lunit_, vunit_, munit_; bool blongids_; std::string this_fname_; @@ -59,8 +59,12 @@ public: MPI_Comm_size(MPI_COMM_WORLD, &num_files_); #endif 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 + lunit_ = cf_.get_value("setup", "BoxLength"); vunit_ = lunit_ / std::sqrt(astart); + munit_ = rhoc * std::pow(cf_.get_value("setup", "BoxLength"), 3); // in 1e10 h^-1 M_sol + blongids_ = cf_.get_value_safe("output", "UseLongids", false); num_simultaneous_writers_ = cf_.get_value_safe("output", "NumSimWriters", num_files_); @@ -133,6 +137,8 @@ public: real_t velocity_unit() const { return vunit_; } + real_t mass_unit() const { return munit_; } + bool has_64bit_reals() const { if (typeid(write_real_t) == typeid(double)) @@ -172,9 +178,10 @@ public: header_.npartTotal[sid] = (uint32_t)(pc.get_global_num_particles()); header_.npartTotalHighWord[sid] = (uint32_t)((pc.get_global_num_particles()) >> 32); - double rhoc = 27.7519737; // in h^2 1e10 M_sol / Mpc^3 - double boxmass = Omega_species * rhoc * std::pow(header_.BoxSize, 3); - header_.mass[sid] = boxmass / pc.get_global_num_particles(); + if( pc.bhas_individual_masses_ ) + header_.mass[sid] = 0.0; + else + header_.mass[sid] = Omega_species * munit_ / pc.get_global_num_particles(); HDFCreateGroup(this_fname_, std::string("PartType") + std::to_string(sid)); @@ -196,6 +203,15 @@ public: else HDFWriteDataset(this_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_); + }else{ + HDFWriteDataset(this_fname_, std::string("PartType") + std::to_string(sid) + std::string("/Masses"), pc.mass32_); + } + } + // std::cout << ">>>A> " << header_.npart[sid] << std::endl; } }; diff --git a/src/plugins/output_generic.cc b/src/plugins/output_generic.cc index 79c2139..49c063e 100644 --- a/src/plugins/output_generic.cc +++ b/src/plugins/output_generic.cc @@ -57,6 +57,8 @@ public: real_t position_unit() const { return 1.0; } real_t velocity_unit() const { return 1.0; } + + real_t mass_unit() const { return 1.0; } void write_grid_data(const Grid_FFT &g, const cosmo_species &s, const fluid_component &c ); }; diff --git a/src/plugins/output_grafic2.cc b/src/plugins/output_grafic2.cc index ca2a81c..b1d707d 100644 --- a/src/plugins/output_grafic2.cc +++ b/src/plugins/output_grafic2.cc @@ -31,7 +31,7 @@ private: protected: header header_; - real_t lunit_, vunit_; + real_t lunit_, vunit_, munit_; uint32_t levelmin_; bool bhavebaryons_; std::vector data_buf_; @@ -82,6 +82,8 @@ public: lunit_ = boxlength; vunit_ = boxlength; + munit_ = 1.0; + #warning need to fix mass unit for grafic output // create directory structure dirname_ = this->fname_; @@ -110,6 +112,8 @@ public: real_t velocity_unit() const { return vunit_; } + real_t mass_unit() const { return munit_; } + void write_grid_data(const Grid_FFT &g, const cosmo_species &s, const fluid_component &c); };