From 9ea23ea581755365cdee3a0de5a14328942aac36 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 4 Feb 2022 00:33:54 +0100 Subject: [PATCH 1/6] first commit of the masked particle load for the upcoming Virgo consortium simulations --- example.conf | 3 +- include/particle_generator.hh | 268 ++++++++++++++++++++++++++++++++-- src/ic_generator.cc | 12 +- 3 files changed, 263 insertions(+), 20 deletions(-) diff --git a/example.conf b/example.conf index 8c72613..9781316 100644 --- a/example.conf +++ b/example.conf @@ -20,7 +20,8 @@ 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 = glass' then specify here where to load the glass distribution from # GlassFileName = glass128.hdf5 diff --git a/include/particle_generator.hh b/include/particle_generator.hh index bc72ae7..3cb3cdb 100644 --- a/include/particle_generator.hh +++ b/include/particle_generator.hh @@ -29,13 +29,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> 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>> lattice_shifts = { // first shift must always be zero! (otherwise set_positions and set_velocities break) @@ -149,11 +156,40 @@ namespace particle private: particle::container particles_; + size_t global_num_particles_; + + static constexpr int masksize_ = 2; + std::array particle_type_mask_; + + inline int get_mask_value( const vec3_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]; + } 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) + /** + * @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) { - if (lattice_type != lattice_glass) + // 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 { music::wlog << "Glass ICs will currently be incorrect due to disabled ghost zone updates! "; @@ -163,6 +199,9 @@ namespace particle 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, 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,7 +227,82 @@ 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) ) + { + // mask everywehere 0, except the last element + for( auto& m : particle_type_mask_) m = 0; + particle_type_mask_[7] = 1; + }else{ + // mask everywehere 0 + 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 ) { glass_ptr_ = std::make_unique( cf, field ); particles_.allocate(glass_ptr_->size(), b64reals, b64ids, false); @@ -206,14 +320,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(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 +349,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 +385,49 @@ 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(); + } + + const real_t pmeanmass = munit / global_num_particles_; + + bool bmass_negative = false; + auto mean_pm = field.mean() * pmeanmass; + auto std_pm = 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) + { + if( this->get_mask_value(field.get_cell_idx_3d(i,j,k)) != lattice_index ) continue; + + // get + const auto pmass = pmeanmass * field.relem(i, j, k); + + // check for negative mass + bmass_negative |= pmass<0.0; + + // set + if (b64reals) particles_.set_mass64(ipcount++, pmass); + else particles_.set_mass32(ipcount++, pmass); + } + } + } + + // 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 +436,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 +468,7 @@ namespace particle { for (size_t k = 0; k < field.size(2); ++k) { - auto pos = field.template get_unit_r_shifted(i, j, k, lattice_shifts[lattice_type][ishift] + (is_second_lattice ? second_lattice_shift[lattice_type] : vec3_t{real_t(0.), real_t(0.), real_t(0.)})); + auto pos = field.template get_unit_r_shifted(i, j, k, lattice_shifts[lattice_type][ishift] + (lattice_index==1 ? second_lattice_shift[lattice_type] : vec3_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 +482,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(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 +534,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 +577,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 ); diff --git a/src/ic_generator.cc b/src/ic_generator.cc index 1561b1a..9e5e907 100644 --- a/src/ic_generator.cc +++ b/src/ic_generator.cc @@ -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] @@ -536,8 +537,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>>( lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), + std::make_unique>>( lattice_type, secondary_lattice, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), bDoBaryons, IDoffset, tmp, the_config ); } @@ -545,7 +549,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(); @@ -568,7 +572,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 ); } From 7155897fd8521b76412c065a1f8bbd54ecb8f61c Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 4 Feb 2022 02:09:47 +0100 Subject: [PATCH 2/6] added new ghost zone routines (that should be robust with MPI) --- include/grid_ghosts.hh | 209 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 209 insertions(+) create mode 100644 include/grid_ghosts.hh diff --git a/include/grid_ghosts.hh b/include/grid_ghosts.hh new file mode 100644 index 0000000..be31ede --- /dev/null +++ b/include/grid_ghosts.hh @@ -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 . +#pragma once + +#include +#include + +#include + +#include +#include + +template +struct grid_with_ghosts +{ + using data_t = typename grid_t::data_t; + using vec3 = std::array; + + 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 boundary_left_, boundary_right_; + std::vector local0starts_; + const grid_t &gridref; + size_t nx_, ny_, nz_, nzp_; + + + //... determine communication offsets + std::vector 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 req; + MPI_Request temp_req; + + if( have_right ){ + for( int itask=0; itask= 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(), 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(), 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= 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(), 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(), 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 &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); + } +}; + From a1e2158f17eaba4f798ba74e16c14a79ba02f9f7 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 4 Feb 2022 02:10:10 +0100 Subject: [PATCH 3/6] subtract mean so that masked particles have correct statistics --- include/particle_generator.hh | 54 ++++++++++++++++++++++++++++++++--- 1 file changed, 50 insertions(+), 4 deletions(-) diff --git a/include/particle_generator.hh b/include/particle_generator.hh index 3cb3cdb..febd478 100644 --- a/include/particle_generator.hh +++ b/include/particle_generator.hh @@ -17,6 +17,7 @@ #pragma once #include +#include #include #if defined(USE_HDF5) @@ -169,6 +170,31 @@ namespace particle return particle_type_mask_[sig]; } + template< typename ggrid_t > + inline real_t get_mean_mask_value( const ggrid_t& gg_field, const vec3_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 gg_field(field); + const real_t pmeanmass = munit / global_num_particles_; bool bmass_negative = false; - auto mean_pm = field.mean() * pmeanmass; - auto std_pm = field.std() * pmeanmass; + 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) @@ -409,10 +437,14 @@ namespace particle { 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; + 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); + const auto pmass = pmeanmass * (field.relem(i, j, k) - mean_mask); // check for negative mass bmass_negative |= pmass<0.0; @@ -420,9 +452,23 @@ namespace particle // 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; From a015ec0f42929cb7c91943170ec3b258b8d61cf8 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 4 Feb 2022 21:16:55 +0100 Subject: [PATCH 4/6] cleaned up some compiler warnings --- .../panphasia_ho/uniform_rand_threefry4x64.c | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/external/panphasia_ho/uniform_rand_threefry4x64.c b/external/panphasia_ho/uniform_rand_threefry4x64.c index 2dcf8dc..683ba85 100644 --- a/external/panphasia_ho/uniform_rand_threefry4x64.c +++ b/external/panphasia_ho/uniform_rand_threefry4x64.c @@ -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(); From 819be9d48d7d68a09ebc67c563dc88320b2cd713 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Wed, 9 Feb 2022 17:56:56 +0100 Subject: [PATCH 5/6] moved warning that was in incorrect place --- include/particle_generator.hh | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/particle_generator.hh b/include/particle_generator.hh index febd478..37bb946 100644 --- a/include/particle_generator.hh +++ b/include/particle_generator.hh @@ -217,8 +217,6 @@ namespace particle if (lattice_type >= 0) // These are the Bravais lattices { - music::wlog << "Glass ICs will currently be incorrect due to disabled ghost zone updates! "; - // 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): @@ -330,6 +328,8 @@ namespace particle } 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( cf, field ); particles_.allocate(glass_ptr_->size(), b64reals, b64ids, false); From 8281f0bed612245f1243a669f1a992644c865f08 Mon Sep 17 00:00:00 2001 From: Oliver Hahn Date: Fri, 25 Mar 2022 18:17:10 +0100 Subject: [PATCH 6/6] added customization to particle masking --- example.conf | 3 +++ include/particle_generator.hh | 30 ++++++++++++++++++++++-------- 2 files changed, 25 insertions(+), 8 deletions(-) diff --git a/example.conf b/example.conf index 9781316..76f3a1d 100644 --- a/example.conf +++ b/example.conf @@ -23,6 +23,9 @@ ParticleLoad = sc # particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc # (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 # GlassTiles = 1 diff --git a/include/particle_generator.hh b/include/particle_generator.hh index 37bb946..b0c1a1b 100644 --- a/include/particle_generator.hh +++ b/include/particle_generator.hh @@ -38,11 +38,11 @@ namespace particle lattice_rsc = 3, // RSC: refined simple cubic }; - const std::vector> 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> 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>> lattice_shifts = { @@ -266,11 +266,25 @@ namespace particle // 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 = 0; - particle_type_mask_[7] = 1; + 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 everywehere 0 + // mask everywhere 0, all particle locations are occupied by CDM for( auto& m : particle_type_mask_) m = 0; }