mirror/monofonIC
1
0
Fork 0
mirror of https://github.com/cosmo-sims/monofonIC.git synced 2024-09-19 17:03:45 +02:00
Code Releases Activity

moved particle lattice creation to separate function in separate header

Browse source
This commit is contained in:
Oliver Hahn 2019-10-19 12:56:43 +02:00
parent e9b11b1bae
commit a58fbc9778
6 changed files with 171 additions and 172 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

4
include/ic_generator.hh
View file

@ -9,10 +9,6 @@
namespace ic_generator{
enum particle_lattice{
lattice_sc, lattice_bcc, lattice_fcc
};
int Run( ConfigFile& the_config );
int Initialise( ConfigFile& the_config );

2
include/output_plugin.hh
View file

@ -44,7 +44,7 @@ public:
virtual ~output_plugin(){}
//! routine to write particle data for a species
virtual void write_particle_data(const particle_container &pc, const cosmo_species &s ) {};
virtual void write_particle_data(const particle::container &pc, const cosmo_species &s ) {};
//! routine to write gridded fluid component data for a species
virtual void write_grid_data(const Grid_FFT<real_t> &g, const cosmo_species &s, const fluid_component &c ) {};

12
include/particle_container.hh
View file

@ -8,17 +8,19 @@
#include <vector>
#include <general.hh>
class particle_container
namespace particle{
class container
{
public:
std::vector<float> positions_, velocities_;
std::vector<int> ids_;
particle_container()
container()
{
}
particle_container(const particle_container &) = delete;
container(const container &) = delete;
const void* get_pos_ptr() const{
return reinterpret_cast<const void*>( &positions_[0] );
@ -82,9 +84,6 @@ public:
std::vector<size_t> nump_p_task(mpi_size,0), off_p_task;
size_t num_p_this_task = this->get_local_num_particles();
// MPI_Allgather(const void *sendbuf, int sendcount, MPI_Datatype sendtype,
// void *recvbuf, int recvcount, MPI_Datatype recvtype,
// MPI_Comm comm)
MPI_Allgather( reinterpret_cast<const void*>(&num_p_this_task), 1, MPI_UNSIGNED_LONG_LONG,
reinterpret_cast<void*>(&nump_p_task[0]), 1, MPI_UNSIGNED_LONG_LONG, MPI_COMM_WORLD );
@ -106,3 +105,4 @@ public:
}
};
}

150
include/particle_generator.hh Normal file
View file

@ -0,0 +1,150 @@
#pragma once
namespace particle {
enum lattice{
lattice_sc=0, lattice_bcc=1, lattice_fcc=2
};
template<typename field_t>
void initialize_lattice( container& particles, lattice lattice_type, const field_t& field ){
const size_t num_p_in_load = field.local_size();
const size_t overload = 1<<lattice_type; // 1 for sc, 2 for bcc, 4 for fcc
particles.allocate( overload * num_p_in_load );
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 ){
particles.set_id( ipcount+iload*num_p_in_load, overload*field.get_cell_idx_1d(i,j,k)+iload );
}
}
}
}
}
// invalidates field, phase shifted to unspecified position after return
template<typename field_t>
void set_positions( container& particles, lattice lattice_type, int idim, real_t lunit, field_t& field )
{
const size_t num_p_in_load = field.local_size();
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){
auto pos = field.template get_unit_r<real_t>(i,j,k);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) );
}
}
}
if( lattice_type == particle::lattice_bcc ){
field.shift_field( 0.5, 0.5, 0.5 );
auto ipcount0 = 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,0.5,0.5,0.5);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) );
}
}
}
}
else if( lattice_type == particle::lattice_fcc ){
// 0.5 0.5 0.0
field.shift_field( 0.5, 0.5, 0.0 );
auto ipcount0 = 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,0.5,0.5,0.0);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) );
}
}
}
// 0.0 0.5 0.5
field.shift_field( -0.5, 0.0, 0.5 );
ipcount0 = 2*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,0.0,0.5,0.5);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) );
}
}
}
// 0.5 0.0 0.5
field.shift_field( 0.5, -0.5, 0.0 );
ipcount0 = 3*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,0.5,0.0,0.5);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + field.relem(i,j,k) );
}
}
}
}
}
template<typename field_t>
void set_velocities( container& particles, lattice lattice_type, int idim, field_t& field )
{
const size_t num_p_in_load = field.local_size();
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){
particles.set_vel( ipcount++, idim, field.relem(i,j,k) );
}
}
}
if( lattice_type == particle::lattice_bcc ){
field.shift_field( 0.5, 0.5, 0.5 );
auto ipcount0 = 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){
particles.set_vel( ipcount++, idim, field.relem(i,j,k) );
}
}
}
}
else if( lattice_type == particle::lattice_fcc ){
// 0.5 0.5 0.0
field.shift_field( 0.5, 0.5, 0.0 );
auto ipcount0 = 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){
particles.set_vel( ipcount++, idim, field.relem(i,j,k) );
}
}
}
// 0.0 0.5 0.5
field.shift_field( -0.5, 0.0, 0.5 );
ipcount0 = 2*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){
particles.set_vel( ipcount++, idim, field.relem(i,j,k) );
}
}
}
// 0.5 0.0 0.5
field.shift_field( 0.5, -0.5, 0.0 );
ipcount0 = 3*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){
particles.set_vel( ipcount++, idim, field.relem(i,j,k) );
}
}
}
}
}
} // end namespace particles

173
src/ic_generator.cc
View file

@ -5,8 +5,8 @@
#include <convolution.hh>
#include <testing.hh>
#include <ic_generator.hh>
#include <particle_generator.hh>
#include <unistd.h> // for unlink
@ -57,7 +57,8 @@ int Run( ConfigFile& the_config )
//--------------------------------------------------------------------------------------------------------
//! initialice particles on a bcc or fcc lattice instead of a standard sc lattice (doubles and quadruples the number of particles)
std::string lattice_str = the_config.GetValueSafe<std::string>("setup","ParticleLoad","sc");
const particle_lattice lattice_type = (lattice_str=="bcc")? lattice_bcc : ((lattice_str=="fcc")? lattice_fcc : lattice_sc);
const particle::lattice lattice_type = (lattice_str=="bcc")? particle::lattice_bcc
: ((lattice_str=="fcc")? particle::lattice_fcc : particle::lattice_sc);
//--------------------------------------------------------------------------------------------------------
//! apply fixing of the complex mode amplitude following Angulo & Pontzen (2016) [https://arxiv.org/abs/1603.05253]
@ -431,60 +432,16 @@ int Run( ConfigFile& the_config )
if( the_output_plugin->write_species_as( this_species ) == output_type::particles
|| the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
{
//======================================================================
//===================================================================================
// we store displacements and velocities here if we compute them
//======================================================================
size_t num_p_in_load = phi.local_size();
particle_container particles;
//===================================================================================
particle::container particles;
// 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
switch( lattice_type ){
default:
case lattice_sc:
// if particles occupy a sc lattice, then there are 2 x N^3 of them...
particles.allocate( num_p_in_load );
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_id( ipcount++, tmp.get_cell_idx_1d(i,j,k) );
}
}
}
break;
case lattice_bcc:
// if particles occupy a bcc lattice, then there are 2 x N^3 of them...
particles.allocate( 2*num_p_in_load );
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_id( ipcount, 2*tmp.get_cell_idx_1d(i,j,k) );
particles.set_id( ipcount+num_p_in_load, 2*tmp.get_cell_idx_1d(i,j,k)+1 );
++ipcount;
}
}
}
break;
case lattice_fcc:
// if particles occupy a fcc lattice, then there are 4 x N^3 of them...
particles.allocate( 4*num_p_in_load );
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_id( ipcount+0*num_p_in_load, 4*tmp.get_cell_idx_1d(i,j,k) );
particles.set_id( ipcount+1*num_p_in_load, 4*tmp.get_cell_idx_1d(i,j,k)+1 );
particles.set_id( ipcount+2*num_p_in_load, 4*tmp.get_cell_idx_1d(i,j,k)+2 );
particles.set_id( ipcount+3*num_p_in_load, 4*tmp.get_cell_idx_1d(i,j,k)+3 );
++ipcount;
}
}
}
break;
}
particle::initialize_lattice( particles, lattice_type, tmp );
}
// write out positions
@ -512,62 +469,7 @@ int Run( ConfigFile& 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 )
{
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
auto pos = tmp.get_unit_r<float>(i,j,k);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
}
}
}
if( lattice_type == lattice_bcc ){
tmp.shift_field( 0.5, 0.5, 0.5 );
auto ipcount0 = num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
auto pos = tmp.get_unit_r_shifted<float>(i,j,k,0.5,0.5,0.5);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
}
}
}
}
else if( lattice_type == lattice_fcc ){
// 0.5 0.5 0.0
tmp.shift_field( 0.5, 0.5, 0.0 );
auto ipcount0 = num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
auto pos = tmp.get_unit_r_shifted<float>(i,j,k,0.5,0.5,0.0);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
}
}
}
// 0.0 0.5 0.5
tmp.shift_field( -0.5, 0.0, 0.5 );
ipcount0 = 2*num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
auto pos = tmp.get_unit_r_shifted<float>(i,j,k,0.0,0.5,0.5);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
}
}
}
// 0.5 0.0 0.5
tmp.shift_field( 0.5, -0.5, 0.0 );
ipcount0 = 3*num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
auto pos = tmp.get_unit_r_shifted<float>(i,j,k,0.5,0.0,0.5);
particles.set_pos( ipcount++, idim, pos[idim]*lunit + tmp.relem(i,j,k) );
}
}
}
}
particle::set_positions( particles, lattice_type, idim, lunit, tmp );
}
// 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 )
@ -611,60 +513,11 @@ int Run( ConfigFile& the_config )
tmp.FourierTransformBackward();
// if we write particle data, store particle data in particle structure
if( the_output_plugin->write_species_as( this_species ) == output_type::particles ){
for( size_t i=0,ipcount=0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
}
}
}
if( lattice_type == lattice_bcc ){
tmp.shift_field( 0.5, 0.5, 0.5 );
auto ipcount0 = num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
}
}
}
}
else if( lattice_type == lattice_fcc ){
// 0.5 0.5 0.0
tmp.shift_field( 0.5, 0.5, 0.0 );
auto ipcount0 = num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
}
}
}
// 0.0 0.5 0.5
tmp.shift_field( -0.5, 0.0, 0.5 );
ipcount0 = 2*num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
}
}
}
// 0.5 0.0 0.5
tmp.shift_field( 0.5, -0.5, 0.0 );
ipcount0 = 3*num_p_in_load;
for( size_t i=0,ipcount=ipcount0; i<tmp.size(0); ++i ){
for( size_t j=0; j<tmp.size(1); ++j){
for( size_t k=0; k<tmp.size(2); ++k){
particles.set_vel( ipcount++, idim, tmp.relem(i,j,k) );
}
}
}
}
}// otherwise write out the grid data directly to the output plugin
if( the_output_plugin->write_species_as( this_species ) == output_type::particles )
{
particle::set_velocities( particles, lattice_type, idim, tmp );
}
// 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 )
{
fluid_component fc = (idim==0)? fluid_component::vx : ((idim==1)? fluid_component::vy : fluid_component::vz );

2
src/plugins/output_gadget2.cc
View file

@ -55,7 +55,7 @@ public:
real_t velocity_unit() const { return vunit_; }
void write_particle_data(const particle_container &pc, const cosmo_species &s )
void write_particle_data(const particle::container &pc, const cosmo_species &s )
{
// fill the Gadget-2 header
memset(reinterpret_cast<void*>(&this_header_),0,sizeof(header));