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

added fcc lattice to possible initial particle loads

Browse source
This commit is contained in:
Oliver Hahn 2019-10-17 17:39:26 +02:00
parent 35910d724b
commit e9b11b1bae
4 changed files with 171 additions and 39 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

16
example.conf
View file

@ -1,12 +1,18 @@
[setup]
GridRes = 64
BoxLength = 100
# number of grid cells per linear dimension for calculations = particles for sc initial load
GridRes = 128
# length of the box in Mpc/h
BoxLength = 250
# starting redshift
zstart = 49.0
LPTorder = 2
SymplecticPT = no
# order of the LPT to be used (1,2 or 3)
LPTorder = 3
# also do baryon ICs?
DoBaryons = no
# do mode fixing à la Angulo&Pontzen
DoFixing = no
BCClattice = no
# particle load, can be 'sc' (1x), 'bcc' (2x) or 'fcc' (4x) (increases number of particles by factor!)
ParticleLoad = sc
[testing]
# enables diagnostic output

37
include/grid_fft.hh
View file

@ -188,6 +188,18 @@ public:
return rr;
}
template <typename ft>
vec3<ft> get_unit_r_shifted(const size_t i, const size_t j, const size_t k, double sx, double sy, double sz) const
{
vec3<ft> rr;
rr[0] = (real_t(i + local_0_start_) + sx) / real_t(n_[0]);
rr[1] = (real_t(j) + sy) / real_t(n_[1]);
rr[2] = (real_t(k) + sz) / real_t(n_[2]);
return rr;
}
void cell_pos(int ilevel, size_t i, size_t j, size_t k, double *x) const
{
x[0] = double(i + local_0_start_) / size(0);
@ -643,16 +655,35 @@ public:
void Write_PDF(std::string ofname, int nbins = 1000, double scale = 1.0, double rhomin = 1e-3, double rhomax = 1e3);
void stagger_field(void)
// void stagger_field(void)
// {
// FourierTransformForward();
// apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
// real_t shift = k[0] * get_dx()[0] + k[1] * get_dx()[1] + k[2] * get_dx()[2];
// return x * std::exp(ccomplex_t(0.0, 0.5 * shift));
// });
// FourierTransformBackward();
// }
void shift_field( double sx, double sy, double sz )
{
FourierTransformForward();
apply_function_k_dep([&](auto x, auto k) -> ccomplex_t {
real_t shift = k[0] * get_dx()[0] + k[1] * get_dx()[1] + k[2] * get_dx()[2];
return x * std::exp(ccomplex_t(0.0, 0.5 * shift));
#ifdef WITH_MPI
real_t shift = sy * k[0] * get_dx()[0] + sx * k[1] * get_dx()[1] + sz * k[2] * get_dx()[2];
#else
real_t shift = sx * k[0] * get_dx()[0] + sy * k[1] * get_dx()[1] + sz * k[2] * get_dx()[2];
#endif
return x * std::exp(ccomplex_t(0.0, shift));
});
FourierTransformBackward();
}
void stagger_field(void)
{
this->shift_field( 0.5, 0.5, 0.5 );
}
void zero_DC_mode(void)
{
if (space_ == kspace_id)

15
include/ic_generator.hh
View file

@ -9,11 +9,16 @@
namespace ic_generator{
int Run( ConfigFile& the_config );
int Initialise( ConfigFile& the_config );
enum particle_lattice{
lattice_sc, lattice_bcc, lattice_fcc
};
extern std::unique_ptr<RNG_plugin> the_random_number_generator;
extern std::unique_ptr<output_plugin> the_output_plugin;
extern std::unique_ptr<CosmologyCalculator> the_cosmo_calc;
int Run( ConfigFile& the_config );
int Initialise( ConfigFile& the_config );
extern std::unique_ptr<RNG_plugin> the_random_number_generator;
extern std::unique_ptr<output_plugin> the_output_plugin;
extern std::unique_ptr<CosmologyCalculator> the_cosmo_calc;
}

142
src/ic_generator.cc
View file

@ -55,8 +55,9 @@ int Run( ConfigFile& the_config )
int LPTorder = the_config.GetValueSafe<double>("setup","LPTorder",100);
//--------------------------------------------------------------------------------------------------------
//! initialice particles on a bcc lattice instead of a standard sc lattice (doubles number of particles)
const bool initial_bcc_lattice = the_config.GetValueSafe<bool>("setup","BCClattice",false);
//! 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);
//--------------------------------------------------------------------------------------------------------
//! apply fixing of the complex mode amplitude following Angulo & Pontzen (2016) [https://arxiv.org/abs/1603.05253]
@ -439,29 +440,50 @@ int Run( ConfigFile& the_config )
// 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 )
{
// if particles occupy a bcc lattice, then there are 2 x N^3 of them...
// generate particle IDs
if( !initial_bcc_lattice ){
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) );
// 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) );
}
}
}
}
}else{
// 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_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;
}
}
@ -499,13 +521,48 @@ int Run( ConfigFile& the_config )
}
}
if( initial_bcc_lattice ){
tmp.stagger_field();
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_staggered<float>(i,j,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) );
}
}
@ -563,8 +620,8 @@ int Run( ConfigFile& the_config )
}
}
if( initial_bcc_lattice ){
tmp.stagger_field();
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){
@ -574,6 +631,39 @@ int Run( ConfigFile& the_config )
}
}
}
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
else if( the_output_plugin->write_species_as( this_species ) == output_type::field_lagrangian )
{