diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
index b947e8f28de4ee22b10b3ebe4aa3de500055f708..836a3311917b00ad5c8cc29d9c0dc1e656a91335 100644
--- a/cpp/full_code/test_tracer_set.cpp
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -30,6 +30,7 @@
#include "particles/rhs/tracer_rhs.hpp"
#include "particles/interpolation/particle_set.hpp"
#include "particles/particle_solver.hpp"
+#include "particles/particles_input_random.hpp"
#include "scope_timer.hpp"
@@ -65,6 +66,10 @@ int test_tracer_set<rnumber>::do_work(void)
this->nx, this->ny, this->nz,
this->comm,
FFTW_ESTIMATE);
+
+ *vec_field = 0.0;
+ vec_field->real_space_representation = true;
+
field_tinterpolator<rnumber, FFTW, THREE, NONE> *fti = new field_tinterpolator<rnumber, FFTW, THREE, NONE>();
tracer_rhs<rnumber, FFTW, NONE> *trhs = new tracer_rhs<rnumber, FFTW, NONE>();
@@ -74,12 +79,16 @@ int test_tracer_set<rnumber>::do_work(void)
this->dky,
this->dkz);
+ particles_input_random<long long int, double, 3> bla(
+ this->comm, 1000, 1, pset.getSpatialLowLimitZ(), pset.getSpatialUpLimitZ());
+
fti->set_field(vec_field);
trhs->setVelocity(fti);
+ pset.init(bla);
- particle_solver psolver(pset);
+ //particle_solver psolver(pset);
- psolver.Euler(0.1, *trhs);
+ //psolver.Euler(0.1, *trhs);
// deallocate
delete trhs;
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 5453e1b55507a58bc16e23278c7729d2bf4d2c1b..6f7c82c1ca3a5f66aff90c5f4dd241d4a53b45e0 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -33,7 +33,8 @@
#include "field_layout.hpp"
#include "field.hpp"
#include "particles/particles_distr_mpi.hpp"
-#include "particles/particles_output_hdf5.hpp"
+#include "particles/abstract_particles_output.hpp"
+#include "particles/abstract_particles_input.hpp"
#include "particles/interpolation/particles_field_computer.hpp"
#include "particles/interpolation/particles_generic_interp.hpp"
#include "particles/interpolation/abstract_particle_set.hpp"
@@ -269,6 +270,44 @@ class particle_set: public abstract_particle_set
&this->local_index);
return EXIT_SUCCESS;
}
+
+ int init(abstract_particles_input<partsize_t, particle_rnumber>& particles_input)
+ {
+ TIMEZONE("particles_system::init");
+
+ this->local_state = particles_input.getMyParticles();
+ this->local_index = particles_input.getMyParticlesIndexes();
+ this->local_number_of_particles = particles_input.getLocalNbParticles();
+
+ //particles_utils::partition_extra_z<partsize_t, state_size>(&this->local_state[0], this->local_number_of_particles, partition_interval_size,
+ // this->local_number_of_particles_per_partition.get(), current_offset_particles_for_partition.get(),
+ //[&](const particle_rnumber& z_pos){
+ // const int partition_level = this->pinterpolator.pbc_field_layer(z_pos, IDXC_Z);
+ // assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second);
+ // return partition_level - current_partition_interval.first;
+ //},
+ //[&](const partsize_t idx1, const partsize_t idx2){
+ // std::swap(this->local_index[idx1], this->local_index[idx2]);
+ // for(int idx_rhs = 0 ; idx_rhs < int(my_particles_rhs.size()) ; ++idx_rhs){
+ // for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
+ // std::swap(my_particles_rhs[idx_rhs][idx1*size_particle_rhs + idx_val],
+ // my_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]);
+ // }
+ // }
+ //});
+
+ return EXIT_SUCCESS;
+ }
+
+ particle_rnumber getSpatialLowLimitZ()
+ {
+ return this->pinterpolator.getSpatialLowLimitZ();
+ }
+
+ particle_rnumber getSpatialUpLimitZ()
+ {
+ return this->pinterpolator.getSpatialUpLimitZ();
+ }
};
#endif//PARTICLE_SET_HPP
diff --git a/cpp/particles/interpolation/particles_field_computer.hpp b/cpp/particles/interpolation/particles_field_computer.hpp
index a80a37ceabedfd9e60f2bbd11ea5b11ab284560b..93e056ffef5dfce1b3a1920cf9d61bbfcc9ea1f8 100644
--- a/cpp/particles/interpolation/particles_field_computer.hpp
+++ b/cpp/particles/interpolation/particles_field_computer.hpp
@@ -84,6 +84,29 @@ public:
0);
}
+ const std::array<real_number, 3> getSpatialBoxWidth() const
+ {
+ return this->spatial_box_width;
+ }
+ const std::array<real_number, 3> getSpatialBoxOffset() const
+ {
+ return this->spatial_box_offset;
+ }
+ const std::array<real_number, 3> getBoxStepWidth() const
+ {
+ return this->box_step_width;
+ }
+
+ const real_number getSpatialLowLimitZ() const
+ {
+ return this->current_partition_interval.first*this->box_step_width[IDXC_Z];
+ }
+
+ const real_number getSpatialUpLimitZ() const
+ {
+ return this->current_partition_interval.second*this->box_step_width[IDXC_Z];
+ }
+
real_number get_norm_pos_in_cell(const real_number in_pos, const int idx_pos) const {
const real_number shifted_pos = in_pos - spatial_box_offset[idx_pos];
const real_number nb_box_repeat = floor(shifted_pos/spatial_box_width[idx_pos]);
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 1bd715c129942a3cbf3a78465eeefba02b9f021e..680daff93c55ad1e6fe6951dfc3705b1229d17de 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -35,7 +35,7 @@
#include "abstract_particles_input.hpp"
#include "base.hpp"
#include "alltoall_exchanger.hpp"
-#include "particles_utils.hpp"
+#include "particles/particles_utils.hpp"
#include "scope_timer.hpp"
@@ -56,40 +56,6 @@ class particles_input_hdf5 : public abstract_particles_input<partsize_t, real_nu
std::unique_ptr<partsize_t[]> my_particles_indexes;
std::vector<std::unique_ptr<real_number[]>> my_particles_rhs;
- static std::vector<real_number> BuildLimitsAllProcesses(
- MPI_Comm mpi_comm,
- const real_number my_spatial_low_limit,
- const real_number my_spatial_up_limit){
- int my_rank;
- int nb_processes;
-
- AssertMpi(MPI_Comm_rank(mpi_comm, &my_rank));
- AssertMpi(MPI_Comm_size(mpi_comm, &nb_processes));
-
- std::vector<real_number> spatial_limit_per_proc(nb_processes*2);
-
- real_number intervalToSend[2] = {my_spatial_low_limit, my_spatial_up_limit};
- AssertMpi(
- MPI_Allgather(
- intervalToSend,
- 2,
- particles_utils::GetMpiType(real_number()),
- spatial_limit_per_proc.data(),
- 2,
- particles_utils::GetMpiType(real_number()),
- mpi_comm));
-
- for(int idx_proc = 0; idx_proc < nb_processes-1 ; ++idx_proc){
- assert(spatial_limit_per_proc[idx_proc*2] <= spatial_limit_per_proc[idx_proc*2+1]);
- assert(spatial_limit_per_proc[idx_proc*2+1] == spatial_limit_per_proc[(idx_proc+1)*2]);
- spatial_limit_per_proc[idx_proc+1] = spatial_limit_per_proc[idx_proc*2+1];
- }
- spatial_limit_per_proc[nb_processes] = spatial_limit_per_proc[(nb_processes-1)*2+1];
- spatial_limit_per_proc.resize(nb_processes+1);
-
- return spatial_limit_per_proc;
- }
-
public:
particles_input_hdf5(
const MPI_Comm in_mpi_comm,
@@ -103,7 +69,7 @@ public:
inFilename,
inDatanameState,
inDatanameRhs,
- BuildLimitsAllProcesses(
+ BuildLimitsAllProcesses<real_number>(
in_mpi_comm,
my_spatial_low_limit,
my_spatial_up_limit)){
diff --git a/cpp/particles/particles_input_random.hpp b/cpp/particles/particles_input_random.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7d5bacefddb3c9a81fe93ecd007c7cdcad657549
--- /dev/null
+++ b/cpp/particles/particles_input_random.hpp
@@ -0,0 +1,138 @@
+/******************************************************************************
+* *
+* Copyright 2020 Max Planck Institute for Dynamics and Self-Organization *
+* *
+* This file is part of TurTLE. *
+* *
+* TurTLE 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. *
+* *
+* TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@ds.mpg.de *
+* *
+******************************************************************************/
+
+
+
+#ifndef PARTICLES_INPUT_RANDOM_HPP
+#define PARTICLES_INPUT_RANDOM_HPP
+
+#include <mpi.h>
+#include <cassert>
+#include <random>
+#include "particles/abstract_particles_input.hpp"
+
+template <class partsize_t, class particle_rnumber, int state_size>
+class particles_input_random: public abstract_particles_input<partsize_t, particle_rnumber>
+{
+ private:
+ MPI_Comm comm;
+ int myrank, nprocs;
+
+ hsize_t total_number_of_particles;
+ hsize_t number_rhs;
+ partsize_t local_number_of_particles;
+
+ std::unique_ptr<particle_rnumber[]> local_particle_state;
+ std::unique_ptr<partsize_t[]> local_particle_index;
+ std::vector<std::unique_ptr<particle_rnumber[]>> local_particle_rhs;
+ public:
+ particles_input_random(
+ const MPI_Comm in_mpi_comm,
+ const int NPARTICLES,
+ const int rseed,
+ const particle_rnumber my_spatial_low_limit,
+ const particle_rnumber my_spatial_up_limit):
+ comm(in_mpi_comm)
+ {
+ assert(state_size >= 3);
+ AssertMpi(MPI_Comm_rank(comm, &myrank));
+ AssertMpi(MPI_Comm_size(comm, &nprocs));
+ std::vector<particle_rnumber> in_spatial_limit_per_prod = BuildLimitsAllProcesses<particle_rnumber>(
+ comm, my_spatial_low_limit, my_spatial_up_limit);
+ assert(int(in_spatial_limit_per_prod.size()) == nprocs+1);
+
+ // initialize a separate random number generator for each MPI process
+ std::mt19937_64 rgen;
+ std::uniform_real_distribution<particle_rnumber> udist(0.0, 1.0);
+ std::normal_distribution<particle_rnumber> gdist;
+ // seed random number generators such that no seed is ever repeated if we change the value of rseed.
+ // basically use a multi-dimensional array indexing technique to come up with actual seed.
+ // Note: this method is not invariant to the number of MPI processes!
+ int current_seed = (
+ rseed * (this->nprocs) +
+ this->myrank);
+ rgen.seed(current_seed);
+
+ this->total_number_of_particles = NPARTICLES;
+ // we want to generate the same number of particles on all processes.
+ // this is not a uniform distribution, but it's close enough for now.
+ // in principle we should generate particles in the entire cube, then use alltoall exchange.
+ // FIXME when there's time
+ this->local_number_of_particles = NPARTICLES/this->nprocs;
+ // add remainder particles to first few MPI processes
+ // FIXME when there's time
+ if (this->myrank < (NPARTICLES%this->nprocs))
+ this->local_number_of_particles++;
+ // allocate arrays
+ this->local_particle_state.reset(new particle_rnumber[this->local_number_of_particles*state_size]);
+ this->local_particle_index.reset(new partsize_t[this->local_number_of_particles]);
+
+ const double twopi = 4*acos(0);
+
+ for (partsize_t idx=0; idx < this->local_number_of_particles; idx++)
+ {
+ this->local_particle_index[idx] = (
+ idx +
+ this->nprocs*(this->total_number_of_particles/this->nprocs) +
+ std::min(partsize_t(this->myrank), partsize_t(this->total_number_of_particles%this->nprocs)));
+ for (int cc=0; cc < 2; cc++)
+ this->local_particle_state[cc] = twopi*udist(rgen);
+ this->local_particle_state[2] = udist(rgen)*(my_spatial_up_limit - my_spatial_low_limit) + my_spatial_low_limit;
+ for (int cc = 3; cc < state_size; cc++)
+ this->local_particle_state[cc] = gdist(rgen);
+ }
+ }
+
+ partsize_t getTotalNbParticles()
+ {
+ return this->total_number_of_particles;
+ }
+ partsize_t getLocalNbParticles()
+ {
+ return this->local_number_of_particles;
+ }
+ int getNbRhs()
+ {
+ return 0;
+ }
+
+ std::unique_ptr<particle_rnumber[]> getMyParticles()
+ {
+ assert(this->local_particle_state != nullptr || this->local_number_of_particles == 0);
+ return std::move(this->local_particle_state);
+ }
+
+ std::unique_ptr<partsize_t[]> getMyParticlesIndexes()
+ {
+ assert(this->local_particle_index != nullptr || this->local_number_of_particles == 0);
+ return std::move(this->local_particle_index);
+ }
+
+ std::vector<std::unique_ptr<particle_rnumber[]>> getMyRhs()
+ {
+ return std::move(std::vector<std::unique_ptr<particle_rnumber[]>>());
+ }
+};
+
+
+#endif
diff --git a/cpp/particles/particles_system_builder.hpp b/cpp/particles/particles_system_builder.hpp
index 2e3b9fa90212602ac803a1a8de0201320e4f467e..038c5774c1fc8fe93d2d2b4bb9476dc3ce6b386c 100644
--- a/cpp/particles/particles_system_builder.hpp
+++ b/cpp/particles/particles_system_builder.hpp
@@ -257,8 +257,16 @@ struct particles_system_build_container {
// TODO P2P load particle data too
// Load particles from hdf5
- particles_input_hdf5<partsize_t, particles_rnumber, size_particle_positions, size_particle_rhs> generator(mpi_comm, fname_input,
- inDatanameState, inDatanameRhs, my_spatial_low_limit_z, my_spatial_up_limit_z);
+ particles_input_hdf5<partsize_t,
+ particles_rnumber,
+ size_particle_positions,
+ size_particle_rhs> generator(
+ mpi_comm,
+ fname_input,
+ inDatanameState,
+ inDatanameRhs,
+ my_spatial_low_limit_z,
+ my_spatial_up_limit_z);
// Ensure parameters match the input file
if(generator.getNbRhs() != nsteps){
diff --git a/cpp/particles/particles_utils.hpp b/cpp/particles/particles_utils.hpp
index 2f7e149568ff304818c8ab70e4023c6258c26f16..2446741a8166a7f87ca53eb868bff3b57d058ea5 100644
--- a/cpp/particles/particles_utils.hpp
+++ b/cpp/particles/particles_utils.hpp
@@ -348,7 +348,41 @@ public:
}
};
+}
+
+template <class real_number>
+std::vector<real_number> BuildLimitsAllProcesses(
+ MPI_Comm mpi_comm,
+ const real_number my_spatial_low_limit,
+ const real_number my_spatial_up_limit){
+ int my_rank;
+ int nb_processes;
+
+ AssertMpi(MPI_Comm_rank(mpi_comm, &my_rank));
+ AssertMpi(MPI_Comm_size(mpi_comm, &nb_processes));
+
+ std::vector<real_number> spatial_limit_per_proc(nb_processes*2);
+
+ real_number intervalToSend[2] = {my_spatial_low_limit, my_spatial_up_limit};
+ AssertMpi(
+ MPI_Allgather(
+ intervalToSend,
+ 2,
+ particles_utils::GetMpiType(real_number()),
+ spatial_limit_per_proc.data(),
+ 2,
+ particles_utils::GetMpiType(real_number()),
+ mpi_comm));
+
+ for(int idx_proc = 0; idx_proc < nb_processes-1 ; ++idx_proc){
+ assert(spatial_limit_per_proc[idx_proc*2] <= spatial_limit_per_proc[idx_proc*2+1]);
+ assert(spatial_limit_per_proc[idx_proc*2+1] == spatial_limit_per_proc[(idx_proc+1)*2]);
+ spatial_limit_per_proc[idx_proc+1] = spatial_limit_per_proc[idx_proc*2+1];
+ }
+ spatial_limit_per_proc[nb_processes] = spatial_limit_per_proc[(nb_processes-1)*2+1];
+ spatial_limit_per_proc.resize(nb_processes+1);
+ return spatial_limit_per_proc;
}
#endif