mirror of
https://github.com/cosmo-sims/monofonIC.git
synced 2024-09-19 17:03:45 +02:00
mpi bugfixes, refactoring of particle creation
This commit is contained in:
Oliver Hahn
parent
226a9303db
commit
9ebe27cb68
3 changed files with 311 additions and 291 deletions
|
|
@ -16,19 +16,8 @@ struct grid_interpolate
|
|||
static constexpr bool is_distributed_trait = grid_t::is_distributed_trait;
|
||||
static constexpr int interpolation_order = interp_order;
|
||||
|
||||
|
||||
#if defined(USE_MPI)
|
||||
const MPI_Datatype MPI_data_t_type =
|
||||
(typeid(data_t) == typeid(float)) ? MPI_FLOAT
|
||||
: (typeid(data_t) == typeid(double)) ? MPI_DOUBLE
|
||||
: (typeid(data_t) == typeid(long double)) ? MPI_LONG_DOUBLE
|
||||
: (typeid(data_t) == typeid(std::complex<float>)) ? MPI_C_FLOAT_COMPLEX
|
||||
: (typeid(data_t) == typeid(std::complex<double>)) ? MPI_C_DOUBLE_COMPLEX
|
||||
: (typeid(data_t) == typeid(std::complex<long double>)) ? MPI_C_LONG_DOUBLE_COMPLEX
|
||||
: MPI_INT;
|
||||
#endif
|
||||
|
||||
std::vector<data_t> boundary_;
|
||||
std::vector<int> local0starts_;
|
||||
const grid_t &gridref;
|
||||
size_t nx_, ny_, nz_;
|
||||
|
||||
|
|
@ -40,6 +29,13 @@ struct grid_interpolate
|
|||
if (is_distributed_trait)
|
||||
{
|
||||
#if defined(USE_MPI)
|
||||
|
||||
int local_0_start = int(gridref.local_0_start_);
|
||||
local0starts_.assign(MPI::get_size(), 0);
|
||||
|
||||
MPI_Allgather(&local_0_start, 1, MPI_INT, &local0starts_[0], 1, MPI_INT, MPI_COMM_WORLD);
|
||||
|
||||
//... exchange boundary
|
||||
size_t nx = interpolation_order + 1;
|
||||
size_t ny = g.n_[1];
|
||||
size_t nz = g.n_[2];
|
||||
|
|
@ -96,6 +92,11 @@ struct grid_interpolate
|
|||
size_t iz1 = (iz + 1) % nz_;
|
||||
|
||||
data_t val{0.0};
|
||||
|
||||
if( get_task(pos) != MPI::get_rank() ){
|
||||
std::cout << "task : " << MPI::get_rank() << " p@(" << pos[0] << ", " << pos[1] << ", " << pos[2] << ") belongs to task " << get_task(pos) << std::endl;
|
||||
abort();
|
||||
}
|
||||
|
||||
if( is_distributed_trait ){
|
||||
ptrdiff_t localix = ix-gridref.local_0_start_;
|
||||
|
|
@ -135,10 +136,10 @@ struct grid_interpolate
|
|||
// {
|
||||
// }
|
||||
|
||||
int get_task(const vec3 &x, const std::vector<int> &local0starts) const noexcept
|
||||
int get_task(const vec3 &x) const noexcept
|
||||
{
|
||||
const auto it = std::upper_bound(local0starts.begin(), local0starts.end(), int(x[0]));
|
||||
return std::distance(local0starts.begin(), it)-1;
|
||||
const auto it = std::upper_bound(local0starts_.begin(), local0starts_.end(), int(x[0]));
|
||||
return std::distance(local0starts_.begin(), it)-1;
|
||||
}
|
||||
|
||||
void domain_decompose_pos(std::vector<vec3> &pos) const noexcept
|
||||
|
|
@ -146,16 +147,12 @@ struct grid_interpolate
|
|||
if (is_distributed_trait)
|
||||
{
|
||||
#if defined(USE_MPI)
|
||||
int local_0_start = int(gridref.local_0_start_);
|
||||
std::vector<int> local0starts(MPI::get_size(), 0);
|
||||
MPI_Alltoall(&local_0_start, 1, MPI_INT, &local0starts[0], 1, MPI_INT, MPI_COMM_WORLD);
|
||||
|
||||
std::sort(pos.begin(), pos.end(), [&](auto x1, auto x2) { return get_task(x1,local0starts) < get_task(x2,local0starts); });
|
||||
std::sort(pos.begin(), pos.end(), [&](auto x1, auto x2) { return get_task(x1) < get_task(x2); });
|
||||
std::vector<int> sendcounts(MPI::get_size(), 0), sendoffsets(MPI::get_size(), 0);
|
||||
std::vector<int> recvcounts(MPI::get_size(), 0), recvoffsets(MPI::get_size(), 0);
|
||||
for (auto x : pos)
|
||||
{
|
||||
sendcounts[get_task(x,local0starts)] += 3;
|
||||
sendcounts[get_task(x)] += 3;
|
||||
}
|
||||
|
||||
MPI_Alltoall(&sendcounts[0], 1, MPI_INT, &recvcounts[0], 1, MPI_INT, MPI_COMM_WORLD);
|
||||
|
|
@ -169,13 +166,12 @@ struct grid_interpolate
|
|||
tot_send += sendcounts[i];
|
||||
}
|
||||
|
||||
std::vector<vec3> recvbuf;
|
||||
recvbuf.assign(tot_receive,{0.,0.,0.});
|
||||
std::vector<vec3> recvbuf(tot_receive/3,{0.,0.,0.});
|
||||
|
||||
MPI_Alltoallv(&pos[0], &sendcounts[0], &sendoffsets[0], MPI_data_t_type,
|
||||
&recvbuf[0], &recvcounts[0], &recvoffsets[0], MPI_data_t_type, MPI_COMM_WORLD);
|
||||
MPI_Alltoallv(&pos[0], &sendcounts[0], &sendoffsets[0], MPI::get_datatype<real_t>(),
|
||||
&recvbuf[0], &recvcounts[0], &recvoffsets[0], MPI::get_datatype<real_t>(), MPI_COMM_WORLD);
|
||||
|
||||
std::swap( pos, recvbuf );
|
||||
pos.swap( recvbuf );
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
|
@ -193,26 +189,4 @@ struct grid_interpolate
|
|||
return std::exp(ccomplex_t(0.0, shift)) / del;
|
||||
}
|
||||
|
||||
void get_at(std::vector<vec3> &pos, std::vector<data_t> &val) const
|
||||
{
|
||||
|
||||
val.assign( pos.size(), data_t{0.0} );
|
||||
|
||||
for( size_t i=0; i<pos.size(); ++i ){
|
||||
const vec3& x = pos[i];
|
||||
|
||||
switch (interpolation_order)
|
||||
{
|
||||
case 0:
|
||||
val[i] = get_ngp_at(x);
|
||||
break;
|
||||
case 1:
|
||||
val[i] = get_cic_at(x);
|
||||
break;
|
||||
// case 2:
|
||||
// val[i] = get_tsc_at(x);
|
||||
// break;
|
||||
};
|
||||
}
|
||||
}
|
||||
};
|
||||
|
|
@ -14,264 +14,298 @@
|
|||
#include "HDF_IO.hh"
|
||||
#endif
|
||||
|
||||
namespace particle {
|
||||
|
||||
enum lattice{
|
||||
lattice_glass = -1,
|
||||
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< std::vector<vec3_t<real_t>> > lattice_shifts =
|
||||
{
|
||||
// first shift must always be zero! (otherwise set_positions and set_velocities break)
|
||||
/* SC : */ {{0.0,0.0,0.0}},
|
||||
/* BCC: */ {{0.0,0.0,0.0},{0.5,0.5,0.5}},
|
||||
/* FCC: */ {{0.0,0.0,0.0},{0.0,0.5,0.5},{0.5,0.0,0.5},{0.5,0.5,0.0}},
|
||||
/* RSC: */ {{0.0,0.0,0.0},{0.0,0.0,0.5},{0.0,0.5,0.0},{0.0,0.5,0.5},{0.5,0.0,0.0},{0.5,0.0,0.5},{0.5,0.5,0.0},{0.5,0.5,0.5}},
|
||||
};
|
||||
|
||||
const std::vector<vec3_t<real_t>> second_lattice_shift =
|
||||
namespace particle
|
||||
{
|
||||
/* SC : */ {0.5, 0.5, 0.5}, // this corresponds to CsCl lattice
|
||||
/* BCC: */ {0.5, 0.5, 0.0}, // is there a diatomic lattice with BCC base?!?
|
||||
/* FCC: */ {0.5, 0.5, 0.5}, // this corresponds to NaCl lattice
|
||||
// /* FCC: */ {0.25, 0.25, 0.25}, // this corresponds to Zincblende/GaAs lattice
|
||||
/* RSC: */ {0.25, 0.25, 0.25},
|
||||
};
|
||||
using vec3 = std::array<real_t,3>;
|
||||
|
||||
template<typename field_t>
|
||||
void initialize_lattice( container& particles, lattice lattice_type, const bool b64reals, const bool b64ids, const size_t IDoffset, const field_t& field, config_file& cf ){
|
||||
if( lattice_type != lattice_glass )
|
||||
enum lattice
|
||||
{
|
||||
// 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
|
||||
particles.allocate( overload * num_p_in_load, b64reals, b64ids );
|
||||
// set particle IDs to the Lagrangian coordinate (1D encoded) with additionally the field shift encoded as well
|
||||
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,++ipcount){
|
||||
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 );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
lattice_glass = -1,
|
||||
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<std::vector<vec3_t<real_t>>> lattice_shifts =
|
||||
{
|
||||
// first shift must always be zero! (otherwise set_positions and set_velocities break)
|
||||
/* SC : */ {{0.0, 0.0, 0.0}},
|
||||
/* BCC: */ {{0.0, 0.0, 0.0}, {0.5, 0.5, 0.5}},
|
||||
/* FCC: */ {{0.0, 0.0, 0.0}, {0.0, 0.5, 0.5}, {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}},
|
||||
/* RSC: */ {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.5}, {0.0, 0.5, 0.0}, {0.0, 0.5, 0.5}, {0.5, 0.0, 0.0}, {0.5, 0.0, 0.5}, {0.5, 0.5, 0.0}, {0.5, 0.5, 0.5}},
|
||||
};
|
||||
|
||||
const std::vector<vec3_t<real_t>> second_lattice_shift =
|
||||
{
|
||||
/* SC : */ {0.5, 0.5, 0.5}, // this corresponds to CsCl lattice
|
||||
/* BCC: */ {0.5, 0.5, 0.0}, // is there a diatomic lattice with BCC base?!?
|
||||
/* FCC: */ {0.5, 0.5, 0.5}, // this corresponds to NaCl lattice
|
||||
// /* FCC: */ {0.25, 0.25, 0.25}, // this corresponds to Zincblende/GaAs lattice
|
||||
/* RSC: */ {0.25, 0.25, 0.25},
|
||||
};
|
||||
|
||||
template <typename field_t>
|
||||
class lattice_generator
|
||||
{
|
||||
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");
|
||||
|
||||
#if defined(USE_HDF5)
|
||||
std::string glass_fname = cf.get_value<std::string>("setup","GlassFileName");
|
||||
|
||||
std::vector<int> glass_dims;
|
||||
HDFGetDatasetExtent( glass_fname, "/PartType1/Coordinates", glass_dims );
|
||||
music::ilog << "Glass file contains " << glass_dims[0] << " particles." << std::endl;
|
||||
|
||||
size_t ntiles = cf.get_value<size_t>("setup","GlassTiles");
|
||||
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;
|
||||
#if defined(USE_MPI)
|
||||
size_t num_p = glass_dims[0] * ntiles*ntiles*ntiles / MPI::get_size();
|
||||
size_t off_p = MPI::get_rank() * num_p;
|
||||
num_p = np_in_file * ntiles * ntiles * ntiles / MPI::get_size();
|
||||
off_p = MPI::get_rank() * num_p;
|
||||
#else
|
||||
size_t num_p = glass_dims[0] * ntiles*ntiles*ntiles;
|
||||
size_t off_p = 0;
|
||||
num_p = np_in_file * ntiles * ntiles * ntiles;
|
||||
off_p = 0;
|
||||
#endif
|
||||
|
||||
particles.allocate( num_p, b64reals, b64ids );
|
||||
glass_posr.assign(num_p, {0.0, 0.0, 0.0});
|
||||
|
||||
for( size_t i=0; i<num_p; ++i ){
|
||||
if( b64ids ){
|
||||
particles.set_id64( i, IDoffset + i + off_p );
|
||||
}else{
|
||||
particles.set_id32( i, IDoffset + i + off_p );
|
||||
}
|
||||
}
|
||||
#else
|
||||
throw std::runtime_error("Class lattice requires HDF5 support. Enable and recompile.");
|
||||
std::array<real_t, 3> ng({real_t(field.n_[0]), real_t(field.n_[1]), real_t(field.n_[2])});
|
||||
|
||||
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];
|
||||
}
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
// invalidates field, phase shifted to unspecified position after return
|
||||
template<typename field_t>
|
||||
void set_positions( container& particles, const lattice lattice_type, bool is_second_lattice, int idim, real_t lunit, 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();
|
||||
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] );
|
||||
data_t get_at( const vec3& x ) const noexcept
|
||||
{
|
||||
return interp_.get_cic_at( x );
|
||||
}
|
||||
// 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){
|
||||
auto pos = field.template get_unit_r_shifted<real_t>(i,j,k,lattice_shifts[lattice_type][ishift]
|
||||
+ (is_second_lattice? second_lattice_shift[lattice_type] : vec3_t<real_t>{0.,0.,0.}) );
|
||||
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) );
|
||||
|
||||
size_t size() const noexcept
|
||||
{
|
||||
return num_p;
|
||||
}
|
||||
|
||||
size_t offset() const noexcept
|
||||
{
|
||||
return off_p;
|
||||
}
|
||||
};
|
||||
|
||||
std::unique_ptr<glass> glass_ptr_;
|
||||
|
||||
private:
|
||||
particle::container particles_;
|
||||
|
||||
public:
|
||||
lattice_generator(lattice lattice_type, const bool b64reals, const bool b64ids, const size_t IDoffset, const field_t &field, config_file &cf)
|
||||
{
|
||||
if (lattice_type != lattice_glass)
|
||||
{
|
||||
// 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
|
||||
particles_.allocate(overload * num_p_in_load, b64reals, b64ids);
|
||||
// set particle IDs to the Lagrangian coordinate (1D encoded) with additionally the field shift encoded as well
|
||||
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, ++ipcount)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}else{
|
||||
#if defined(USE_HDF5)
|
||||
std::string glass_fname = cf.get_value<std::string>("setup","GlassFileName");
|
||||
size_t ntiles = cf.get_value<size_t>("setup","GlassTiles");
|
||||
|
||||
real_t lglassbox = 1.0;
|
||||
HDFReadGroupAttribute( glass_fname, "Header", "BoxSize", lglassbox );
|
||||
else
|
||||
{
|
||||
glass_ptr_ = std::make_unique<glass>( cf, field );
|
||||
particles_.allocate(glass_ptr_->size(), b64reals, b64ids);
|
||||
|
||||
std::vector<real_t> glass_pos;
|
||||
HDFReadDataset( glass_fname, "/PartType1/Coordinates", glass_pos );
|
||||
size_t np_in_file = glass_pos.size()/3;
|
||||
#if defined(USE_MPI)
|
||||
size_t num_p = np_in_file * ntiles*ntiles*ntiles / MPI::get_size();
|
||||
size_t off_p = MPI::get_rank() * num_p;
|
||||
#else
|
||||
size_t num_p = np_in_file * ntiles*ntiles*ntiles;
|
||||
size_t off_p = 0;
|
||||
#endif
|
||||
|
||||
std::vector< std::array<real_t,3> > glass_posr(num_p,{0.0,0.0,0.0});
|
||||
|
||||
std::array<real_t,3> ng({real_t(field.n_[0]),real_t(field.n_[1]),real_t(field.n_[2])});
|
||||
|
||||
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]);
|
||||
}
|
||||
|
||||
grid_interpolate<1,field_t> interp( field );
|
||||
|
||||
interp.domain_decompose_pos( glass_posr );
|
||||
|
||||
for( size_t i=0; i<num_p; ++i ){
|
||||
auto pos = glass_posr[i];
|
||||
real_t disp = interp.get_cic_at( pos );
|
||||
if( b64reals ){
|
||||
particles.set_pos64( i, idim, pos[idim]/ng[idim]*lunit + disp );
|
||||
}else{
|
||||
particles.set_pos32( i, idim, pos[idim]/ng[idim]*lunit + disp );
|
||||
}
|
||||
}
|
||||
|
||||
#else
|
||||
throw std::runtime_error("Class lattice requires HDF5 support. Enable and recompile.");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
template <typename field_t>
|
||||
void set_velocities(container &particles, 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();
|
||||
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_vel64( ipcount++, idim, field.relem(i,j,k) );
|
||||
}else{
|
||||
particles.set_vel32( ipcount++, idim, field.relem(i,j,k) );
|
||||
}
|
||||
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());
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}else{
|
||||
#if defined(USE_HDF5)
|
||||
std::string glass_fname = cf.get_value<std::string>("setup","GlassFileName");
|
||||
size_t ntiles = cf.get_value<size_t>("setup","GlassTiles");
|
||||
|
||||
real_t lglassbox = 1.0;
|
||||
HDFReadGroupAttribute( glass_fname, "Header", "BoxSize", lglassbox );
|
||||
|
||||
std::vector<real_t> glass_pos;
|
||||
HDFReadDataset( glass_fname, "/PartType1/Coordinates", glass_pos );
|
||||
size_t np_in_file = glass_pos.size()/3;
|
||||
#if defined(USE_MPI)
|
||||
size_t num_p = np_in_file * ntiles*ntiles*ntiles / MPI::get_size();
|
||||
size_t off_p = MPI::get_rank() * num_p;
|
||||
#else
|
||||
size_t num_p = np_in_file * ntiles*ntiles*ntiles;
|
||||
size_t off_p = 0;
|
||||
#endif
|
||||
// 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)
|
||||
{
|
||||
// works only for Bravais types
|
||||
if (lattice_type >= 0)
|
||||
{
|
||||
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)
|
||||
{
|
||||
field.shift_field(second_lattice_shift[lattice_type]);
|
||||
}
|
||||
|
||||
|
||||
std::vector< std::array<real_t,3> > glass_posr(num_p,{0.0,0.0,0.0});
|
||||
|
||||
std::array<real_t,3> ng({real_t(field.n_[0]),real_t(field.n_[1]),real_t(field.n_[2])});
|
||||
|
||||
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]);
|
||||
}
|
||||
|
||||
grid_interpolate<1,field_t> interp( field );
|
||||
|
||||
interp.domain_decompose_pos( glass_posr );
|
||||
|
||||
for( size_t i=0; i<num_p; ++i ){
|
||||
auto pos = glass_posr[i];
|
||||
real_t vel = interp.get_cic_at( pos );
|
||||
if( b64reals ){
|
||||
particles.set_vel64( i, idim, vel );
|
||||
}else{
|
||||
particles.set_vel32( i, idim, vel );
|
||||
// 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)
|
||||
{
|
||||
auto pos = field.template get_unit_r_shifted<real_t>(i, j, k, lattice_shifts[lattice_type][ishift] + (is_second_lattice ? second_lattice_shift[lattice_type] : vec3_t<real_t>{0., 0., 0.}));
|
||||
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
|
||||
{
|
||||
for (size_t i = 0; i < this->glass_ptr_->size(); ++i)
|
||||
{
|
||||
auto pos = this->glass_ptr_->glass_posr[i];
|
||||
real_t disp = this->glass_ptr_->get_at(pos);
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
#else
|
||||
throw std::runtime_error("Class lattice requires HDF5 support. Enable and recompile.");
|
||||
#endif
|
||||
}
|
||||
}
|
||||
|
||||
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();
|
||||
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_vel64(ipcount++, idim, field.relem(i, j, k));
|
||||
}
|
||||
else
|
||||
{
|
||||
particles_.set_vel32(ipcount++, idim, field.relem(i, j, k));
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // end namespace particles
|
||||
const particle::container& get_particles() const noexcept{
|
||||
return particles_;
|
||||
}
|
||||
|
||||
}; // struct lattice
|
||||
|
||||
} // namespace particle
|
||||
|
|
|
|||
|
|
@ -452,6 +452,19 @@ int Run( config_file& the_config )
|
|||
// temporary storage of data
|
||||
Grid_FFT<real_t> tmp({ngrid, ngrid, ngrid}, {boxlen, boxlen, boxlen});
|
||||
|
||||
std::unique_ptr<particle::lattice_generator<Grid_FFT<real_t>>> particle_lattice_generator_ptr;
|
||||
|
||||
// if output plugin wants particles, then we need to store them, along with their IDs
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
// somewhat arbitrarily, start baryon particle IDs from 2**31 if we have 32bit and from 2**56 if we have 64 bits
|
||||
size_t IDoffset = (this_species == cosmo_species::baryon)? ((the_output_plugin->has_64bit_ids())? 1ul<<56 : 1ul<<31): 0 ;
|
||||
|
||||
// allocate particle structure and generate particle IDs
|
||||
particle_lattice_generator_ptr =
|
||||
std::make_unique<particle::lattice_generator<Grid_FFT<real_t>>>( lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), IDoffset, tmp, 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 )
|
||||
|
|
@ -542,22 +555,21 @@ int Run( config_file& the_config )
|
|||
//===================================================================================
|
||||
// we store displacements and velocities here if we compute them
|
||||
//===================================================================================
|
||||
particle::container particles;
|
||||
|
||||
|
||||
bool shifted_lattice = (this_species == cosmo_species::baryon &&
|
||||
the_output_plugin->write_species_as(this_species) == output_type::particles) ? true : false;
|
||||
|
||||
// somewhat arbitrarily, start baryon particle IDs from 2**31 if we have 32bit and from 2**56 if we have 64 bits
|
||||
size_t IDoffset = (this_species == cosmo_species::baryon)? ((the_output_plugin->has_64bit_ids())? 1ul<<56 : 1ul<<31): 0 ;
|
||||
|
||||
|
||||
grid_interpolate<1,Grid_FFT<real_t>> interp( tmp );
|
||||
|
||||
// if output plugin wants particles, then we need to store them, along with their IDs
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
// allocate particle structure and generate particle IDs
|
||||
particle::initialize_lattice( particles, lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), IDoffset, tmp, the_config );
|
||||
}
|
||||
// if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
// {
|
||||
// // allocate particle structure and generate particle IDs
|
||||
// particle::initialize_lattice( particles, lattice_type, the_output_plugin->has_64bit_reals(), the_output_plugin->has_64bit_ids(), IDoffset, tmp, the_config );
|
||||
// }
|
||||
|
||||
// write out positions
|
||||
for( int idim=0; idim<3; ++idim ){
|
||||
|
|
@ -604,7 +616,7 @@ int Run( config_file& the_config )
|
|||
// if we write particle data, store particle data in particle structure
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
particle::set_positions( particles, lattice_type, shifted_lattice, idim, lunit, the_output_plugin->has_64bit_reals(), tmp, the_config );
|
||||
particle_lattice_generator_ptr->set_positions( lattice_type, shifted_lattice, idim, lunit, the_output_plugin->has_64bit_reals(), tmp, the_config );
|
||||
}
|
||||
// otherwise write out the grid data directly to the output plugin
|
||||
// else if( the_output_plugin->write_species_as( cosmo_species::dm ) == output_type::field_lagrangian )
|
||||
|
|
@ -672,7 +684,7 @@ int Run( config_file& the_config )
|
|||
// if we write particle data, store particle data in particle structure
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
particle::set_velocities( particles, lattice_type, shifted_lattice, idim, the_output_plugin->has_64bit_reals(), tmp, the_config );
|
||||
particle_lattice_generator_ptr->set_velocities( lattice_type, shifted_lattice, idim, the_output_plugin->has_64bit_reals(), tmp, the_config );
|
||||
}
|
||||
// otherwise write out the grid data directly to the output plugin
|
||||
else if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
|
||||
|
|
@ -684,7 +696,7 @@ int Run( config_file& the_config )
|
|||
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
|
||||
{
|
||||
the_output_plugin->write_particle_data( particles, this_species, Omega[this_species] );
|
||||
the_output_plugin->write_particle_data( particle_lattice_generator_ptr->get_particles(), this_species, Omega[this_species] );
|
||||
}
|
||||
|
||||
if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
|
||||
|
|
|
|||
Loading…
Reference in a new issue