2020-08-21 17:03:33 +02:00
// This file is part of monofonIC (MUSIC2)
// A software package to generate ICs for cosmological simulations
// Copyright (C) 2020 by Oliver Hahn
//
// monofonIC is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// monofonIC is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
2019-10-19 12:56:43 +02:00
# pragma once
2020-04-04 23:59:13 +02:00
# include <math/vec3.hh>
2022-02-04 02:10:10 +01:00
# include <grid_ghosts.hh>
2020-05-03 00:14:47 +02:00
# include <grid_interpolate.hh>
2019-11-01 04:47:02 +01:00
2020-05-03 04:20:12 +02:00
# if defined(USE_HDF5)
# include "HDF_IO.hh"
# endif
2020-05-04 00:49:11 +02:00
namespace particle
{
using vec3 = std : : array < real_t , 3 > ;
2019-10-19 12:56:43 +02:00
2020-05-04 00:49:11 +02:00
enum lattice
{
2022-02-04 00:33:54 +01:00
lattice_masked = - 2 , // masked lattices are SC lattices with a specifiable mask leaving out particles
lattice_glass = - 1 , // glass: needs a glass file
2020-05-04 00:49:11 +02:00
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
} ;
2019-11-01 04:47:02 +01:00
2022-03-25 18:17:10 +01:00
// const std::vector<std::vector<bool>> lattice_masks =
// {
// // mask from Richings et al. https://arxiv.org/pdf/2005.14495.pdf
// {true,true,true,true,true,true,true,false},
// };
2022-02-04 00:33:54 +01:00
2020-05-04 00:49:11 +02:00
const std : : vector < std : : vector < vec3_t < real_t > > > lattice_shifts =
{
// first shift must always be zero! (otherwise set_positions and set_velocities break)
2020-10-29 19:50:05 +01:00
/* SC : */ { { real_t { 0.0 } , real_t { 0.0 } , real_t { 0.0 } } } ,
/* BCC: */ { { real_t { 0.0 } , real_t { 0.0 } , real_t { 0.0 } } , { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.5 } } } ,
/* FCC: */ { { real_t { 0.0 } , real_t { 0.0 } , real_t { 0.0 } } , { real_t { 0.0 } , real_t { 0.5 } , real_t { 0.5 } } , { real_t { 0.5 } , real_t { 0.0 } , real_t { 0.5 } } , { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.0 } } } ,
/* RSC: */ { { real_t { 0.0 } , real_t { 0.0 } , real_t { 0.0 } } , { real_t { 0.0 } , real_t { 0.0 } , real_t { 0.5 } } , { real_t { 0.0 } , real_t { 0.5 } , real_t { 0.0 } } , { real_t { 0.0 } , real_t { 0.5 } , real_t { 0.5 } } , { real_t { 0.5 } , real_t { 0.0 } , real_t { 0.0 } } , { real_t { 0.5 } , real_t { 0.0 } , real_t { 0.5 } } , { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.0 } } , { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.5 } } } ,
2020-05-04 00:49:11 +02:00
} ;
2019-10-19 12:56:43 +02:00
2020-05-04 00:49:11 +02:00
const std : : vector < vec3_t < real_t > > second_lattice_shift =
{
2020-10-29 19:50:05 +01:00
/* SC : */ { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.5 } } , // this corresponds to CsCl lattice
/* BCC: */ { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.0 } } , // is there a diatomic lattice with BCC base?!?
/* FCC: */ { real_t { 0.5 } , real_t { 0.5 } , real_t { 0.5 } } , // this corresponds to NaCl lattice
// /* FCC: */ {real_t{0.25}, real_t{0.25}, real_t{0.25}}, // this corresponds to Zincblende/GaAs lattice
/* RSC: */ { real_t { 0.25 } , real_t { 0.25 } , real_t { 0.25 } } ,
2020-05-04 00:49:11 +02:00
} ;
2020-01-24 15:00:32 +01:00
2020-05-04 00:49:11 +02:00
template < typename field_t >
class lattice_generator
2020-05-03 04:20:12 +02:00
{
2020-05-04 00:49:11 +02:00
protected :
struct glass
{
using data_t = typename field_t : : data_t ;
size_t num_p , off_p ;
grid_interpolate < 1 , field_t > interp_ ;
std : : vector < vec3 > glass_posr ;
glass ( config_file & cf , const field_t & field )
: num_p ( 0 ) , off_p ( 0 ) , interp_ ( field )
{
std : : vector < real_t > glass_pos ;
real_t lglassbox = 1.0 ;
std : : string glass_fname = cf . get_value < std : : string > ( " setup " , " GlassFileName " ) ;
size_t ntiles = cf . get_value < size_t > ( " setup " , " GlassTiles " ) ;
2020-05-03 04:20:12 +02:00
# if defined(USE_HDF5)
2020-05-04 00:49:11 +02:00
HDFReadGroupAttribute ( glass_fname , " Header " , " BoxSize " , lglassbox ) ;
HDFReadDataset ( glass_fname , " /PartType1/Coordinates " , glass_pos ) ;
# else
throw std : : runtime_error ( " Class lattice requires HDF5 support. Enable and recompile. " ) ;
# endif
size_t np_in_file = glass_pos . size ( ) / 3 ;
2020-05-03 16:15:40 +02:00
# if defined(USE_MPI)
2020-05-04 00:49:11 +02:00
num_p = np_in_file * ntiles * ntiles * ntiles / MPI : : get_size ( ) ;
off_p = MPI : : get_rank ( ) * num_p ;
2020-05-03 16:15:40 +02:00
# else
2020-05-04 00:49:11 +02:00
num_p = np_in_file * ntiles * ntiles * ntiles ;
off_p = 0 ;
2020-05-03 16:15:40 +02:00
# endif
2020-05-03 04:20:12 +02:00
2020-05-04 01:13:07 +02:00
music : : ilog < < " Glass file contains " < < np_in_file < < " particles. " < < std : : endl ;
2020-05-04 00:49:11 +02:00
glass_posr . assign ( num_p , { 0.0 , 0.0 , 0.0 } ) ;
2020-05-03 04:20:12 +02:00
2020-05-04 00:49:11 +02:00
std : : array < real_t , 3 > ng ( { real_t ( field . n_ [ 0 ] ) , real_t ( field . n_ [ 1 ] ) , real_t ( field . n_ [ 2 ] ) } ) ;
2020-05-04 02:09:42 +02:00
# pragma omp parallel for
2020-05-04 00:49:11 +02:00
for ( size_t i = 0 ; i < num_p ; + + i )
{
size_t idxpart = off_p + i ;
size_t idx_in_glass = idxpart % np_in_file ;
size_t idxtile = idxpart / np_in_file ;
size_t tile_z = idxtile % ( ntiles * ntiles ) ;
size_t tile_y = ( ( idxtile - tile_z ) / ntiles ) % ntiles ;
size_t tile_x = ( ( ( idxtile - tile_z ) / ntiles ) - tile_y ) / ntiles ;
glass_posr [ i ] [ 0 ] = std : : fmod ( ( glass_pos [ 3 * idx_in_glass + 0 ] / lglassbox + real_t ( tile_x ) ) / ntiles * ng [ 0 ] + ng [ 0 ] , ng [ 0 ] ) ;
glass_posr [ i ] [ 1 ] = std : : fmod ( ( glass_pos [ 3 * idx_in_glass + 1 ] / lglassbox + real_t ( tile_y ) ) / ntiles * ng [ 1 ] + ng [ 1 ] , ng [ 1 ] ) ;
glass_posr [ i ] [ 2 ] = std : : fmod ( ( glass_pos [ 3 * idx_in_glass + 2 ] / lglassbox + real_t ( tile_z ) ) / ntiles * ng [ 2 ] + ng [ 2 ] , ng [ 2 ] ) ;
}
# if defined(USE_MPI)
interp_ . domain_decompose_pos ( glass_posr ) ;
num_p = glass_posr . size ( ) ;
std : : vector < size_t > all_num_p ( MPI : : get_size ( ) , 0 ) ;
MPI_Allgather ( & num_p , 1 , MPI_UNSIGNED_LONG_LONG , & all_num_p [ 0 ] , 1 , MPI_UNSIGNED_LONG_LONG , MPI_COMM_WORLD ) ;
off_p = 0 ;
for ( int itask = 0 ; itask < = MPI : : get_rank ( ) ; + + itask ) {
off_p + = all_num_p [ itask ] ;
}
2020-05-03 04:20:12 +02:00
# endif
2020-05-04 00:49:11 +02:00
}
2019-10-19 12:56:43 +02:00
2020-05-04 01:13:07 +02:00
void update_ghosts ( const field_t & field )
{
interp_ . update_ghosts ( field ) ;
}
2020-05-04 00:49:11 +02:00
data_t get_at ( const vec3 & x ) const noexcept
{
return interp_ . get_cic_at ( x ) ;
}
size_t size ( ) const noexcept
{
return num_p ;
2020-05-03 00:14:47 +02:00
}
2020-01-24 15:00:32 +01:00
2020-05-04 00:49:11 +02:00
size_t offset ( ) const noexcept
{
return off_p ;
2020-05-03 00:14:47 +02:00
}
2020-05-04 00:49:11 +02:00
} ;
std : : unique_ptr < glass > glass_ptr_ ;
private :
particle : : container particles_ ;
2022-02-04 00:33:54 +01:00
size_t global_num_particles_ ;
static constexpr int masksize_ = 2 ;
std : : array < int , masksize_ * masksize_ * masksize_ > particle_type_mask_ ;
inline int get_mask_value ( const vec3_t < size_t > & global_idx_3d ) const
{
int sig = ( ( global_idx_3d [ 0 ] % masksize_ ) * masksize_
+ global_idx_3d [ 1 ] % masksize_ ) * masksize_
+ global_idx_3d [ 2 ] % masksize_ ;
return particle_type_mask_ [ sig ] ;
}
2020-05-04 00:49:11 +02:00
2022-02-04 02:10:10 +01:00
template < typename ggrid_t >
inline real_t get_mean_mask_value ( const ggrid_t & gg_field , const vec3_t < size_t > & global_idx_3d , size_t i , size_t j , size_t k , int lattice_index ) const
{
ptrdiff_t ox = ( ptrdiff_t ) i - ( ptrdiff_t ) ( global_idx_3d [ 0 ] % masksize_ ) ;
ptrdiff_t oy = ( ptrdiff_t ) j - ( ptrdiff_t ) ( global_idx_3d [ 1 ] % masksize_ ) ;
ptrdiff_t oz = ( ptrdiff_t ) k - ( ptrdiff_t ) ( global_idx_3d [ 2 ] % masksize_ ) ;
size_t count_full { 0 } , count_masked { 0 } ;
real_t mean_full { 0.0 } , mean_masked { 0.0 } ;
for ( size_t i = 0 ; i < masksize_ ; + + i ) {
for ( size_t j = 0 ; j < masksize_ ; + + j ) {
for ( size_t k = 0 ; k < masksize_ ; + + k ) {
mean_full + = gg_field . relem ( ox + i , oy + j , oz + k ) ;
+ + count_full ;
if ( particle_type_mask_ [ 4 * i + 2 * j + k ] = = lattice_index ) {
mean_masked + = gg_field . relem ( ox + i , oy + j , oz + k ) ;
+ + count_masked ;
}
}
}
}
mean_full / = count_full ;
mean_masked / = count_masked ;
return mean_masked - mean_full ;
}
2020-05-04 00:49:11 +02:00
public :
2022-02-04 00:33:54 +01:00
/**
* @ 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 64 bit floating point shall be used to store positions / velocities / masses
* @ param b64ids Boolean whether 64 bit 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 )
2020-05-04 00:49:11 +02:00
{
2022-02-04 00:33:54 +01:00
// 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
2020-05-04 00:49:11 +02:00
{
// number of modes present in the field
const size_t num_p_in_load = field . local_size ( ) ;
// unless SC lattice is used, particle number is a multiple of the number of modes (=num_p_in_load):
const size_t overload = 1ull < < std : : max < int > ( 0 , lattice_type ) ; // 1 for sc, 2 for bcc, 4 for fcc, 8 for rsc
// allocate memory for all local particles
2020-07-02 01:27:53 +02:00
particles_ . allocate ( overload * num_p_in_load , b64reals , b64ids , bwithmasses ) ;
2022-02-04 00:33:54 +01:00
// set the global number of particles for this lattice_type and lattice_index
global_num_particles_ = field . global_size ( ) * overload ;
2020-05-04 00:49:11 +02:00
// set particle IDs to the Lagrangian coordinate (1D encoded) with additionally the field shift encoded as well
2020-05-04 02:09:42 +02:00
2020-05-24 13:01:14 +02:00
IDoffset = IDoffset * overload * field . global_size ( ) ;
2020-05-25 18:55:40 +02:00
for ( size_t i = 0 , ipcount = 0 ; i < field . rsize ( 0 ) ; + + i )
2020-05-04 00:49:11 +02:00
{
2020-05-25 18:55:40 +02:00
for ( size_t j = 0 ; j < field . rsize ( 1 ) ; + + j )
2020-05-04 00:49:11 +02:00
{
2020-05-25 18:55:40 +02:00
for ( size_t k = 0 ; k < field . rsize ( 2 ) ; + + k , + + ipcount )
2020-05-04 00:49:11 +02:00
{
for ( size_t iload = 0 ; iload < overload ; + + iload )
{
if ( b64ids )
{
particles_ . set_id64 ( ipcount + iload * num_p_in_load , IDoffset + overload * field . get_cell_idx_1d ( i , j , k ) + iload ) ;
}
else
{
particles_ . set_id32 ( ipcount + iload * num_p_in_load , IDoffset + overload * field . get_cell_idx_1d ( i , j , k ) + iload ) ;
}
}
2020-05-03 00:14:47 +02:00
}
2019-11-01 12:03:02 +01:00
}
2019-11-01 04:47:02 +01:00
}
}
2022-02-04 00:33:54 +01:00
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 ) )
{
2022-03-25 18:17:10 +01:00
unsigned mask_type = cf . get_value_safe ( " setup " , " ParticleMaskType " , 3 ) ;
2022-02-04 00:33:54 +01:00
// mask everywehere 0, except the last element
2022-03-25 18:17:10 +01:00
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)
2022-02-04 00:33:54 +01:00
} else {
2022-03-25 18:17:10 +01:00
// mask everywhere 0, all particle locations are occupied by CDM
2022-02-04 00:33:54 +01:00
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 )
2020-05-04 00:49:11 +02:00
{
2022-02-09 17:56:56 +01:00
music : : wlog < < " Glass ICs will currently be incorrect due to disabled ghost zone updates! " ;
2020-05-04 00:49:11 +02:00
glass_ptr_ = std : : make_unique < glass > ( cf , field ) ;
2020-07-02 01:27:53 +02:00
particles_ . allocate ( glass_ptr_ - > size ( ) , b64reals , b64ids , false ) ;
2020-05-03 04:20:12 +02:00
2020-05-04 02:09:42 +02:00
# pragma omp parallel for
2020-05-04 00:49:11 +02:00
for ( size_t i = 0 ; i < glass_ptr_ - > size ( ) ; + + i )
{
if ( b64ids )
{
particles_ . set_id64 ( i , IDoffset + i + glass_ptr_ - > offset ( ) ) ;
}
else
{
particles_ . set_id32 ( i , IDoffset + i + glass_ptr_ - > offset ( ) ) ;
}
}
}
2022-02-04 00:33:54 +01:00
music : : ilog < < " Created Particles [ " < < lattice_index < < " ] : " < < global_num_particles_ < < std : : endl ;
2020-05-03 04:20:12 +02:00
}
2020-05-04 00:49:11 +02:00
// invalidates field, phase shifted to unspecified position after return
2022-02-04 00:33:54 +01:00
void set_masses ( const lattice lattice_type , int lattice_index , const real_t munit , const bool b64reals , field_t & field , config_file & cf )
2020-05-04 00:49:11 +02:00
{
2022-02-04 00:33:54 +01:00
// works only for Bravais types and masked type
assert ( lattice_type > = 0 | | lattice_type = = lattice_masked ) ;
if ( lattice_type > = 0 ) // Bravais lattices
2020-07-02 01:27:53 +02:00
{
2022-02-04 00:33:54 +01:00
if ( lattice_index > 1 | | lattice_index < 0 ) {
music : : elog < < " For Bravais lattice type, lattice index must be 0 or 1. " < < std : : endl ;
abort ( ) ;
}
2020-07-05 14:26:53 +02:00
const size_t overload = 1ull < < std : : max < int > ( 0 , lattice_type ) ; // 1 for sc, 2 for bcc, 4 for fcc, 8 for rsc
2020-07-02 01:27:53 +02:00
const size_t num_p_in_load = field . local_size ( ) ;
2020-07-05 14:26:53 +02:00
const real_t pmeanmass = munit / real_t ( field . global_size ( ) * overload ) ;
2020-07-02 01:27:53 +02:00
2021-10-05 23:00:26 +02:00
bool bmass_negative = false ;
auto mean_pm = field . mean ( ) * pmeanmass ;
auto std_pm = field . std ( ) * pmeanmass ;
2020-07-02 01:27:53 +02:00
for ( int ishift = 0 ; ishift < ( 1 < < lattice_type ) ; + + ishift )
{
// if we are dealing with the secondary lattice, apply a global shift
2022-02-04 00:33:54 +01:00
if ( ishift = = 0 & & lattice_index > 0 )
2020-07-02 01:27:53 +02:00
{
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 )
{
2021-10-05 23:00:26 +02:00
// 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 ) ;
2020-07-02 01:27:53 +02:00
}
}
}
}
2021-10-05 23:00:26 +02:00
// 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 ;
2022-02-04 00:33:54 +01:00
} 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 ( ) ;
}
2022-02-04 02:10:10 +01:00
grid_with_ghosts < 1 , true , true , field_t > gg_field ( field ) ;
2022-02-04 00:33:54 +01:00
const real_t pmeanmass = munit / global_num_particles_ ;
bool bmass_negative = false ;
2022-02-04 02:10:10 +01:00
real_t mean_pm = 0.0 ; //field.mean() * pmeanmass;
real_t std_pm = 0.0 ; //field.std() * pmeanmass;
2022-02-04 00:33:54 +01:00
// 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 )
{
2022-02-04 02:10:10 +01:00
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 ) ;
2022-02-04 00:33:54 +01:00
// get
2022-02-04 02:10:10 +01:00
const auto pmass = pmeanmass * ( field . relem ( i , j , k ) - mean_mask ) ;
2021-10-05 23:00:26 +02:00
2022-02-04 00:33:54 +01:00
// check for negative mass
bmass_negative | = pmass < 0.0 ;
// set
if ( b64reals ) particles_ . set_mass64 ( ipcount + + , pmass ) ;
else particles_ . set_mass32 ( ipcount + + , pmass ) ;
2022-02-04 02:10:10 +01:00
// statistics
mean_pm + = pmass ;
std_pm + = pmass * pmass ;
2022-02-04 00:33:54 +01:00
}
}
}
2022-02-04 02:10:10 +01:00
# 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 ;
2022-02-04 00:33:54 +01:00
// 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 ;
2020-07-02 01:27:53 +02:00
} 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
2022-02-04 00:33:54 +01:00
void set_positions ( const lattice lattice_type , int lattice_index , int idim , real_t lunit , const bool b64reals , field_t & field , config_file & cf )
2020-07-02 01:27:53 +02:00
{
2022-02-04 00:33:54 +01:00
if ( lattice_type > = 0 ) // Bravais lattice
2020-05-04 00:49:11 +02:00
{
2022-02-04 00:33:54 +01:00
if ( lattice_index > 1 | | lattice_index < 0 ) {
music : : elog < < " For Bravais lattice type, lattice index must be 0 or 1. " < < std : : endl ;
abort ( ) ;
}
2020-05-04 00:49:11 +02:00
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
2022-02-04 00:33:54 +01:00
if ( ishift = = 0 & & lattice_index = = 1 )
2020-05-04 00:49:11 +02:00
{
field . shift_field ( second_lattice_shift [ lattice_type ] ) ;
}
2020-05-03 04:20:12 +02:00
2020-05-04 00:49:11 +02:00
// 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 )
{
2022-02-04 00:33:54 +01:00
auto pos = field . template get_unit_r_shifted < real_t > ( i , j , k , lattice_shifts [ lattice_type ] [ ishift ] + ( lattice_index = = 1 ? second_lattice_shift [ lattice_type ] : vec3_t < real_t > { real_t ( 0. ) , real_t ( 0. ) , real_t ( 0. ) } ) ) ;
2020-05-04 00:49:11 +02:00
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 ) ) ;
}
}
}
}
}
}
2022-02-04 00:33:54 +01:00
else if ( lattice_type = = lattice_masked )
{
if ( field . global_size ( ) % 8 ! = 0 ) {
music : : elog < < " For masked lattice type, linear field resolution must be a multiple of two. " < < std : : endl ;
abort ( ) ;
}
if ( lattice_index > 1 | | lattice_index < 0 ) {
music : : elog < < " For masked lattice type, lattice index must be 0 or 1. " < < std : : endl ;
abort ( ) ;
}
for ( size_t i = 0 , ipcount = 0 ; i < field . size ( 0 ) ; + + i )
{
for ( size_t j = 0 ; j < field . size ( 1 ) ; + + j )
{
for ( size_t k = 0 ; k < field . size ( 2 ) ; + + k )
{
if ( this - > get_mask_value ( field . get_cell_idx_3d ( i , j , k ) ) ! = lattice_index ) continue ;
// get position (in box units) of the current cell of 3d array 'field'
auto pos = field . template get_unit_r < real_t > ( i , j , k ) ;
// add the displacement to get the particle position
if ( b64reals ) {
particles_ . set_pos64 ( ipcount + + , idim , pos [ idim ] * lunit + field . relem ( i , j , k ) ) ;
} else {
particles_ . set_pos32 ( ipcount + + , idim , pos [ idim ] * lunit + field . relem ( i , j , k ) ) ;
}
}
}
}
}
2020-05-04 00:49:11 +02:00
else
{
2020-05-04 01:13:07 +02:00
glass_ptr_ - > update_ghosts ( field ) ;
2020-05-04 02:09:42 +02:00
# pragma omp parallel for
2020-05-04 01:13:07 +02:00
for ( size_t i = 0 ; i < glass_ptr_ - > size ( ) ; + + i )
2020-05-04 00:49:11 +02:00
{
2020-05-04 01:13:07 +02:00
auto pos = glass_ptr_ - > glass_posr [ i ] ;
real_t disp = glass_ptr_ - > get_at ( pos ) ;
2020-05-04 00:49:11 +02:00
if ( b64reals )
{
particles_ . set_pos64 ( i , idim , pos [ idim ] / field . n_ [ idim ] * lunit + disp ) ;
}
else
{
particles_ . set_pos32 ( i , idim , pos [ idim ] / field . n_ [ idim ] * lunit + disp ) ;
}
}
2020-05-03 04:20:12 +02:00
}
}
2019-11-01 04:47:02 +01:00
2022-02-04 00:33:54 +01:00
void set_velocities ( lattice lattice_type , int lattice_index , int idim , const bool b64reals , field_t & field , config_file & cf )
2020-05-04 00:49:11 +02:00
{
2022-02-04 00:33:54 +01:00
if ( lattice_type > = 0 ) // Bravais lattice
2020-05-04 00:49:11 +02:00
{
2022-02-04 00:33:54 +01:00
if ( lattice_index > 1 | | lattice_index < 0 ) {
music : : elog < < " For Bravais lattice type, lattice index must be 0 or 1. " < < std : : endl ;
abort ( ) ;
}
2020-05-04 00:49:11 +02:00
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
2022-02-04 00:33:54 +01:00
if ( ishift = = 0 & & lattice_index = = 1 )
2020-05-04 00:49:11 +02:00
{
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_vel64 ( ipcount + + , idim , field . relem ( i , j , k ) ) ;
}
else
{
particles_ . set_vel32 ( ipcount + + , idim , field . relem ( i , j , k ) ) ;
}
}
2020-05-03 04:20:12 +02:00
}
2019-11-01 12:03:02 +01:00
}
2019-11-01 04:47:02 +01:00
}
}
2022-02-04 00:33:54 +01:00
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 ) ) ;
}
}
}
}
}
2020-05-04 00:49:11 +02:00
else
{
2020-05-04 01:13:07 +02:00
glass_ptr_ - > update_ghosts ( field ) ;
2020-05-04 02:09:42 +02:00
# pragma omp parallel for
2020-05-04 00:49:11 +02:00
for ( size_t i = 0 ; i < glass_ptr_ - > size ( ) ; + + i )
{
auto pos = glass_ptr_ - > glass_posr [ i ] ;
real_t vel = glass_ptr_ - > get_at ( pos ) ;
if ( b64reals )
{
particles_ . set_vel64 ( i , idim , vel ) ;
}
else
{
particles_ . set_vel32 ( i , idim , vel ) ;
}
}
}
2019-11-01 04:47:02 +01:00
}
2020-05-03 04:20:12 +02:00
2020-05-04 00:49:11 +02:00
const particle : : container & get_particles ( ) const noexcept {
return particles_ ;
2020-05-03 04:20:12 +02:00
}
2019-11-01 04:47:02 +01:00
2020-05-04 00:49:11 +02:00
} ; // struct lattice
2019-11-01 04:47:02 +01:00
2020-05-04 00:49:11 +02:00
} // namespace particle