From 92985fce9a8903c4c050c42a020d557d0daf2810 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Tue, 9 Nov 2021 11:49:10 +0100
Subject: [PATCH 01/37] [wip] makes particle output more flexible
This adds a template parameter to `abstract_particles_output`, allowing
to output data of a different type than the particle's positions.
At the moment it's useless because the implementation in
`particles_output_sampling_hdf5` can only handle float or double, but
this will change in the future commit.
---
TurTLE/test/particle_set/NSVEparticle_set.cpp | 2 +-
TurTLE/test/particle_set/NSVEparticle_set.hpp | 2 +-
cpp/full_code/NSVE_Stokes_particles.cpp | 2 +-
cpp/full_code/NSVE_Stokes_particles.hpp | 2 +-
cpp/full_code/NSVEcomplex_particles.cpp | 2 +-
cpp/full_code/NSVEcomplex_particles.hpp | 4 +-
cpp/full_code/NSVEparticles.cpp | 2 +-
cpp/full_code/NSVEparticles.hpp | 4 +-
cpp/full_code/code_base.hpp | 2 +-
cpp/full_code/direct_numerical_simulation.hpp | 2 +-
cpp/full_code/kraichnan_field.cpp | 2 +-
cpp/full_code/kraichnan_field.hpp | 2 +-
cpp/full_code/static_field.cpp | 2 +-
cpp/full_code/static_field.hpp | 4 +-
cpp/full_code/test_interpolation.cpp | 2 +-
cpp/full_code/test_interpolation.hpp | 2 +-
cpp/full_code/test_tracer_set.cpp | 4 +-
cpp/particles/abstract_particle_rhs.hpp | 2 +-
cpp/particles/abstract_particles_input.hpp | 2 +-
cpp/particles/abstract_particles_output.hpp | 50 ++++++++------
cpp/particles/abstract_particles_system.hpp | 2 +-
.../interpolation/abstract_particle_set.hpp | 66 +++++++++++++++++--
cpp/particles/interpolation/particle_set.hpp | 2 +-
cpp/particles/lock_free_bool_array.hpp | 4 +-
cpp/particles/particles_distr_mpi.hpp | 2 +-
cpp/particles/particles_input_hdf5.hpp | 2 +-
cpp/particles/particles_output_hdf5.hpp | 5 +-
cpp/particles/particles_output_mpiio.hpp | 8 +--
.../particles_output_sampling_hdf5.hpp | 27 ++++----
cpp/particles/particles_sampling.hpp | 8 +--
cpp/particles/particles_system.hpp | 2 +-
31 files changed, 146 insertions(+), 78 deletions(-)
diff --git a/TurTLE/test/particle_set/NSVEparticle_set.cpp b/TurTLE/test/particle_set/NSVEparticle_set.cpp
index 8a32ebc0..4e98159a 100644
--- a/TurTLE/test/particle_set/NSVEparticle_set.cpp
+++ b/TurTLE/test/particle_set/NSVEparticle_set.cpp
@@ -87,7 +87,7 @@ int NSVEparticle_set<rnumber>::initialize(void)
nparticles,
tracers0_integration_steps);
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->pset->getTotalNumberOfParticles(),
(this->simname + "_particles.h5"),
diff --git a/TurTLE/test/particle_set/NSVEparticle_set.hpp b/TurTLE/test/particle_set/NSVEparticle_set.hpp
index 71f58c81..4e297d34 100644
--- a/TurTLE/test/particle_set/NSVEparticle_set.hpp
+++ b/TurTLE/test/particle_set/NSVEparticle_set.hpp
@@ -76,7 +76,7 @@ class NSVEparticle_set: public NSVE<rnumber>
field<rnumber, FFTW, THREE> *previous_velocity;
particles_output_hdf5<long long int, double, 3> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
NSVEparticle_set(
const MPI_Comm COMMUNICATOR,
diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp
index c77b9103..29eeca39 100644
--- a/cpp/full_code/NSVE_Stokes_particles.cpp
+++ b/cpp/full_code/NSVE_Stokes_particles.cpp
@@ -70,7 +70,7 @@ int NSVE_Stokes_particles<rnumber>::initialize(void)
tracers0_integration_steps);
this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/NSVE_Stokes_particles.hpp b/cpp/full_code/NSVE_Stokes_particles.hpp
index ecc2d135..fb8497c9 100644
--- a/cpp/full_code/NSVE_Stokes_particles.hpp
+++ b/cpp/full_code/NSVE_Stokes_particles.hpp
@@ -65,7 +65,7 @@ class NSVE_Stokes_particles: public NSVE<rnumber>
field<rnumber, FFTW, ONE> *pressure;
particles_output_hdf5<long long int, double,6> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
NSVE_Stokes_particles(
diff --git a/cpp/full_code/NSVEcomplex_particles.cpp b/cpp/full_code/NSVEcomplex_particles.cpp
index a242cc64..7b37476d 100644
--- a/cpp/full_code/NSVEcomplex_particles.cpp
+++ b/cpp/full_code/NSVEcomplex_particles.cpp
@@ -68,7 +68,7 @@ int NSVEcomplex_particles<rnumber>::initialize(void)
nparticles,
tracers0_integration_steps);
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/NSVEcomplex_particles.hpp b/cpp/full_code/NSVEcomplex_particles.hpp
index a63995f9..3f6adaac 100644
--- a/cpp/full_code/NSVEcomplex_particles.hpp
+++ b/cpp/full_code/NSVEcomplex_particles.hpp
@@ -70,7 +70,7 @@ class NSVEcomplex_particles: public NSVE<rnumber>
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
// TODO P2P use a reader with particle data
particles_output_hdf5<long long int, double,6> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
// field for sampling velocity gradient
field<rnumber, FFTW, THREExTHREE> *nabla_u;
@@ -82,7 +82,7 @@ class NSVEcomplex_particles: public NSVE<rnumber>
COMMUNICATOR,
simulation_name),
cutoff(10), inner_v0(1), lambda(1.0), enable_p2p(true), enable_inner(true), enable_vorticity_omega(true){}
- ~NSVEcomplex_particles(){}
+ ~NSVEcomplex_particles() noexcept(false){}
int initialize(void);
int step(void);
diff --git a/cpp/full_code/NSVEparticles.cpp b/cpp/full_code/NSVEparticles.cpp
index 170f0046..ab05dc41 100644
--- a/cpp/full_code/NSVEparticles.cpp
+++ b/cpp/full_code/NSVEparticles.cpp
@@ -115,7 +115,7 @@ int NSVEparticles<rnumber>::initialize(void)
// initialize sample object
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/NSVEparticles.hpp b/cpp/full_code/NSVEparticles.hpp
index bbb8ca9f..dfa5598a 100644
--- a/cpp/full_code/NSVEparticles.hpp
+++ b/cpp/full_code/NSVEparticles.hpp
@@ -65,7 +65,7 @@ class NSVEparticles: public NSVE<rnumber>
field<rnumber, FFTW, ONE> *pressure;
particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
NSVEparticles(
@@ -74,7 +74,7 @@ class NSVEparticles: public NSVE<rnumber>
NSVE<rnumber>(
COMMUNICATOR,
simulation_name){}
- ~NSVEparticles(){}
+ ~NSVEparticles() noexcept(false){}
int initialize(void);
int step(void);
diff --git a/cpp/full_code/code_base.hpp b/cpp/full_code/code_base.hpp
index bd9c30e1..b794347a 100644
--- a/cpp/full_code/code_base.hpp
+++ b/cpp/full_code/code_base.hpp
@@ -71,7 +71,7 @@ class code_base
code_base(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name);
- virtual ~code_base(){}
+ virtual ~code_base() noexcept(false){}
int check_stopping_condition(void);
diff --git a/cpp/full_code/direct_numerical_simulation.hpp b/cpp/full_code/direct_numerical_simulation.hpp
index 668fa53f..87f19078 100644
--- a/cpp/full_code/direct_numerical_simulation.hpp
+++ b/cpp/full_code/direct_numerical_simulation.hpp
@@ -50,7 +50,7 @@ class direct_numerical_simulation: public code_base
code_base(
COMMUNICATOR,
simulation_name){}
- virtual ~direct_numerical_simulation(){}
+ virtual ~direct_numerical_simulation() noexcept(false){}
virtual int read_parameters(void);
virtual int write_checkpoint(void) = 0;
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 661ed6fd..26887a72 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -93,7 +93,7 @@ int kraichnan_field<rnumber>::initialize(void)
tracers0_integration_steps);
this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
index 0a9d4a3f..020cd0cd 100644
--- a/cpp/full_code/kraichnan_field.hpp
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -72,7 +72,7 @@ class kraichnan_field: public direct_numerical_simulation
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
kraichnan_field(
diff --git a/cpp/full_code/static_field.cpp b/cpp/full_code/static_field.cpp
index 48396bc0..118721a9 100644
--- a/cpp/full_code/static_field.cpp
+++ b/cpp/full_code/static_field.cpp
@@ -117,7 +117,7 @@ int static_field<rnumber>::initialize(void)
tracers0_integration_steps);
this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/static_field.hpp b/cpp/full_code/static_field.hpp
index 6393ae21..90b563d1 100644
--- a/cpp/full_code/static_field.hpp
+++ b/cpp/full_code/static_field.hpp
@@ -65,7 +65,7 @@ class static_field: public direct_numerical_simulation
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
static_field(
@@ -74,7 +74,7 @@ class static_field: public direct_numerical_simulation
direct_numerical_simulation(
COMMUNICATOR,
simulation_name){}
- ~static_field(){}
+ ~static_field() noexcept(false){}
int initialize(void);
int step(void);
diff --git a/cpp/full_code/test_interpolation.cpp b/cpp/full_code/test_interpolation.cpp
index 2bbb1f6e..fd17ca4f 100644
--- a/cpp/full_code/test_interpolation.cpp
+++ b/cpp/full_code/test_interpolation.cpp
@@ -100,7 +100,7 @@ int test_interpolation<rnumber>::initialize(void)
nparticles,
this->tracers0_integration_steps);
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
this->ps->getGlobalNbParticles(),
(this->simname + "_output.h5"),
diff --git a/cpp/full_code/test_interpolation.hpp b/cpp/full_code/test_interpolation.hpp
index e312eddc..5debb710 100644
--- a/cpp/full_code/test_interpolation.hpp
+++ b/cpp/full_code/test_interpolation.hpp
@@ -52,7 +52,7 @@ class test_interpolation: public test
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
- particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
field<rnumber, FFTW, THREE> *velocity, *vorticity;
field<rnumber, FFTW, THREExTHREE> *nabla_u;
diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
index df779411..5c2b1fc5 100644
--- a/cpp/full_code/test_tracer_set.cpp
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -108,8 +108,8 @@ int test_tracer_set<rnumber>::do_work(void)
"tracers0",
pset.getTotalNumberOfParticles(),
0);
- particles_output_sampling_hdf5<long long int, double, 3> *particle_sample_writer = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particle_sample_writer = new particles_output_sampling_hdf5<
+ long long int, double, double, 3>(
MPI_COMM_WORLD,
pset.getTotalNumberOfParticles(),
"test_particle_sample.h5",
diff --git a/cpp/particles/abstract_particle_rhs.hpp b/cpp/particles/abstract_particle_rhs.hpp
index d1a47bdc..947e1524 100644
--- a/cpp/particles/abstract_particle_rhs.hpp
+++ b/cpp/particles/abstract_particle_rhs.hpp
@@ -45,7 +45,7 @@ class abstract_particle_rhs
using particle_rnumber = abstract_particle_set::particle_rnumber;
public:
// destructor
- virtual ~abstract_particle_rhs(){}
+ virtual ~abstract_particle_rhs() noexcept(false){}
/** Probe the dimension of the ODE
*/
diff --git a/cpp/particles/abstract_particles_input.hpp b/cpp/particles/abstract_particles_input.hpp
index ab078a42..2b581b61 100644
--- a/cpp/particles/abstract_particles_input.hpp
+++ b/cpp/particles/abstract_particles_input.hpp
@@ -31,7 +31,7 @@
template <class partsize_t, class real_number>
class abstract_particles_input {
public:
- virtual ~abstract_particles_input(){}
+ virtual ~abstract_particles_input() noexcept(false){}
virtual partsize_t getTotalNbParticles() = 0;
virtual partsize_t getLocalNbParticles() = 0;
diff --git a/cpp/particles/abstract_particles_output.hpp b/cpp/particles/abstract_particles_output.hpp
index 26749461..de926680 100644
--- a/cpp/particles/abstract_particles_output.hpp
+++ b/cpp/particles/abstract_particles_output.hpp
@@ -38,7 +38,15 @@
#include "scope_timer.hpp"
#include "env_utils.hpp"
-template <class partsize_t, class real_number, int size_particle_positions>
+/** \brief Class to handle distributed output of particle data
+ *
+ * \tparam partsize_t data type of particle index
+ * \tparam position_type data type of particle positions
+ * \tparam output_type data type of output (typically same as position_type, but it may be different)
+ * \tparam size_particle_positions int, size of particle state
+ *
+ */
+template <class partsize_t, class position_type, class output_type, int size_particle_positions>
class abstract_particles_output {
MPI_Comm mpi_com;
MPI_Comm mpi_com_writer;
@@ -50,13 +58,13 @@ class abstract_particles_output {
const int nb_rhs;
std::unique_ptr<std::pair<partsize_t,partsize_t>[]> buffer_indexes_send;
- std::unique_ptr<real_number[]> buffer_particles_positions_send;
- std::vector<std::unique_ptr<real_number[]>> buffer_particles_rhs_send;
+ std::unique_ptr<position_type[]> buffer_particles_positions_send;
+ std::vector<std::unique_ptr<output_type[]>> buffer_particles_rhs_send;
partsize_t size_buffers_send;
int buffers_size_particle_rhs_send;
- std::unique_ptr<real_number[]> buffer_particles_positions_recv;
- std::vector<std::unique_ptr<real_number[]>> buffer_particles_rhs_recv;
+ std::unique_ptr<position_type[]> buffer_particles_positions_recv;
+ std::vector<std::unique_ptr<output_type[]>> buffer_particles_rhs_recv;
std::unique_ptr<partsize_t[]> buffer_indexes_recv;
partsize_t size_buffers_recv;
int buffers_size_particle_rhs_recv;
@@ -107,7 +115,7 @@ public:
const int MaxProcessesInvolved = std::min(nb_processes, env_utils::GetValue<int>("BFPS_PO_MAX_PROCESSES", 128));
// We split the processes using positions size only
- const size_t totalBytesForPositions = total_nb_particles*size_particle_positions*sizeof(real_number);
+ const size_t totalBytesForPositions = total_nb_particles*size_particle_positions*sizeof(position_type);
if(MinBytesPerProcess*MaxProcessesInvolved < totalBytesForPositions){
@@ -116,13 +124,13 @@ public:
extraChunkBytes += 1;
}
const size_t bytesPerProcess = (MinBytesPerProcess+extraChunkBytes*ChunkBytes);
- particles_chunk_per_process = partsize_t((bytesPerProcess+sizeof(real_number)*size_particle_positions-1)/(sizeof(real_number)*size_particle_positions));
+ particles_chunk_per_process = partsize_t((bytesPerProcess+sizeof(position_type)*size_particle_positions-1)/(sizeof(position_type)*size_particle_positions));
nb_processes_involved = int((total_nb_particles+particles_chunk_per_process-1)/particles_chunk_per_process);
}
// else limit based on minBytesPerProcess
else{
nb_processes_involved = std::max(1,std::min(MaxProcessesInvolved,int((totalBytesForPositions+MinBytesPerProcess-1)/MinBytesPerProcess)));
- particles_chunk_per_process = partsize_t((MinBytesPerProcess+sizeof(real_number)*size_particle_positions-1)/(sizeof(real_number)*size_particle_positions));
+ particles_chunk_per_process = partsize_t((MinBytesPerProcess+sizeof(position_type)*size_particle_positions-1)/(sizeof(position_type)*size_particle_positions));
}
// Print out
@@ -178,8 +186,8 @@ public:
template <int size_particle_rhs>
void save(
- const real_number input_particles_positions[],
- const std::unique_ptr<real_number[]> input_particles_rhs[],
+ const position_type input_particles_positions[],
+ const std::unique_ptr<output_type[]> input_particles_rhs[],
const partsize_t index_particles[],
const partsize_t nb_particles,
const int idx_time_step){
@@ -192,20 +200,20 @@ public:
if(size_buffers_send < nb_particles){
size_buffers_send = nb_particles;
buffer_indexes_send.reset(new std::pair<partsize_t,partsize_t>[size_buffers_send]);
- buffer_particles_positions_send.reset(new real_number[size_buffers_send*size_particle_positions]);
+ buffer_particles_positions_send.reset(new position_type[size_buffers_send*size_particle_positions]);
if(buffers_size_particle_rhs_send < size_particle_rhs){
buffers_size_particle_rhs_send = size_particle_rhs;
}
for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
- buffer_particles_rhs_send[idx_rhs].reset(new real_number[size_buffers_send*buffers_size_particle_rhs_send]);
+ buffer_particles_rhs_send[idx_rhs].reset(new output_type[size_buffers_send*buffers_size_particle_rhs_send]);
}
}
else if(buffers_size_particle_rhs_send < size_particle_rhs){
buffers_size_particle_rhs_send = size_particle_rhs;
if(size_buffers_send > 0){
for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
- buffer_particles_rhs_send[idx_rhs].reset(new real_number[size_buffers_send*buffers_size_particle_rhs_send]);
+ buffer_particles_rhs_send[idx_rhs].reset(new output_type[size_buffers_send*buffers_size_particle_rhs_send]);
}
}
}
@@ -255,18 +263,18 @@ public:
if(size_buffers_recv < nb_to_receive){
size_buffers_recv = nb_to_receive;
buffer_indexes_recv.reset(new partsize_t[size_buffers_recv]);
- buffer_particles_positions_recv.reset(new real_number[size_buffers_recv*size_particle_positions]);
+ buffer_particles_positions_recv.reset(new position_type[size_buffers_recv*size_particle_positions]);
buffers_size_particle_rhs_recv = size_particle_rhs;
for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
- buffer_particles_rhs_recv[idx_rhs].reset(new real_number[size_buffers_recv*buffers_size_particle_rhs_recv]);
+ buffer_particles_rhs_recv[idx_rhs].reset(new output_type[size_buffers_recv*buffers_size_particle_rhs_recv]);
}
}
else if(buffers_size_particle_rhs_recv < size_particle_rhs){
buffers_size_particle_rhs_recv = size_particle_rhs;
if(size_buffers_recv > 0){
for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
- buffer_particles_rhs_recv[idx_rhs].reset(new real_number[size_buffers_recv*buffers_size_particle_rhs_recv]);
+ buffer_particles_rhs_recv[idx_rhs].reset(new output_type[size_buffers_recv*buffers_size_particle_rhs_recv]);
}
}
}
@@ -275,9 +283,9 @@ public:
TIMEZONE("exchange");
// Could be done with multiple asynchronous coms
exchanger.alltoallv<partsize_t>(buffer_indexes_send_tmp, buffer_indexes_recv.get());
- exchanger.alltoallv<real_number>(buffer_particles_positions_send.get(), buffer_particles_positions_recv.get(), size_particle_positions);
+ exchanger.alltoallv<position_type>(buffer_particles_positions_send.get(), buffer_particles_positions_recv.get(), size_particle_positions);
for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
- exchanger.alltoallv<real_number>(buffer_particles_rhs_send[idx_rhs].get(), buffer_particles_rhs_recv[idx_rhs].get(), size_particle_rhs);
+ exchanger.alltoallv<output_type>(buffer_particles_rhs_send[idx_rhs].get(), buffer_particles_rhs_recv[idx_rhs].get(), size_particle_rhs);
}
}
@@ -290,10 +298,10 @@ public:
if(size_buffers_send < nb_to_receive){
size_buffers_send = nb_to_receive;
buffer_indexes_send.reset(new std::pair<partsize_t,partsize_t>[size_buffers_send]);
- buffer_particles_positions_send.reset(new real_number[size_buffers_send*size_particle_positions]);
+ buffer_particles_positions_send.reset(new position_type[size_buffers_send*size_particle_positions]);
buffers_size_particle_rhs_send = size_particle_rhs;
for(int idx_rhs = 0 ; idx_rhs < nb_rhs ; ++idx_rhs){
- buffer_particles_rhs_send[idx_rhs].reset(new real_number[size_buffers_send*buffers_size_particle_rhs_send]);
+ buffer_particles_rhs_send[idx_rhs].reset(new output_type[size_buffers_send*buffers_size_particle_rhs_send]);
}
}
@@ -322,7 +330,7 @@ public:
nb_to_receive, particles_chunk_current_offset, size_particle_rhs);
}
- virtual void write(const int idx_time_step, const real_number* positions, const std::unique_ptr<real_number[]>* rhs,
+ virtual void write(const int idx_time_step, const position_type* positions, const std::unique_ptr<output_type[]>* rhs,
const partsize_t nb_particles, const partsize_t particles_idx_offset, const int size_particle_rhs) = 0;
};
diff --git a/cpp/particles/abstract_particles_system.hpp b/cpp/particles/abstract_particles_system.hpp
index abcb8e68..4edf20c7 100644
--- a/cpp/particles/abstract_particles_system.hpp
+++ b/cpp/particles/abstract_particles_system.hpp
@@ -38,7 +38,7 @@
template <class partsize_t, class real_number>
class abstract_particles_system {
public:
- virtual ~abstract_particles_system(){}
+ virtual ~abstract_particles_system() noexcept(false){}
virtual void delete_particles_from_indexes(const std::vector<partsize_t>& inIndexToDelete) = 0;
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index 085956f4..a8343b86 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -30,13 +30,14 @@
#include <cassert>
#include <vector>
#include <memory>
+#include <string>
#include "field.hpp"
#include "particles/p2p/p2p_distr_mpi.hpp"
#include "particles/particles_sampling.hpp"
-/** Brings together particle information with interpolation functionality.
+/** \brief Brings together particle information with interpolation functionality.
*
* This is an abstract class that defines the functionality required by the
* particle solver code.
@@ -56,7 +57,7 @@ class abstract_particle_set
using particle_rnumber = double;
// virtual destructor
- virtual ~abstract_particle_set(){}
+ virtual ~abstract_particle_set() noexcept(false){}
// extract particle information
virtual particle_rnumber* getParticleState() const = 0;
@@ -187,12 +188,56 @@ class abstract_particle_set
return EXIT_SUCCESS;
}
+ template <typename field_rnumber,
+ field_backend be,
+ field_components fc>
+ int writeParticleIndex(
+ const std::string file_name,
+ const std::string species_name,
+ const std::string field_name,
+ const int iteration)
+ {
+ TIMEZONE("abstract_particle_set::writeParticleIndex");
+ particles_output_sampling_hdf5<partsize_t, particle_rnumber, partsize_t, 3> *particle_sample_writer = new particles_output_sampling_hdf5<partsize_t, double, partsize_t, 3>(
+ MPI_COMM_WORLD,
+ this->getTotalNumberOfParticles(),
+ file_name,
+ species_name,
+ field_name + std::string("/") + std::to_string(iteration));
+ // set file layout
+ particle_sample_writer->setParticleFileLayout(this->getParticleFileLayout());
+ // allocate position array
+ std::unique_ptr<particle_rnumber[]> xx = this->extractFromParticleState(0, 3);
+ // allocate temporary array
+ std::unique_ptr<partsize_t[]> pdata(new partsize_t[this->getLocalNumberOfParticles()]);
+ // clean up temporary array
+ std::copy(this->getParticleIndex(),
+ this->getParticleIndex() + this->getLocalNumberOfParticles(),
+ pdata.get());
+ particle_sample_writer->template save_dataset<1>(
+ species_name,
+ field_name,
+ xx.get(),
+ &pdata,
+ this->getParticleIndex(),
+ this->getLocalNumberOfParticles(),
+ iteration);
+ // deallocate temporary array
+ delete[] pdata.release();
+ // deallocate position array
+ delete[] xx.release();
+ return EXIT_SUCCESS;
+ }
+
template <typename field_rnumber,
field_backend be,
field_components fc>
int writeSample(
field<field_rnumber, be, fc> *field_to_sample,
- particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
+ particles_output_sampling_hdf5<partsize_t,
+ particle_rnumber,
+ particle_rnumber,
+ 3> *particle_sample_writer,
const std::string species_name,
const std::string field_name,
const int iteration)
@@ -257,7 +302,10 @@ class abstract_particle_set
int writeStateChunk(
const int i0,
const int contiguous_state_chunk,
- particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
+ particles_output_sampling_hdf5<partsize_t,
+ particle_rnumber,
+ particle_rnumber,
+ 3> *particle_sample_writer,
const std::string species_name,
const std::string field_name,
const int iteration)
@@ -316,7 +364,10 @@ class abstract_particle_set
int writeStateTriplet(
const int i0,
- particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
+ particles_output_sampling_hdf5<partsize_t,
+ particle_rnumber,
+ particle_rnumber,
+ 3> *particle_sample_writer,
const std::string species_name,
const std::string field_name,
const int iteration)
@@ -332,7 +383,10 @@ class abstract_particle_set
int writeStateComponent(
const int i0,
- particles_output_sampling_hdf5<partsize_t, particle_rnumber, 3> *particle_sample_writer,
+ particles_output_sampling_hdf5<partsize_t,
+ particle_rnumber,
+ particle_rnumber,
+ 3> *particle_sample_writer,
const std::string species_name,
const std::string field_name,
const int iteration)
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 206e9da5..9abb2c1b 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -167,7 +167,7 @@ class particle_set: public abstract_particle_set
this->offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
}
- ~particle_set(){}
+ ~particle_set() noexcept(false){}
particle_rnumber* getParticleState() const
{
diff --git a/cpp/particles/lock_free_bool_array.hpp b/cpp/particles/lock_free_bool_array.hpp
index bfcd4888..9c0ab5de 100644
--- a/cpp/particles/lock_free_bool_array.hpp
+++ b/cpp/particles/lock_free_bool_array.hpp
@@ -42,7 +42,7 @@ class lock_free_bool_array{
Locker(){
omp_init_nest_lock(&lock);
}
- ~Locker(){
+ ~Locker() noexcept(false){
omp_destroy_nest_lock(&lock);
}
omp_nest_lock_t lock;
@@ -68,7 +68,7 @@ public:
}
}
#ifndef NDEBUG
- ~lock_free_bool_array(){
+ ~lock_free_bool_array() noexcept(false){
for(auto& k : keys){
#ifdef __INTEL_COMPILER
#else
diff --git a/cpp/particles/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index b7d0f9f8..1681ea3d 100644
--- a/cpp/particles/particles_distr_mpi.hpp
+++ b/cpp/particles/particles_distr_mpi.hpp
@@ -160,7 +160,7 @@ public:
assert(int(field_grid_dim[IDXC_Z]) == partition_interval_offset_per_proc[nb_processes_involved]);
}
- virtual ~particles_distr_mpi(){}
+ virtual ~particles_distr_mpi() noexcept(false){}
////////////////////////////////////////////////////////////////////////////
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 680daff9..28740ba8 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -320,7 +320,7 @@ public:
}
}
- ~particles_input_hdf5(){
+ ~particles_input_hdf5() noexcept(false){
}
partsize_t getTotalNbParticles() final{
diff --git a/cpp/particles/particles_output_hdf5.hpp b/cpp/particles/particles_output_hdf5.hpp
index 236504e1..0ef809e6 100644
--- a/cpp/particles/particles_output_hdf5.hpp
+++ b/cpp/particles/particles_output_hdf5.hpp
@@ -38,9 +38,11 @@ template <class partsize_t,
class real_number,
int size_particle_positions>
class particles_output_hdf5 : public abstract_particles_output<partsize_t,
+ real_number,
real_number,
size_particle_positions>{
using Parent = abstract_particles_output<partsize_t,
+ real_number,
real_number,
size_particle_positions>;
@@ -62,6 +64,7 @@ public:
const int in_nb_rhs,
const bool in_use_collective_io = false)
: abstract_particles_output<partsize_t,
+ real_number,
real_number,
size_particle_positions>(
in_mpi_com,
@@ -112,7 +115,7 @@ public:
return EXIT_SUCCESS;
}
- ~particles_output_hdf5(){}
+ ~particles_output_hdf5() noexcept(false){}
void update_particle_species_name(
const std::string new_name)
diff --git a/cpp/particles/particles_output_mpiio.hpp b/cpp/particles/particles_output_mpiio.hpp
index b05579e2..ee9bd9ce 100644
--- a/cpp/particles/particles_output_mpiio.hpp
+++ b/cpp/particles/particles_output_mpiio.hpp
@@ -36,8 +36,8 @@
#include "particles_utils.hpp"
template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
-class particles_output_mpiio : public abstract_particles_output<partsize_t, real_number, size_particle_positions>{
- using Parent = abstract_particles_output<partsize_t, real_number, size_particle_positions>;
+class particles_output_mpiio : public abstract_particles_output<partsize_t, real_number, real_number, size_particle_positions>{
+ using Parent = abstract_particles_output<partsize_t, real_number, real_number, size_particle_positions>;
const std::string filename;
const int nb_step_prealloc;
@@ -49,7 +49,7 @@ class particles_output_mpiio : public abstract_particles_output<partsize_t, real
public:
particles_output_mpiio(MPI_Comm in_mpi_com, const std::string in_filename, const partsize_t inTotalNbParticles,
const int in_nb_rhs, const int in_nb_step_prealloc = -1)
- : abstract_particles_output<partsize_t, real_number, size_particle_positions>(in_mpi_com, inTotalNbParticles, in_nb_rhs),
+ : abstract_particles_output<partsize_t, real_number, real_number, size_particle_positions>(in_mpi_com, inTotalNbParticles, in_nb_rhs),
filename(in_filename), nb_step_prealloc(in_nb_step_prealloc), current_step_in_file(0){
if(Parent::isInvolved()){
{
@@ -65,7 +65,7 @@ public:
}
}
- ~particles_output_mpiio(){
+ ~particles_output_mpiio() noexcept(false){
if(Parent::isInvolved()){
TIMEZONE("particles_output_mpiio::MPI_File_close");
AssertMpi(MPI_File_close(&mpi_file));
diff --git a/cpp/particles/particles_output_sampling_hdf5.hpp b/cpp/particles/particles_output_sampling_hdf5.hpp
index 11a7dcc1..0e0d1fc6 100644
--- a/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -31,14 +31,17 @@
#include <hdf5.h>
template <class partsize_t,
- class real_number,
+ class position_type,
+ class output_type,
int size_particle_positions>
class particles_output_sampling_hdf5 : public abstract_particles_output<
partsize_t,
- real_number,
+ position_type,
+ output_type,
size_particle_positions>{
using Parent = abstract_particles_output<partsize_t,
- real_number,
+ position_type,
+ output_type,
size_particle_positions>;
hid_t file_id, pgroup_id;
@@ -116,7 +119,7 @@ public:
}
}
- ~particles_output_sampling_hdf5(){
+ ~particles_output_sampling_hdf5() noexcept(false){
if(Parent::isInvolved()){
// close group
int retTest = H5Gclose(pgroup_id);
@@ -151,8 +154,8 @@ public:
int save_dataset(
const std::string& in_groupname,
const std::string& in_dataset_name,
- const real_number input_particles_positions[],
- const std::unique_ptr<real_number[]> input_particles_rhs[],
+ const position_type input_particles_positions[],
+ const std::unique_ptr<output_type[]> input_particles_rhs[],
const partsize_t index_particles[],
const partsize_t nb_particles,
const int idx_time_step)
@@ -183,8 +186,8 @@ public:
void write(
const int /*idx_time_step*/,
- const real_number* /*particles_positions*/,
- const std::unique_ptr<real_number[]>* particles_rhs,
+ const position_type* /*particles_positions*/,
+ const std::unique_ptr<output_type[]>* particles_rhs,
const partsize_t nb_particles,
const partsize_t particles_idx_offset,
const int size_particle_rhs) final{
@@ -197,10 +200,10 @@ public:
nb_particles == 0));
assert(particles_idx_offset+nb_particles <= Parent::getTotalNbParticles());
- static_assert(std::is_same<real_number, double>::value ||
- std::is_same<real_number, float>::value,
- "real_number must be double or float");
- const hid_t type_id = (sizeof(real_number) == 8 ?
+ static_assert(std::is_same<output_type, double>::value ||
+ std::is_same<output_type, float>::value,
+ "output_type must be double or float");
+ const hid_t type_id = (sizeof(output_type) == 8 ?
H5T_NATIVE_DOUBLE :
H5T_NATIVE_FLOAT);
diff --git a/cpp/particles/particles_sampling.hpp b/cpp/particles/particles_sampling.hpp
index 1606b022..8e83a5a5 100644
--- a/cpp/particles/particles_sampling.hpp
+++ b/cpp/particles/particles_sampling.hpp
@@ -46,7 +46,7 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
const int size_particle_rhs = ncomp(fc);
// Stop here if already exists
- if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3>::DatasetExistsCol(MPI_COMM_WORLD,
+ if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, particles_rnumber, 3>::DatasetExistsCol(MPI_COMM_WORLD,
filename,
parent_groupname,
datasetname)){
@@ -61,7 +61,7 @@ void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a p
- particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3> outputclass(MPI_COMM_WORLD,
+ particles_output_sampling_hdf5<partsize_t, particles_rnumber, particles_rnumber, 3> outputclass(MPI_COMM_WORLD,
ps->getGlobalNbParticles(),
filename,
parent_groupname,
@@ -82,7 +82,7 @@ void sample_particles_system_position(
const std::string datasetname = fname + std::string("/") + std::to_string(ps->get_step_idx());
// Stop here if already exists
- if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3>::DatasetExistsCol(MPI_COMM_WORLD,
+ if(particles_output_sampling_hdf5<partsize_t, particles_rnumber, particles_rnumber, 3>::DatasetExistsCol(MPI_COMM_WORLD,
filename,
parent_groupname,
datasetname)){
@@ -93,7 +93,7 @@ void sample_particles_system_position(
std::unique_ptr<particles_rnumber[]> sample_rhs(new particles_rnumber[3*nb_particles]);
std::copy(ps->getParticlesState(), ps->getParticlesState() + 3*nb_particles, sample_rhs.get());
- particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3> outputclass(MPI_COMM_WORLD,
+ particles_output_sampling_hdf5<partsize_t, particles_rnumber, particles_rnumber, 3> outputclass(MPI_COMM_WORLD,
ps->getGlobalNbParticles(),
filename,
parent_groupname,
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index b4776302..9589e4dc 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -137,7 +137,7 @@ public:
current_offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
}
- ~particles_system(){
+ ~particles_system() noexcept(false){
}
void init(abstract_particles_input<partsize_t, real_number>& particles_input) {
--
GitLab
From 7013fdd2ecbb0afc50bb236ffcd03ba692919551 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Tue, 9 Nov 2021 16:52:21 +0100
Subject: [PATCH 02/37] cleans up HDF5 type selection
---
cpp/particles/particles_input_hdf5.hpp | 2 +-
cpp/particles/particles_output_sampling_hdf5.hpp | 9 ++++-----
2 files changed, 5 insertions(+), 6 deletions(-)
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 28740ba8..3cb04ce3 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -166,7 +166,7 @@ public:
static_assert(std::is_same<real_number, double>::value
|| std::is_same<real_number, float>::value, "real_number must be double or float");
- const hid_t type_id = (sizeof(real_number) == 8?H5T_NATIVE_DOUBLE:H5T_NATIVE_FLOAT);
+ const hid_t type_id = hdf5_tools::hdf5_type_id<real_number>();
/// Load the data
std::unique_ptr<real_number[]> split_particles_positions;
diff --git a/cpp/particles/particles_output_sampling_hdf5.hpp b/cpp/particles/particles_output_sampling_hdf5.hpp
index 0e0d1fc6..73bc253f 100644
--- a/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -201,11 +201,10 @@ public:
assert(particles_idx_offset+nb_particles <= Parent::getTotalNbParticles());
static_assert(std::is_same<output_type, double>::value ||
- std::is_same<output_type, float>::value,
- "output_type must be double or float");
- const hid_t type_id = (sizeof(output_type) == 8 ?
- H5T_NATIVE_DOUBLE :
- H5T_NATIVE_FLOAT);
+ std::is_same<output_type, float>::value ||
+ std::is_same<output_type, long long int>::value,
+ "output_type must be double or float or long long int");
+ const hid_t type_id = hdf5_tools::hdf5_type_id<output_type>();
hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
assert(plist_id >= 0);
--
GitLab
From 8f0a8f83ebe6d4f382e113650d537b926bf470f8 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 10 Nov 2021 07:16:39 +0100
Subject: [PATCH 03/37] fixes license header
---
TurTLE/test/particle_set/NSVEparticle_set.cpp | 8 ++++----
TurTLE/test/particle_set/NSVEparticle_set.hpp | 8 ++++----
TurTLE/test/profiler/Kraichnan_scalar_v1.cpp | 8 ++++----
TurTLE/test/profiler/Kraichnan_scalar_v1.hpp | 8 ++++----
TurTLE/test/profiler/scalar_evolution.cpp | 8 ++++----
TurTLE/test/profiler/scalar_evolution.hpp | 8 ++++----
6 files changed, 24 insertions(+), 24 deletions(-)
diff --git a/TurTLE/test/particle_set/NSVEparticle_set.cpp b/TurTLE/test/particle_set/NSVEparticle_set.cpp
index 4e98159a..e2196a85 100644
--- a/TurTLE/test/particle_set/NSVEparticle_set.cpp
+++ b/TurTLE/test/particle_set/NSVEparticle_set.cpp
@@ -3,20 +3,20 @@
* Copyright 2019 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* 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. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/TurTLE/test/particle_set/NSVEparticle_set.hpp b/TurTLE/test/particle_set/NSVEparticle_set.hpp
index 4e297d34..7d07bbf8 100644
--- a/TurTLE/test/particle_set/NSVEparticle_set.hpp
+++ b/TurTLE/test/particle_set/NSVEparticle_set.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* 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. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp b/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
index 1b33f36f..7e9a647d 100644
--- a/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
+++ b/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* 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. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp b/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp
index cb502906..97aa9cdb 100644
--- a/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp
+++ b/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* 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. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/TurTLE/test/profiler/scalar_evolution.cpp b/TurTLE/test/profiler/scalar_evolution.cpp
index fa52286e..5e112ea3 100644
--- a/TurTLE/test/profiler/scalar_evolution.cpp
+++ b/TurTLE/test/profiler/scalar_evolution.cpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* 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. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/TurTLE/test/profiler/scalar_evolution.hpp b/TurTLE/test/profiler/scalar_evolution.hpp
index 48eca028..0eea8131 100644
--- a/TurTLE/test/profiler/scalar_evolution.hpp
+++ b/TurTLE/test/profiler/scalar_evolution.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* 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. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
--
GitLab
From a0aa6e59f65d6a50e95a11984699d1db670626d5 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sun, 28 Nov 2021 10:59:40 +0100
Subject: [PATCH 04/37] adds plain test of particle index output
---
CMakeLists.txt | 5 +++++
cpp/particles/interpolation/abstract_particle_set.hpp | 3 ---
setup.py.in | 3 ++-
3 files changed, 7 insertions(+), 4 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 15b83e02..96b3530c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -556,6 +556,11 @@ if (BUILD_TESTING)
NAME test_Heun_p2p
COMMAND turtle.test_Heun_p2p -n 32 --np 2 --ntpp 2 --simname tracer_set_Heun_p2p --src-simname dns_nsveparticles --src-iteration 0
WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ ### check whether particle deletion works
+ add_test(
+ NAME test_particle_deleter
+ COMMAND turtle.test_particle_deleter -n 32 --np 2 --ntpp 2 --simname tracer_set_particle_deleter --src-simname dns_nsveparticles --src-iteration 0
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
### simple runs of post-processing tools
add_test(
NAME test_pp_single_to_double
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index a8343b86..eea52c17 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -188,9 +188,6 @@ class abstract_particle_set
return EXIT_SUCCESS;
}
- template <typename field_rnumber,
- field_backend be,
- field_components fc>
int writeParticleIndex(
const std::string file_name,
const std::string species_name,
diff --git a/setup.py.in b/setup.py.in
index 369f5143..548534e7 100644
--- a/setup.py.in
+++ b/setup.py.in
@@ -47,7 +47,8 @@ setup(
'turtle.test_particles = TurTLE.test.test_particles:main',
'turtle.test_Parseval = TurTLE.test.test_Parseval:main',
'turtle.test_fftw = TurTLE.test.test_fftw:main',
- 'turtle.test_Heun_p2p = TurTLE.test.test_Heun_p2p:main'],
+ 'turtle.test_Heun_p2p = TurTLE.test.test_Heun_p2p:main',
+ 'turtle.test_particle_deleter = TurTLE.test.test_particle_deleter:main'],
},
version = VERSION,
########################################################################
--
GitLab
From 5bc82b1c2ae9d0652f55a340ce8258516de3d4db Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sun, 28 Nov 2021 11:00:46 +0100
Subject: [PATCH 05/37] adds files needed for test
---
TurTLE/test/particle_set/particle_deleter.cpp | 213 ++++++++++
TurTLE/test/particle_set/particle_deleter.hpp | 98 +++++
TurTLE/test/test_particle_deleter.py | 364 ++++++++++++++++++
3 files changed, 675 insertions(+)
create mode 100644 TurTLE/test/particle_set/particle_deleter.cpp
create mode 100644 TurTLE/test/particle_set/particle_deleter.hpp
create mode 100644 TurTLE/test/test_particle_deleter.py
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
new file mode 100644
index 00000000..a4e1ebf3
--- /dev/null
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -0,0 +1,213 @@
+/**********************************************************************
+* *
+* Copyright 2021 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 *
+* *
+**********************************************************************/
+
+
+
+#include <string>
+#include <cmath>
+//#include "particle_deleter.hpp"
+#include "scope_timer.hpp"
+
+template <typename rnumber>
+int particle_deleter<rnumber>::initialize(void)
+{
+ TIMEZONE("particle_deleter::intialize");
+ this->NSVE<rnumber>::initialize();
+
+
+ // allocate previous velocity
+ this->previous_velocity = new field<rnumber, FFTW, THREE>(
+ this->fs->cvelocity->rlayout->sizes[2],
+ this->fs->cvelocity->rlayout->sizes[1],
+ this->fs->cvelocity->rlayout->sizes[0],
+ this->fs->cvelocity->rlayout->comm,
+ this->fs->cvelocity->fftw_plan_rigor);
+ this->fs->compute_velocity(this->fs->cvorticity);
+ *this->previous_velocity = *this->fs->cvelocity;
+ this->previous_velocity->ift();
+
+ this->fti = new field_tinterpolator<rnumber, FFTW, THREE, LINEAR>();
+ this->trhs = new tracer_with_collision_counter_rhs<rnumber, FFTW, LINEAR>();
+
+ // We're not using Adams-Bashforth in this code.
+ // This assert is there to ensure
+ // user knows what's happening.
+ assert(tracers0_integration_steps == 0);
+ // neighbours and smoothness are second and third template parameters of particle_set
+ assert(tracers0_neighbours == 3);
+ assert(tracers0_smoothness == 2);
+ this->pset = new particle_set<3, 3, 2>(
+ this->fs->cvelocity->rlayout,
+ this->dkx,
+ this->dky,
+ this->dkz,
+ 0.5);
+
+ // set two fields for the temporal interpolation
+ this->fti->set_field(this->previous_velocity, 0);
+ this->fti->set_field(this->fs->cvelocity, 1);
+ this->trhs->setVelocity(this->fti);
+
+ particles_input_hdf5<long long int, double, 3, 3> particle_reader(
+ this->comm,
+ this->fs->get_current_fname(),
+ std::string("/tracers0/state/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ std::string("/tracers0/rhs/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ this->pset->getSpatialLowLimitZ(),
+ this->pset->getSpatialUpLimitZ());
+ this->pset->init(particle_reader);
+
+ this->psolver = new particle_solver(*(this->pset), 0);
+
+ this->particles_output_writer_mpi = new particles_output_hdf5<
+ long long int, double, 3>(
+ MPI_COMM_WORLD,
+ "tracers0",
+ nparticles,
+ tracers0_integration_steps);
+ this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+ long long int, double, double, 3>(
+ MPI_COMM_WORLD,
+ this->pset->getTotalNumberOfParticles(),
+ (this->simname + "_particles.h5"),
+ "tracers0",
+ "position/0");
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int particle_deleter<rnumber>::step(void)
+{
+ TIMEZONE("particle_deleter::step");
+ // compute next step Navier-Stokes
+ this->NSVE<rnumber>::step();
+ // update velocity 1
+ this->fs->compute_velocity(this->fs->cvorticity);
+ this->fs->cvelocity->ift();
+ // compute particle step using velocity 0 and velocity 1
+ this->psolver->Heun(this->dt, *(this->trhs));
+ // update velocity 0
+ *this->previous_velocity = *this->fs->cvelocity;
+ this->psolver->setIteration(this->fs->iteration);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int particle_deleter<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("particle_deleter::write_checkpoint");
+ this->NSVE<rnumber>::write_checkpoint();
+ this->psolver->setIteration(this->fs->iteration);
+ this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+ this->psolver->template writeCheckpoint<3>(this->particles_output_writer_mpi);
+ this->particles_output_writer_mpi->close_file();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int particle_deleter<rnumber>::finalize(void)
+{
+ TIMEZONE("particle_deleter::finalize");
+ delete this->particles_output_writer_mpi;
+ delete this->particles_sample_writer_mpi;
+ delete this->psolver;
+ delete this->pset;
+ delete this->trhs;
+ delete this->fti;
+ this->NSVE<rnumber>::finalize();
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample fields at particle locations.
+ */
+
+template <typename rnumber>
+int particle_deleter<rnumber>::do_stats()
+{
+ TIMEZONE("particle_deleter::do_stats");
+ /// fluid stats go here
+ this->NSVE<rnumber>::do_stats();
+
+
+ /// either one of two conditions suffices to compute statistics:
+ /// 1) current iteration is a multiple of niter_part
+ /// 2) we are within niter_part_fine_duration/2 of a multiple of niter_part_fine_period
+ if (!(this->iteration % this->niter_part == 0 ||
+ ((this->iteration + this->niter_part_fine_duration/2) % this->niter_part_fine_period <=
+ this->niter_part_fine_duration)))
+ return EXIT_SUCCESS;
+
+ // sample position
+ this->pset->writeStateTriplet(
+ 0,
+ this->particles_sample_writer_mpi,
+ "tracers0",
+ "position",
+ psolver->getIteration());
+
+ // sample velocity
+ // first ensure velocity field is computed, and is in real space
+ if (!(this->iteration % this->niter_stat == 0))
+ {
+ // we need to compute velocity field manually, because it didn't happen in NSVE::do_stats()
+ this->fs->compute_velocity(this->fs->cvorticity);
+ *this->tmp_vec_field = this->fs->cvelocity->get_cdata();
+ this->tmp_vec_field->ift();
+ }
+ this->pset->writeSample(
+ this->tmp_vec_field,
+ this->particles_sample_writer_mpi,
+ "tracers0",
+ "velocity",
+ psolver->getIteration());
+
+ this->pset->writeParticleIndex(
+ "particle_index.h5",
+ "tracers0",
+ "index",
+ psolver->getIteration());
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int particle_deleter<rnumber>::read_parameters(void)
+{
+ TIMEZONE("particle_deleter::read_parameters");
+ this->NSVE<rnumber>::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+ this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+ this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+ this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+ this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
+ this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
+ this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
+ H5Fclose(parameter_file);
+ return EXIT_SUCCESS;
+}
+
+template class particle_deleter<float>;
+template class particle_deleter<double>;
+
diff --git a/TurTLE/test/particle_set/particle_deleter.hpp b/TurTLE/test/particle_set/particle_deleter.hpp
new file mode 100644
index 00000000..0561b56d
--- /dev/null
+++ b/TurTLE/test/particle_set/particle_deleter.hpp
@@ -0,0 +1,98 @@
+/**********************************************************************
+* *
+* Copyright 2021 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 PARTICLE_DELETER_HPP
+#define PARTICLE_DELETER_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/interpolation/field_tinterpolator.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "particles/particle_solver.hpp"
+#include "particles/rhs/tracer_with_collision_counter_rhs.hpp"
+
+/** \brief Test of particle deletion functionality.
+ *
+ * Child of Navier Stokes vorticity equation solver, this class calls all the
+ * methods from `NSVE`, and in addition:
+ * --- integrates simple Lagrangian tracers in the resulting velocity field.
+ * --- removes arbitrary particles from the set of tracers at each time step.
+ * It also executes a p2p kernel that doesn't affect the particle trajectories
+ * in any way.
+ *
+ */
+
+template <typename rnumber>
+class particle_deleter: public NSVE<rnumber>
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ int niter_part;
+ int niter_part_fine_period;
+ int niter_part_fine_duration;
+ int nparticles;
+ int tracers0_integration_steps;
+ int tracers0_neighbours;
+ int tracers0_smoothness;
+
+ /* other stuff */
+ field_tinterpolator<rnumber, FFTW, THREE, LINEAR> *fti;
+ tracer_with_collision_counter_rhs<rnumber, FFTW, LINEAR> *trhs;
+ particle_set<3, 3, 2> *pset;
+ particle_solver *psolver;
+
+ field<rnumber, FFTW, THREE> *previous_velocity;
+
+ particles_output_hdf5<long long int, double, 3> *particles_output_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, double, 3> *particles_sample_writer_mpi;
+
+ particle_deleter(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ ~particle_deleter(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//PARTICLE_DELETER_HPP
+
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
new file mode 100644
index 00000000..9e98271e
--- /dev/null
+++ b/TurTLE/test/test_particle_deleter.py
@@ -0,0 +1,364 @@
+import os
+import sys
+import argparse
+import getpass
+
+import numpy as np
+import h5py
+
+import TurTLE
+import TurTLE._base
+import TurTLE.DNS
+from TurTLE._code import _code
+
+cpp_location = os.path.join(
+ TurTLE.data_dir, 'particle_set')
+
+
+class ADNS(TurTLE.DNS):
+ def __init__(
+ self,
+ **kwargs):
+ TurTLE.DNS.__init__(self, **kwargs)
+ return None
+ def write_src(
+ self):
+ self.version_message = (
+ '/***********************************************************************\n' +
+ '* this code automatically generated by TurTLE\n' +
+ '* version {0}\n'.format(TurTLE.__version__) +
+ '***********************************************************************/\n\n\n')
+ self.include_list = [
+ '"base.hpp"',
+ '"scope_timer.hpp"',
+ '"fftw_interface.hpp"',
+ '"full_code/main_code.hpp"',
+ '"full_code/NSVEparticles.hpp"',
+ '<cmath>',
+ '<iostream>',
+ '<hdf5.h>',
+ '<string>',
+ '<cstring>',
+ '<fftw3-mpi.h>',
+ '<omp.h>',
+ '<cfenv>',
+ '<cstdlib>']
+ self.main = """
+ int main(int argc, char *argv[])
+ {{
+ bool fpe = (
+ (getenv("BFPS_FPE_OFF") == nullptr) ||
+ (getenv("BFPS_FPE_OFF") != std::string("TRUE")));
+ return main_code< {0} >(argc, argv, fpe);
+ }}
+ """.format(self.dns_type + '<{0}>'.format(self.C_field_dtype))
+ self.includes = '\n'.join(
+ ['#include ' + hh
+ for hh in self.include_list])
+ self.name = 'particle_deleter'
+ self.dns_type = 'particle_deleter'
+ self.definitions += open(
+ os.path.join(
+ cpp_location, self.dns_type + '.hpp'), 'r').read()
+ self.definitions += open(
+ os.path.join(
+ cpp_location, self.dns_type + '.cpp'), 'r').read()
+ with open(self.name + '.cpp', 'w') as outfile:
+ outfile.write(self.version_message + '\n\n')
+ outfile.write(self.includes + '\n\n')
+ outfile.write(self.definitions + '\n\n')
+ outfile.write(self.main + '\n')
+ self.check_current_vorticity_exists = True
+ return None
+ def generate_default_parameters(self):
+ # these parameters are relevant for all DNS classes
+ self.parameters['fftw_plan_rigor'] = 'FFTW_ESTIMATE'
+ self.parameters['dealias_type'] = int(1)
+ self.parameters['dkx'] = float(1.0)
+ self.parameters['dky'] = float(1.0)
+ self.parameters['dkz'] = float(1.0)
+ self.parameters['niter_todo'] = int(200)
+ self.parameters['niter_stat'] = int(20)
+ self.parameters['niter_out'] = int(200)
+ self.parameters['checkpoints_per_file'] = int(1)
+ self.parameters['dt'] = float(0.0001)
+ self.parameters['histogram_bins'] = int(64)
+ self.parameters['max_velocity_estimate'] = float(1)
+ self.parameters['max_vorticity_estimate'] = float(1)
+ self.parameters['niter_part'] = int(1)
+ self.parameters['niter_part_fine_period'] = int(2)
+ self.parameters['niter_part_fine_duration'] = int(0)
+ self.parameters['nparticles'] = int(1000)
+ self.parameters['tracers0_integration_steps'] = int(0)
+ self.parameters['tracers0_neighbours'] = int(3)
+ self.parameters['tracers0_smoothness'] = int(2)
+ self.parameters['tracers0_enable_p2p'] = int(0)
+ self.parameters['tracers0_enable_inner'] = int(0)
+ self.parameters['tracers0_enable_vorticity_omega'] = int(0)
+ self.parameters['tracers0_cutoff'] = float(0.920272)
+ self.parameters['tracers0_inner_v0'] = float(1)
+ self.parameters['tracers0_lambda'] = float(1)
+
+ self.parameters['nu'] = float(0.1)
+ self.parameters['fmode'] = int(1)
+ self.parameters['famplitude'] = float(0.5)
+ self.parameters['friction_coefficient'] = float(0.5)
+ self.parameters['injection_rate'] = float(0.4)
+ self.parameters['fk0'] = float(2.0)
+ self.parameters['fk1'] = float(4.0)
+ self.parameters['energy'] = float(0.5)
+ self.parameters['forcing_type'] = 'fixed_energy_injection_rate'
+ return None
+ def generate_tracer_state(
+ self,
+ rseed = None,
+ species = 0,
+ integration_steps = None,
+ ncomponents = 3):
+ try:
+ if type(integration_steps) == type(None):
+ integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps']
+ if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys():
+ integration_steps = self.parameters['tracers{0}_integration_steps'.format(species)]
+ with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
+ nn = self.parameters['nparticles']
+ if not 'tracers{0}'.format(species) in data_file.keys():
+ data_file.create_group('tracers{0}'.format(species))
+ data_file.create_group('tracers{0}/rhs'.format(species))
+ data_file.create_group('tracers{0}/state'.format(species))
+ data_file['tracers{0}/rhs'.format(species)].create_dataset(
+ '0',
+ shape = (integration_steps, nn, ncomponents,),
+ dtype = np.float)
+ dset = data_file['tracers{0}/state'.format(species)].create_dataset(
+ '0',
+ shape = (nn, ncomponents,),
+ dtype = np.float)
+ if not type(rseed) == type(None):
+ np.random.seed(rseed)
+ cc = int(0)
+ batch_size = int(1e6)
+ def get_random_phases(npoints):
+ return np.random.random(
+ (npoints, 3))*2*np.pi
+ def get_random_versors(npoints):
+ bla = np.random.normal(
+ size = (npoints, 3))
+ bla /= np.sum(bla**2, axis = 1)[:, None]**.5
+ return bla
+ while nn > 0:
+ #TODO introduce initial vel parameter as in established Stokes class, now set to zero
+ if nn > batch_size:
+ dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
+ if dset.shape[1] == 6:
+ dset[cc*batch_size:(cc+1)*batch_size, 3:] = 0.0* get_random_versors(batch_size)
+ nn -= batch_size
+ else:
+ dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
+ if dset.shape[1] == 6:
+ dset[cc*batch_size:cc*batch_size+nn, 3:] = 0.0* get_random_versors(nn)
+ nn = 0
+ cc += 1
+ except Exception as e:
+ print(e)
+ return None
+ def prepare_launch(
+ self,
+ args = [],
+ **kwargs):
+ parser = argparse.ArgumentParser('turtle ' + type(self).__name__)
+ self.job_parser_arguments(parser)
+ self.simulation_parser_arguments(parser)
+ self.particle_parser_arguments(parser)
+ self.parameters_to_parser_arguments(parser)
+ opt = parser.parse_args(args)
+ self.simname=opt.simname
+ self.set_precision(opt.precision)
+ self.dns_type = 'particle_deleter'
+ self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
+ if ((self.parameters['niter_todo'] % self.parameters['niter_out']) != 0):
+ self.parameters['niter_out'] = self.parameters['niter_todo']
+ if len(opt.src_work_dir) == 0:
+ opt.src_work_dir = os.path.realpath(opt.work_dir)
+ if type(opt.dkx) == type(None):
+ opt.dkx = 2. / opt.Lx
+ if type(opt.dky) == type(None):
+ opt.dky = 2. / opt.Ly
+ if type(opt.dkz) == type(None):
+ opt.dkz = 2. / opt.Lz
+ if type(opt.nx) == type(None):
+ opt.nx = opt.n
+ if type(opt.ny) == type(None):
+ opt.ny = opt.n
+ if type(opt.nz) == type(None):
+ opt.nz = opt.n
+ if type(opt.fk0) == type(None):
+ opt.fk0 = self.parameters['fk0']
+ if type(opt.fk1) == type(None):
+ opt.fk1 = self.parameters['fk1']
+ if type(opt.injection_rate) == type(None):
+ opt.injection_rate = self.parameters['injection_rate']
+ if type(opt.dealias_type) == type(None):
+ opt.dealias_type = self.parameters['dealias_type']
+ if (opt.nx > opt.n or
+ opt.ny > opt.n or
+ opt.nz > opt.n):
+ opt.n = min(opt.nx, opt.ny, opt.nz)
+ print("Warning: '-n' parameter changed to minimum of nx, ny, nz. This affects the computation of nu.")
+ self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
+ # check value of kMax
+ kM = opt.n * 0.5
+ if opt.dealias_type == 1:
+ kM *= 0.8
+ # tweak forcing/viscosity based on forcint type
+ if opt.forcing_type == 'linear':
+ # custom famplitude for 288 and 576
+ if opt.n == 288:
+ self.parameters['famplitude'] = 0.45
+ elif opt.n == 576:
+ self.parameters['famplitude'] = 0.47
+ elif opt.forcing_type == 'fixed_energy_injection_rate':
+ # use the fact that mean dissipation rate is equal to injection rate
+ self.parameters['nu'] = (
+ opt.injection_rate *
+ (opt.kMeta / kM)**4)**(1./3)
+ elif opt.forcing_type == 'fixed_energy':
+ kf = 1. / (1./opt.fk0 +
+ 1./opt.fk1)
+ self.parameters['nu'] = (
+ (opt.kMeta / kM)**(4./3) *
+ (np.pi / kf)**(1./3) *
+ (2*self.parameters['energy'] / 3)**0.5)
+ if type(opt.checkpoints_per_file) == type(None):
+ # hardcoded FFTW complex representation size
+ field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
+ checkpoint_size = field_size
+ self.set_host_info(TurTLE.host_info)
+ if type(opt.environment) != type(None):
+ self.host_info['environment'] = opt.environment
+ self.pars_from_namespace(opt)
+ return opt
+
+ def write_par(
+ self,
+ iter0 = 0):
+ assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
+ assert (self.parameters['niter_todo'] % self.parameters['niter_out'] == 0)
+ assert (self.parameters['niter_out'] % self.parameters['niter_stat'] == 0)
+ assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
+ assert (self.parameters['niter_out'] % self.parameters['niter_part'] == 0)
+ _code.write_par(self, iter0 = iter0)
+ with h5py.File(self.get_data_file_name(), 'r+') as ofile:
+ ofile['code_info/exec_name'] = self.name
+ kspace = self.get_kspace()
+ for k in kspace.keys():
+ ofile['kspace/' + k] = kspace[k]
+ nshells = kspace['nshell'].shape[0]
+ kspace = self.get_kspace()
+ nshells = kspace['nshell'].shape[0]
+ vec_stat_datasets = ['velocity', 'vorticity']
+ scal_stat_datasets = []
+ for k in vec_stat_datasets:
+ time_chunk = 2**20//(8*3*3*nshells)
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+ (1, nshells, 3, 3),
+ chunks = (time_chunk, nshells, 3, 3),
+ maxshape = (None, nshells, 3, 3),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*4*10)
+ time_chunk = max(time_chunk, 1)
+ a = ofile.create_dataset('statistics/moments/' + k,
+ (1, 10, 4),
+ chunks = (time_chunk, 10, 4),
+ maxshape = (None, 10, 4),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/histograms/' + k,
+ (1,
+ self.parameters['histogram_bins'],
+ 4),
+ chunks = (time_chunk,
+ self.parameters['histogram_bins'],
+ 4),
+ maxshape = (None,
+ self.parameters['histogram_bins'],
+ 4),
+ dtype = np.int64)
+ #TODO move collision stats to particle stats
+ for i in range(self.parameters['niter_out']):
+ ofile.create_dataset('statistics/collisions/'+str(i),
+ shape=(1,),
+ dtype = np.int64)
+ ofile['checkpoint'] = int(0)
+ if (self.dns_type in ['NSVE', 'NSVE_no_output']):
+ return None
+ return None
+
+ def launch(
+ self,
+ args = [],
+ **kwargs):
+ opt = self.prepare_launch(args = args)
+ if not os.path.exists(self.get_data_file_name()):
+ self.write_par()
+ self.generate_initial_condition(opt = opt)
+ with h5py.File(self.get_particle_file_name(), 'w') as particle_file:
+ particle_file.create_group('tracers0/position')
+ particle_file.create_group('tracers0/velocity')
+ particle_file.create_group('collisions0')
+ with h5py.File('particle_index.h5', 'w') as particle_file:
+ particle_file.create_group('tracers0/index')
+ self.generate_tracer_state(integration_steps = self.parameters['tracers0_integration_steps'],
+ rseed = opt.particle_rand_seed,
+ ncomponents = 3)
+ self.run(
+ nb_processes = opt.nb_processes,
+ nb_threads_per_process = opt.nb_threads_per_process,
+ njobs = opt.njobs,
+ hours = opt.minutes // 60,
+ minutes = opt.minutes % 60,
+ no_submit = opt.no_submit,
+ no_debug = opt.no_debug)
+ return None
+
+def main():
+ bla = ADNS()
+ bla.launch(sys.argv[1:])
+ assert(bla.parameters['niter_part'] == 1)
+ df = bla.get_particle_file()
+ x = []
+ v = []
+ for ii in range(bla.parameters['niter_out']):
+ x.append(df['tracers0/position/{0}'.format(ii)][()])
+ v.append(df['tracers0/velocity/{0}'.format(ii)][()])
+ df.close()
+ x = np.array(x)
+ v = np.array(v)
+ d1x = (x[2:] - x[:-2]) / (2*bla.parameters['dt'])
+ vv = v[1:-1]
+ diff = vv - d1x
+ relative_error = (np.sum(diff**2, axis = 2) / np.sum(vv**2, axis = 2))**0.5
+ tindex, trajindex = np.unravel_index(
+ np.argmax(relative_error),
+ relative_error.shape)
+ try:
+ import matplotlib.pyplot as plt
+ f = plt.figure()
+ a = f.add_subplot(211)
+ a.plot(vv[:, trajindex])
+ a.plot(d1x[:, trajindex], dashes = (2, 2))
+ a = f.add_subplot(212)
+ a.plot(x[:, trajindex])
+ a.set_title('maximum error for (tindex, trajindex) = ({0}, {1})'.format(tindex, trajindex))
+ f.tight_layout()
+ f.savefig('traj.pdf')
+ except ImportError:
+ print('couldn\'t import matplotlib, no figure')
+ print(np.max(relative_error))
+ assert(np.max(relative_error) < 1e-4)
+ return None
+
+if __name__ == '__main__':
+ main()
--
GitLab
From 210fc9a70bcfbac1fce5ebdfa543957fa0909832 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sun, 28 Nov 2021 16:58:57 +0100
Subject: [PATCH 06/37] [wip] still debugging particle_delete functionality
---
TurTLE/test/particle_set/particle_deleter.cpp | 26 ++++
.../interpolation/abstract_particle_set.hpp | 3 +
cpp/particles/interpolation/particle_set.hpp | 123 ++++++++++++++++++
3 files changed, 152 insertions(+)
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
index a4e1ebf3..af7d93ca 100644
--- a/TurTLE/test/particle_set/particle_deleter.cpp
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -188,6 +188,32 @@ int particle_deleter<rnumber>::do_stats()
"index",
psolver->getIteration());
+ std::vector<long long int> indices_to_delete(10);
+ indices_to_delete[0] = 1;
+ indices_to_delete[1] = 3;
+ indices_to_delete[2] = 19;
+ indices_to_delete[3] = 13;
+ indices_to_delete[4] = 101;
+ indices_to_delete[5] = 91;
+ indices_to_delete[6] = 141;
+ indices_to_delete[7] = 143;
+ indices_to_delete[8] = 145;
+ indices_to_delete[9] = 149;
+
+ DEBUG_MSG("before delete particles");
+ this->pset->_print_debug_info();
+ this->pset->deleteParticles(indices_to_delete);
+ DEBUG_MSG("after delete particles");
+ this->pset->_print_debug_info();
+
+ DEBUG_MSG("before writeParticleIndex");
+ this->pset->writeParticleIndex(
+ "particle_index.h5",
+ "tracers0",
+ "index",
+ psolver->getIteration()+1);
+ DEBUG_MSG("after writeParticleIndex");
+
return EXIT_SUCCESS;
}
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index eea52c17..c0a1128e 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -188,6 +188,9 @@ class abstract_particle_set
return EXIT_SUCCESS;
}
+ virtual int deleteParticles(
+ const std::vector<partsize_t> &inIndicesToDelete) = 0;
+
int writeParticleIndex(
const std::string file_name,
const std::string species_name,
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 9abb2c1b..5d41c1c6 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -29,6 +29,7 @@
#include <array>
#include <cassert>
#include <vector>
+#include <set>
#include <memory>
#include "field_layout.hpp"
#include "field.hpp"
@@ -167,6 +168,58 @@ class particle_set: public abstract_particle_set
this->offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
}
+ int _print_debug_info()
+ {
+ int myrank;
+ AssertMpi(MPI_Comm_rank(mpi_comm, &myrank));
+ int nprocs;
+ AssertMpi(MPI_Comm_size(mpi_comm, &nprocs));
+
+ for (int rr = 0; rr < nprocs; rr++)
+ {
+ if (myrank == rr)
+ {
+ DEBUG_MSG("##########################################\n");
+ DEBUG_MSG("# particle_set debug info follows\n");
+
+ DEBUG_MSG("current partition interval %d %d\n",
+ this->current_partition_interval.first,
+ this->current_partition_interval.second);
+ DEBUG_MSG("partition_interval_size %d\n",
+ this->partition_interval_size);
+
+ for (int i=0; i<this->partition_interval_size; i++)
+ {
+ DEBUG_MSG("number_particles_per_partition[%d] = %d\n",
+ i,
+ number_particles_per_partition[i]);
+ }
+
+ for (int i=0; i<this->partition_interval_size+1; i++)
+ {
+ DEBUG_MSG("offset_particles_for_partition[%d] = %d\n",
+ i,
+ offset_particles_for_partition[i]);
+ }
+
+ DEBUG_MSG("local_number_of_particles = %d\n",
+ this->local_number_of_particles);
+ DEBUG_MSG("total_number_of_particles = %d\n",
+ this->total_number_of_particles);
+
+ for (int i=0; i<this->particle_file_layout.size(); i++)
+ DEBUG_MSG("particle_file_layout[%d] = %d\n",
+ i,
+ this->particle_file_layout[i]);
+
+ DEBUG_MSG("# particle_set debug info ends\n");
+ DEBUG_MSG("##########################################\n");
+ }
+ MPI_Barrier(this->mpi_comm);
+ }
+ return EXIT_SUCCESS;
+ }
+
~particle_set() noexcept(false){}
particle_rnumber* getParticleState() const
@@ -431,6 +484,76 @@ class particle_set: public abstract_particle_set
{
return &(this->p2pComputer);
}
+
+ std::vector<partsize_t> selectLocalParticleFromFewIndices(
+ const std::vector<partsize_t> &inGlobalIndicesToDelete)
+ {
+ std::vector<partsize_t> local_indices;
+ for (partsize_t i=0; i<this->getLocalNumberOfParticles(); i++)
+ for (partsize_t ii=0; ii<inGlobalIndicesToDelete.size(); ii++)
+ {
+ if (this->getParticleIndex()[i] == inGlobalIndicesToDelete[ii])
+ local_indices.push_back(inGlobalIndicesToDelete[ii]);
+ }
+ return local_indices;
+ }
+
+ /** \brief Delete particles from arbitrary list
+ *
+ * \note List of particles given to this method may contain indices to
+ * particles at arbitrary locations.
+ */
+ int deleteParticles(
+ const std::vector<partsize_t> &inGlobalIndicesToDelete)
+ {
+ std::vector<partsize_t> inLocalIndicesToDelete =
+ this->selectLocalParticleFromFewIndices(inGlobalIndicesToDelete);
+ std::set<partsize_t> setIndicesToDelete(inLocalIndicesToDelete.begin(), inLocalIndicesToDelete.end());
+ if (setIndicesToDelete.size() > 0)
+ {
+ partsize_t number_of_deleted_particles = 0;
+
+ for(int idxPartition = 0 ; idxPartition < partition_interval_size ; ++idxPartition){
+ const partsize_t nbDeletedParticlesBefore = number_of_deleted_particles;
+ for (partsize_t idx = this->offset_particles_for_partition[idxPartition];
+ idx < this->offset_particles_for_partition[idxPartition+1];
+ ++idx)
+ {
+ if(setIndicesToDelete.find(this->local_index[idx]) != setIndicesToDelete.end())
+ {
+ number_of_deleted_particles += 1;
+ }
+ else if(number_of_deleted_particles)
+ {
+ this->local_index[idx-number_of_deleted_particles] = this->local_index[idx];
+ for(int idx_pos = 0 ; idx_pos < state_size; ++idx_pos)
+ {
+ this->local_state[(idx-number_of_deleted_particles)*state_size+idx_pos] =
+ this->local_state[idx*state_size+idx_pos];
+ }
+ }
+ }
+
+ this->offset_particles_for_partition[idxPartition] -= nbDeletedParticlesBefore;
+ }
+
+ this->offset_particles_for_partition[partition_interval_size] -= number_of_deleted_particles;
+
+ this->local_number_of_particles -= number_of_deleted_particles;
+ assert(this->local_number_of_particles == this->offset_particles_for_partition[partition_interval_size]);
+ }
+
+ AssertMpi(MPI_Allreduce(
+ const_cast<partsize_t*>(&(this->local_number_of_particles)),
+ &(this->total_number_of_particles),
+ 1,
+ particles_utils::GetMpiType(partsize_t()),
+ MPI_SUM,
+ this->mpi_comm));
+ this->particle_file_layout.resize(1);
+ this->particle_file_layout[0] = this->total_number_of_particles;
+ return EXIT_SUCCESS;
+ }
};
#endif//PARTICLE_SET_HPP
--
GitLab
From 5a37f9826b9c2dc58130effb543d6e381afb61e2 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Mon, 29 Nov 2021 11:49:59 +0100
Subject: [PATCH 07/37] [wip] adds debug message
clearly, the issue is that particle labels are used to place data in the
output buffer. So there's an assumption about labels being smaller than
the size of the buffer.
I still have to think of a solution to this.
---
cpp/particles/abstract_particles_output.hpp | 2 ++
1 file changed, 2 insertions(+)
diff --git a/cpp/particles/abstract_particles_output.hpp b/cpp/particles/abstract_particles_output.hpp
index de926680..cadc65ce 100644
--- a/cpp/particles/abstract_particles_output.hpp
+++ b/cpp/particles/abstract_particles_output.hpp
@@ -311,6 +311,8 @@ public:
const partsize_t src_idx = idx_part;
const partsize_t dst_idx = buffer_indexes_recv[idx_part]-particles_chunk_current_offset;
assert(0 <= dst_idx);
+ DEBUG_MSG("dst_idx is %d, particles_chunk_current_size is %d\n",
+ dst_idx, particles_chunk_current_size);
assert(dst_idx < particles_chunk_current_size);
for(int idx_val = 0 ; idx_val < size_particle_positions ; ++idx_val){
--
GitLab
From cc2912b4d0d37726642dee4d88fd2de8b544bab0 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 20:13:37 +0100
Subject: [PATCH 08/37] mark destructors with noexcept(false)
---
cpp/field.hpp | 2 +-
cpp/field_binary_IO.hpp | 2 +-
cpp/field_layout.hpp | 4 ++--
cpp/kspace.hpp | 2 +-
cpp/scope_timer.hpp | 2 +-
cpp/spectrum_function.hpp | 2 +-
cpp/vorticity_equation.hpp | 2 +-
7 files changed, 8 insertions(+), 8 deletions(-)
diff --git a/cpp/field.hpp b/cpp/field.hpp
index a42d4801..8da33703 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -87,7 +87,7 @@ class field
const int nz,
const MPI_Comm COMM_TO_USE,
const unsigned FFTW_PLAN_RIGOR = DEFAULT_FFTW_FLAG);
- ~field();
+ ~field() noexcept(false);
int io(
const std::string fname,
diff --git a/cpp/field_binary_IO.hpp b/cpp/field_binary_IO.hpp
index 66234ab3..e1c839cc 100644
--- a/cpp/field_binary_IO.hpp
+++ b/cpp/field_binary_IO.hpp
@@ -64,7 +64,7 @@ class field_binary_IO:public field_layout<fc>
const hsize_t *SUBSIZES,
const hsize_t *STARTS,
const MPI_Comm COMM_TO_USE);
- ~field_binary_IO();
+ ~field_binary_IO() noexcept(false);
int read(
const std::string fname,
diff --git a/cpp/field_layout.hpp b/cpp/field_layout.hpp
index edc094a9..27d0ea9a 100644
--- a/cpp/field_layout.hpp
+++ b/cpp/field_layout.hpp
@@ -52,7 +52,7 @@ constexpr unsigned int ndim(
class abstract_field_layout
{
public:
- virtual ~abstract_field_layout(){}
+ virtual ~abstract_field_layout() noexcept(false){}
virtual MPI_Comm getMPIComm() const = 0;
virtual int getMPIRank() const = 0;
@@ -85,7 +85,7 @@ class field_layout: public abstract_field_layout
const hsize_t *SUBSIZES,
const hsize_t *STARTS,
const MPI_Comm COMM_TO_USE);
- ~field_layout(){}
+ ~field_layout() noexcept(false){}
MPI_Comm getMPIComm() const
{
diff --git a/cpp/kspace.hpp b/cpp/kspace.hpp
index f082f11e..912129f6 100644
--- a/cpp/kspace.hpp
+++ b/cpp/kspace.hpp
@@ -77,7 +77,7 @@ class kspace
const double DKX = 1.0,
const double DKY = 1.0,
const double DKZ = 1.0);
- ~kspace();
+ ~kspace() noexcept(false);
int store(hid_t stat_file);
diff --git a/cpp/scope_timer.hpp b/cpp/scope_timer.hpp
index 9e822e65..51fc2d38 100644
--- a/cpp/scope_timer.hpp
+++ b/cpp/scope_timer.hpp
@@ -243,7 +243,7 @@ public:
omp_init_lock(&m_recordsLock);
}
- ~EventManager() throw() {
+ ~EventManager() noexcept(false) {
assert(m_currentEventsStackPerThread[0].size() == 1);
assert(m_currentEventsStackPerThread[0].top().size() == 1);
diff --git a/cpp/spectrum_function.hpp b/cpp/spectrum_function.hpp
index fcc729ff..2df2ecab 100644
--- a/cpp/spectrum_function.hpp
+++ b/cpp/spectrum_function.hpp
@@ -20,7 +20,7 @@ class spectrum_function
{
assert(this->values.size() == this->kk->nshells);
}
- ~spectrum_function(){}
+ ~spectrum_function() noexcept(false){}
double operator()(double kvalue)
{
diff --git a/cpp/vorticity_equation.hpp b/cpp/vorticity_equation.hpp
index 8bf2ea8e..2ba11293 100644
--- a/cpp/vorticity_equation.hpp
+++ b/cpp/vorticity_equation.hpp
@@ -87,7 +87,7 @@ class vorticity_equation
double DKY = 1.0,
double DKZ = 1.0,
unsigned FFTW_PLAN_RIGOR = FFTW_MEASURE);
- virtual ~vorticity_equation(void);
+ virtual ~vorticity_equation(void) noexcept(false);
/* solver essential methods */
virtual void omega_nonlin(int src);
--
GitLab
From 95278bcf4cfac16b8d1b735c026883300aa611f2 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 20:17:23 +0100
Subject: [PATCH 09/37] marks more destructors as noexcept(false)
---
cpp/full_code/NSVE.hpp | 2 +-
cpp/full_code/NSVE_field_stats.hpp | 2 +-
cpp/full_code/NSVE_no_output.hpp | 2 +-
cpp/full_code/dealias_test.hpp | 2 +-
cpp/full_code/field_single_to_double.hpp | 2 +-
cpp/full_code/native_binary_to_hdf5.hpp | 2 +-
cpp/full_code/phase_shift_test.hpp | 2 +-
cpp/full_code/postprocess.hpp | 2 +-
cpp/full_code/pressure_stats.hpp | 2 +-
cpp/full_code/symmetrize_test.hpp | 2 +-
cpp/full_code/write_rpressure.hpp | 2 +-
11 files changed, 11 insertions(+), 11 deletions(-)
diff --git a/cpp/full_code/NSVE.hpp b/cpp/full_code/NSVE.hpp
index df8ffc5d..2f29f91f 100644
--- a/cpp/full_code/NSVE.hpp
+++ b/cpp/full_code/NSVE.hpp
@@ -75,7 +75,7 @@ class NSVE: public direct_numerical_simulation
direct_numerical_simulation(
COMMUNICATOR,
simulation_name){}
- ~NSVE(){}
+ ~NSVE() noexcept(false){}
int initialize(void);
int step(void);
diff --git a/cpp/full_code/NSVE_field_stats.hpp b/cpp/full_code/NSVE_field_stats.hpp
index de2078ee..a0421538 100644
--- a/cpp/full_code/NSVE_field_stats.hpp
+++ b/cpp/full_code/NSVE_field_stats.hpp
@@ -52,7 +52,7 @@ class NSVE_field_stats: public postprocess
postprocess(
COMMUNICATOR,
simulation_name){}
- virtual ~NSVE_field_stats(){}
+ virtual ~NSVE_field_stats() noexcept(false){}
virtual int initialize(void);
virtual int work_on_current_iteration(void);
diff --git a/cpp/full_code/NSVE_no_output.hpp b/cpp/full_code/NSVE_no_output.hpp
index d912aec4..d739745c 100644
--- a/cpp/full_code/NSVE_no_output.hpp
+++ b/cpp/full_code/NSVE_no_output.hpp
@@ -39,7 +39,7 @@ class NSVE_no_output: public NSVE<rnumber>
NSVE<rnumber>(
COMMUNICATOR,
simulation_name){}
- ~NSVE_no_output(){}
+ ~NSVE_no_output() noexcept(false){}
int write_checkpoint(void)
{
TIMEZONE("NSVE_no_output::write_checkpoint");
diff --git a/cpp/full_code/dealias_test.hpp b/cpp/full_code/dealias_test.hpp
index 17578928..ba765c45 100644
--- a/cpp/full_code/dealias_test.hpp
+++ b/cpp/full_code/dealias_test.hpp
@@ -65,7 +65,7 @@ class dealias_test: public test
test(
COMMUNICATOR,
simulation_name){}
- ~dealias_test(){}
+ ~dealias_test() noexcept(false){}
int initialize(void);
int do_work(void);
diff --git a/cpp/full_code/field_single_to_double.hpp b/cpp/full_code/field_single_to_double.hpp
index 04f56357..de98e63e 100644
--- a/cpp/full_code/field_single_to_double.hpp
+++ b/cpp/full_code/field_single_to_double.hpp
@@ -52,7 +52,7 @@ class field_single_to_double: public NSVE_field_stats<rnumber>
NSVE_field_stats<rnumber>(
COMMUNICATOR,
simulation_name){}
- virtual ~field_single_to_double(){}
+ virtual ~field_single_to_double() noexcept(false){}
int initialize(void);
int work_on_current_iteration(void);
diff --git a/cpp/full_code/native_binary_to_hdf5.hpp b/cpp/full_code/native_binary_to_hdf5.hpp
index d8ded73b..3baaf609 100644
--- a/cpp/full_code/native_binary_to_hdf5.hpp
+++ b/cpp/full_code/native_binary_to_hdf5.hpp
@@ -50,7 +50,7 @@ class native_binary_to_hdf5: public postprocess
postprocess(
COMMUNICATOR,
simulation_name){}
- virtual ~native_binary_to_hdf5(){}
+ virtual ~native_binary_to_hdf5() noexcept(false){}
int initialize(void);
int work_on_current_iteration(void);
diff --git a/cpp/full_code/phase_shift_test.hpp b/cpp/full_code/phase_shift_test.hpp
index 8e805300..e98f76cd 100644
--- a/cpp/full_code/phase_shift_test.hpp
+++ b/cpp/full_code/phase_shift_test.hpp
@@ -51,7 +51,7 @@ class phase_shift_test: public test
test(
COMMUNICATOR,
simulation_name){}
- ~phase_shift_test(){}
+ ~phase_shift_test() noexcept(false){}
int initialize(void);
int do_work(void);
diff --git a/cpp/full_code/postprocess.hpp b/cpp/full_code/postprocess.hpp
index c232a6b3..da22c91e 100644
--- a/cpp/full_code/postprocess.hpp
+++ b/cpp/full_code/postprocess.hpp
@@ -63,7 +63,7 @@ class postprocess: public code_base
code_base(
COMMUNICATOR,
simulation_name){}
- virtual ~postprocess(){}
+ virtual ~postprocess() noexcept(false){}
virtual int initialize(void) = 0;
virtual int work_on_current_iteration(void) = 0;
diff --git a/cpp/full_code/pressure_stats.hpp b/cpp/full_code/pressure_stats.hpp
index d8a54494..f20ba08e 100644
--- a/cpp/full_code/pressure_stats.hpp
+++ b/cpp/full_code/pressure_stats.hpp
@@ -64,7 +64,7 @@ class pressure_stats: public NSVE_field_stats<rnumber>
NSVE_field_stats<rnumber>(
COMMUNICATOR,
simulation_name){}
- virtual ~pressure_stats(){}
+ virtual ~pressure_stats() noexcept(false){}
int initialize(void);
int work_on_current_iteration(void);
diff --git a/cpp/full_code/symmetrize_test.hpp b/cpp/full_code/symmetrize_test.hpp
index 52fbac28..69443d1a 100644
--- a/cpp/full_code/symmetrize_test.hpp
+++ b/cpp/full_code/symmetrize_test.hpp
@@ -51,7 +51,7 @@ class symmetrize_test: public test
test(
COMMUNICATOR,
simulation_name){}
- ~symmetrize_test(){}
+ ~symmetrize_test() noexcept(false){}
int initialize(void);
int do_work(void);
diff --git a/cpp/full_code/write_rpressure.hpp b/cpp/full_code/write_rpressure.hpp
index 64fbcfc7..2a14dcb6 100644
--- a/cpp/full_code/write_rpressure.hpp
+++ b/cpp/full_code/write_rpressure.hpp
@@ -56,7 +56,7 @@ class write_rpressure: public NSVE_field_stats<rnumber>
NSVE_field_stats<rnumber>(
COMMUNICATOR,
simulation_name){}
- virtual ~write_rpressure(){}
+ virtual ~write_rpressure() noexcept(false){}
int initialize(void);
int work_on_current_iteration(void);
--
GitLab
From 5f76a5c1dff7031e35cfdc1677d94fcd8291ac9c Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 20:19:00 +0100
Subject: [PATCH 10/37] marks more destructors with noexcept(false)
---
cpp/particles/particle_solver.hpp | 2 +-
cpp/particles/rhs/deformation_tensor_rhs.hpp | 2 +-
cpp/particles/rhs/tracer_rhs.hpp | 2 +-
cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp | 2 +-
4 files changed, 4 insertions(+), 4 deletions(-)
diff --git a/cpp/particles/particle_solver.hpp b/cpp/particles/particle_solver.hpp
index e389add1..87443ed4 100644
--- a/cpp/particles/particle_solver.hpp
+++ b/cpp/particles/particle_solver.hpp
@@ -74,7 +74,7 @@ class particle_solver
this->additional_states.resize(number_of_additional_states);
this->particle_species_name = "tracers0";
}
- ~particle_solver(){}
+ ~particle_solver() noexcept(false){}
int setIteration(const int i)
{
diff --git a/cpp/particles/rhs/deformation_tensor_rhs.hpp b/cpp/particles/rhs/deformation_tensor_rhs.hpp
index cbf43cbc..0835f986 100644
--- a/cpp/particles/rhs/deformation_tensor_rhs.hpp
+++ b/cpp/particles/rhs/deformation_tensor_rhs.hpp
@@ -46,7 +46,7 @@ class deformation_tensor_rhs: public abstract_particle_rhs
public:
deformation_tensor_rhs(){}
- ~deformation_tensor_rhs(){}
+ ~deformation_tensor_rhs() noexcept(false){}
const int getStateSize() const
{
diff --git a/cpp/particles/rhs/tracer_rhs.hpp b/cpp/particles/rhs/tracer_rhs.hpp
index 03bfab35..680b46fd 100644
--- a/cpp/particles/rhs/tracer_rhs.hpp
+++ b/cpp/particles/rhs/tracer_rhs.hpp
@@ -41,7 +41,7 @@ class tracer_rhs: public abstract_particle_rhs
public:
tracer_rhs(){}
- ~tracer_rhs(){}
+ ~tracer_rhs() noexcept(false){}
const int getStateSize() const
{
diff --git a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
index 9ba41c30..5aa15213 100644
--- a/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
+++ b/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
@@ -42,7 +42,7 @@ class tracer_with_collision_counter_rhs: public tracer_rhs<rnumber, be, tt>
public:
tracer_with_collision_counter_rhs(){}
- ~tracer_with_collision_counter_rhs(){}
+ ~tracer_with_collision_counter_rhs() noexcept(false){}
virtual const int getStateSize() const
{
--
GitLab
From 76ef8450c0c54defa7ace020d0dcd1bcef378d14 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 20:20:06 +0100
Subject: [PATCH 11/37] [wip] adds comment on desired openmp link change
---
CMakeLists.txt | 13 ++++++++-----
1 file changed, 8 insertions(+), 5 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 96b3530c..3dc4417c 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -173,10 +173,6 @@ set(CMAKE_CXX_STANDARD_REQUIRED ON)
find_package(OpenMP REQUIRED)
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
-list(APPEND TURTLE_LIBS "${OpenMP_CXX_LIB_NAMES}")
-
#####################################################################################
## HDF5
@@ -234,7 +230,7 @@ endif()
#####################################################################################
## Extra flags
-set(CMAKE_CXX_COMPILE_FLAGS "${CMAKE_CXX_COMPILE_FLAGS} $ENV{TURTLE_COMPILATION_FLAGS} -Wall -g -fopenmp")
+set(CMAKE_CXX_COMPILE_FLAGS "${CMAKE_CXX_COMPILE_FLAGS} $ENV{TURTLE_COMPILATION_FLAGS} -Wall -g -fsanitize=address")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CMAKE_CXX_COMPILE_FLAGS}")
#####################################################################################
@@ -468,6 +464,13 @@ LIST(APPEND source_files ${hpp_for_lib} ${cpp_for_lib})
add_library(TurTLE ${source_files})
+# TODO: we should be using the bit below, but we need to fix everything else about TURTLE_LIBS before we can do that
+#target_link_libraries(TurTLE PUBLIC OpenMP::OpenMP_CXX)
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")
+list(APPEND TURTLE_LIBS "${OpenMP_CXX_LIB_NAMES}")
+
target_link_libraries(TurTLE ${TURTLE_LIBS})
install(TARGETS TurTLE EXPORT TURTLE_EXPORT DESTINATION lib/ )
--
GitLab
From 613ab51e80d90301844fcda121575b4adbfedb09 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 21:23:52 +0100
Subject: [PATCH 12/37] gets rid of warnings
---
cpp/field.cpp | 6 ++----
1 file changed, 2 insertions(+), 4 deletions(-)
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 3efc8ee2..cead918e 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -103,7 +103,7 @@ field<rnumber, be, fc>::field(
sizes, subsizes, starts, this->comm);
this->data = fftw_interface<rnumber>::alloc_real(
this->rmemlayout->local_size);
- char *plan_information = NULL;
+ //char *plan_information = NULL;
this->c2r_plan = fftw_interface<rnumber>::mpi_plan_many_dft_c2r(
3, nfftw, ncomp(fc),
FFTW_MPI_DEFAULT_BLOCK, FFTW_MPI_DEFAULT_BLOCK,
@@ -1812,8 +1812,6 @@ void field<rnumber, be, fc>::symmetrize_alternate()
ptrdiff_t cindex1 = this->get_cindex(0, 0, this->clayout->sizes[1] - iz);
for (int cc = 0; cc < int(ncomp(fc)); cc++)
{
- double re = ((*(cdata + cc + ncomp(fc)*cindex0))[0] + (*(cdata + cc + ncomp(fc)*cindex1))[0])/2;
- double im = ((*(cdata + cc + ncomp(fc)*cindex0))[1] - (*(cdata + cc + ncomp(fc)*cindex1))[1])/2;
const double ampp = sqrt(
(*(cdata + cc + ncomp(fc)*cindex0))[0]*(*(cdata + cc + ncomp(fc)*cindex0))[0] +
(*(cdata + cc + ncomp(fc)*cindex0))[1]*(*(cdata + cc + ncomp(fc)*cindex0))[1]);
@@ -3085,7 +3083,7 @@ int generate_random_phase_field(
}
}
}
- if (xindex == output_field->clayout->sizes[2]-1)
+ if (xindex == ptrdiff_t(output_field->clayout->sizes[2])-ptrdiff_t(1))
{
output_field->cval(cindex, 0) = 0;
output_field->cval(cindex, 1) = 0;
--
GitLab
From 441624c8bd0146c22c1a4b3ccc07c3ccb5f8a7a7 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 21:26:07 +0100
Subject: [PATCH 13/37] marks destructors as noexcept(false)
---
cpp/particles/interpolation/field_tinterpolator.hpp | 2 +-
cpp/particles/p2p/p2p_distr_mpi.hpp | 2 +-
2 files changed, 2 insertions(+), 2 deletions(-)
diff --git a/cpp/particles/interpolation/field_tinterpolator.hpp b/cpp/particles/interpolation/field_tinterpolator.hpp
index 22b830f8..c788f7c5 100644
--- a/cpp/particles/interpolation/field_tinterpolator.hpp
+++ b/cpp/particles/interpolation/field_tinterpolator.hpp
@@ -56,7 +56,7 @@ class field_tinterpolator
this->field_list[3] = NULL;
}
- ~field_tinterpolator(){}
+ ~field_tinterpolator() noexcept(false){}
int set_field(
field<rnumber, be, fc> *field_src = NULL,
diff --git a/cpp/particles/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index ecf6aa4d..b8fa7ec2 100644
--- a/cpp/particles/p2p/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p/p2p_distr_mpi.hpp
@@ -197,7 +197,7 @@ public:
nb_cell_levels[IDXC_Z] = nb_cells_factor;
}
- virtual ~p2p_distr_mpi(){}
+ virtual ~p2p_distr_mpi() noexcept(false){}
////////////////////////////////////////////////////////////////////////////
--
GitLab
From 9d207049f6241e6a1c7e54972c412d0de186b505 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 21:28:40 +0100
Subject: [PATCH 14/37] adds particle labels and indices instead of indexes
---
cpp/full_code/NSVEparticles.cpp | 2 +-
cpp/particles/abstract_particles_input.hpp | 1 +
cpp/particles/abstract_particles_output.hpp | 29 +-
.../interpolation/abstract_particle_set.hpp | 26 +-
cpp/particles/interpolation/particle_set.hpp | 313 +++++++++++++-----
cpp/particles/particle_solver.cpp | 2 +-
cpp/particles/particles_input_hdf5.hpp | 12 +
cpp/particles/particles_input_random.hpp | 13 +
8 files changed, 298 insertions(+), 100 deletions(-)
diff --git a/cpp/full_code/NSVEparticles.cpp b/cpp/full_code/NSVEparticles.cpp
index ab05dc41..4d1ca6dc 100644
--- a/cpp/full_code/NSVEparticles.cpp
+++ b/cpp/full_code/NSVEparticles.cpp
@@ -92,7 +92,7 @@ int NSVEparticles<rnumber>::initialize(void)
this->particles_output_writer_mpi->template save<3>(
pset.getParticleState(),
rhs_data.data(),
- pset.getParticleIndex(),
+ pset.getParticleIndices(),
pset.getLocalNumberOfParticles(),
0);
this->particles_output_writer_mpi->close_file();
diff --git a/cpp/particles/abstract_particles_input.hpp b/cpp/particles/abstract_particles_input.hpp
index 2b581b61..20884283 100644
--- a/cpp/particles/abstract_particles_input.hpp
+++ b/cpp/particles/abstract_particles_input.hpp
@@ -39,6 +39,7 @@ public:
virtual std::unique_ptr<real_number[]> getMyParticles() = 0;
virtual std::unique_ptr<partsize_t[]> getMyParticlesIndexes() = 0;
+ virtual std::unique_ptr<partsize_t[]> getMyParticlesLabels() = 0;
virtual std::vector<std::unique_ptr<real_number[]>> getMyRhs() = 0;
virtual std::vector<hsize_t> getParticleFileLayout() = 0;
};
diff --git a/cpp/particles/abstract_particles_output.hpp b/cpp/particles/abstract_particles_output.hpp
index cadc65ce..af690dcb 100644
--- a/cpp/particles/abstract_particles_output.hpp
+++ b/cpp/particles/abstract_particles_output.hpp
@@ -97,7 +97,7 @@ protected:
}
public:
- abstract_particles_output(MPI_Comm in_mpi_com, const partsize_t inTotalNbParticles, const int in_nb_rhs) throw()
+ abstract_particles_output(MPI_Comm in_mpi_com, const partsize_t inTotalNbParticles, const int in_nb_rhs) noexcept(false)
: mpi_com(in_mpi_com), my_rank(-1), nb_processes(-1),
total_nb_particles(inTotalNbParticles), nb_rhs(in_nb_rhs),
buffer_particles_rhs_send(in_nb_rhs), size_buffers_send(0),
@@ -223,10 +223,14 @@ public:
buffer_indexes_send[idx_part].second = index_particles[idx_part];
}
- std::sort(&buffer_indexes_send[0], &buffer_indexes_send[nb_particles], [](const std::pair<partsize_t,partsize_t>& p1,
- const std::pair<partsize_t,partsize_t>& p2){
- return p1.second < p2.second;
- });
+ std::sort(
+ &buffer_indexes_send[0],
+ &buffer_indexes_send[nb_particles],
+ [](const std::pair<partsize_t,partsize_t>& p1,
+ const std::pair<partsize_t,partsize_t>& p2)
+ {
+ return p1.second < p2.second;
+ });
for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
const partsize_t src_idx = buffer_indexes_send[idx_part].first;
@@ -258,6 +262,8 @@ public:
// nb_particles_to_send is invalid after here
const int nb_to_receive = exchanger.getTotalToRecv();
+ //DEBUG_MSG("nb_to_receive = %d, particles_chunk_current_size = %d\n",
+ // nb_to_receive, particles_chunk_current_size);
assert(nb_to_receive == particles_chunk_current_size);
if(size_buffers_recv < nb_to_receive){
@@ -306,13 +312,22 @@ public:
}
{
+ // Local process has received particles in `buffer_particles_positions_recv`.
+ // It will now place this data (which is shuffled), in order, in the `buffer_particles_positions_send` array.
+ // However, to place the data it uses the particle indices as memory addresses.
+ // This is fine in general, but if we use particle indices as labels, and we sometimes delete particles, the particle indices are no longer valid memory addresses.
+ // So we need to fix this.
TIMEZONE("copy-local-order");
for(partsize_t idx_part = 0 ; idx_part < nb_to_receive ; ++idx_part){
const partsize_t src_idx = idx_part;
+ //DEBUG_MSG("idx_part = %d, buffer_indexes_recv[idx_part] = %d, particles_chunk_current_offset = %d\n",
+ // idx_part,
+ // buffer_indexes_recv[idx_part],
+ // particles_chunk_current_size);
const partsize_t dst_idx = buffer_indexes_recv[idx_part]-particles_chunk_current_offset;
+ //DEBUG_MSG("dst_idx is %d, particles_chunk_current_size is %d\n",
+ // dst_idx, particles_chunk_current_size);
assert(0 <= dst_idx);
- DEBUG_MSG("dst_idx is %d, particles_chunk_current_size is %d\n",
- dst_idx, particles_chunk_current_size);
assert(dst_idx < particles_chunk_current_size);
for(int idx_val = 0 ; idx_val < size_particle_positions ; ++idx_val){
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index c0a1128e..f7188e62 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -61,7 +61,8 @@ class abstract_particle_set
// extract particle information
virtual particle_rnumber* getParticleState() const = 0;
- virtual partsize_t* getParticleIndex() const = 0;
+ virtual partsize_t* getParticleIndices() const = 0;
+ virtual partsize_t* getParticleLabels() const = 0;
virtual int setParticleState(particle_rnumber *) = 0;
virtual int getParticleState(particle_rnumber *) = 0;
@@ -184,20 +185,17 @@ class abstract_particle_set
this->getParticleState(),
additional_data.data(),
int(additional_data.size()),
- this->getParticleIndex());
+ this->getParticleIndices());
return EXIT_SUCCESS;
}
- virtual int deleteParticles(
- const std::vector<partsize_t> &inIndicesToDelete) = 0;
-
- int writeParticleIndex(
+ int writeParticleLabels(
const std::string file_name,
const std::string species_name,
const std::string field_name,
const int iteration)
{
- TIMEZONE("abstract_particle_set::writeParticleIndex");
+ TIMEZONE("abstract_particle_set::writeParticleLabels");
particles_output_sampling_hdf5<partsize_t, particle_rnumber, partsize_t, 3> *particle_sample_writer = new particles_output_sampling_hdf5<partsize_t, double, partsize_t, 3>(
MPI_COMM_WORLD,
this->getTotalNumberOfParticles(),
@@ -211,15 +209,15 @@ class abstract_particle_set
// allocate temporary array
std::unique_ptr<partsize_t[]> pdata(new partsize_t[this->getLocalNumberOfParticles()]);
// clean up temporary array
- std::copy(this->getParticleIndex(),
- this->getParticleIndex() + this->getLocalNumberOfParticles(),
+ std::copy(this->getParticleLabels(),
+ this->getParticleLabels() + this->getLocalNumberOfParticles(),
pdata.get());
particle_sample_writer->template save_dataset<1>(
species_name,
field_name,
xx.get(),
&pdata,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
// deallocate temporary array
@@ -257,7 +255,7 @@ class abstract_particle_set
field_name,
xx.get(),
&pdata,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
// deallocate temporary array
@@ -327,7 +325,7 @@ class abstract_particle_set
field_name,
xx.get(),
&yy,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
break;
@@ -337,7 +335,7 @@ class abstract_particle_set
field_name,
xx.get(),
&yy,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
break;
@@ -347,7 +345,7 @@ class abstract_particle_set
field_name,
xx.get(),
&yy,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
break;
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 5d41c1c6..14ce618a 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -74,7 +74,8 @@ class particle_set: public abstract_particle_set
std::unique_ptr<partsize_t[]> offset_particles_for_partition;
std::unique_ptr<particle_rnumber[]> local_state;
- std::unique_ptr<partsize_t[]> local_index;
+ std::unique_ptr<partsize_t[]> local_indices;
+ std::unique_ptr<partsize_t[]> local_labels;
partsize_t local_number_of_particles;
partsize_t total_number_of_particles;
@@ -207,10 +208,20 @@ class particle_set: public abstract_particle_set
DEBUG_MSG("total_number_of_particles = %d\n",
this->total_number_of_particles);
- for (int i=0; i<this->particle_file_layout.size(); i++)
+ for (int i=0; i<int(this->particle_file_layout.size()); i++)
DEBUG_MSG("particle_file_layout[%d] = %d\n",
i,
this->particle_file_layout[i]);
+ if (this->total_number_of_particles < 100)
+ {
+ for (partsize_t ii=0; ii < this->local_number_of_particles; ii++)
+ {
+ DEBUG_MSG("local_counter = %d, index = %d, label = %d\n",
+ ii,
+ this->local_indices[ii],
+ this->local_labels[ii]);
+ }
+ }
DEBUG_MSG("# particle_set debug info ends\n");
DEBUG_MSG("##########################################\n");
@@ -245,9 +256,14 @@ class particle_set: public abstract_particle_set
return EXIT_SUCCESS;
}
- partsize_t* getParticleIndex() const
+ partsize_t* getParticleIndices() const
{
- return this->local_index.get();
+ return this->local_indices.get();
+ }
+
+ partsize_t* getParticleLabels() const
+ {
+ return this->local_labels.get();
}
partsize_t getLocalNumberOfParticles() const
@@ -367,7 +383,7 @@ class particle_set: public abstract_particle_set
&this->local_state,
additional_data.data(),
int(additional_data.size()),
- &this->local_index);
+ &this->local_indices);
return EXIT_SUCCESS;
}
@@ -376,10 +392,12 @@ class particle_set: public abstract_particle_set
TIMEZONE("particle_set::init");
this->local_state = particles_input.getMyParticles();
- this->local_index = particles_input.getMyParticlesIndexes();
+ this->local_indices = particles_input.getMyParticlesIndexes();
+ this->local_labels = particles_input.getMyParticlesLabels();
this->local_number_of_particles = particles_input.getLocalNbParticles();
this->total_number_of_particles = particles_input.getTotalNbParticles();
+
particles_utils::partition_extra_z<partsize_t, state_size>(
&this->local_state[0],
this->local_number_of_particles,
@@ -392,7 +410,8 @@ class particle_set: public abstract_particle_set
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]);
+ std::swap(this->local_indices[idx1], this->local_indices[idx2]);
+ std::swap(this->local_labels[idx1], this->local_labels[idx2]);
});
int file_layout_result = this->setParticleFileLayout(particles_input.getParticleFileLayout());
@@ -409,7 +428,8 @@ class particle_set: public abstract_particle_set
TIMEZONE("particle_set::init_as_subset_of");
assert(indices_to_copy.size() > 0);
particle_rnumber *tmp_local_state = new particle_rnumber[state_size * src.getLocalNbParticles()];
- partsize_t *tmp_local_index = new partsize_t[src.getLocalNbParticles()];
+ partsize_t *tmp_local_indices = new partsize_t[src.getLocalNbParticles()];
+ partsize_t *tmp_local_labels = new partsize_t[src.getLocalNbParticles()];
this->local_number_of_particles = 0;
this->total_number_of_particles = indices_to_copy.size();
@@ -421,13 +441,13 @@ class particle_set: public abstract_particle_set
// dumb selection of interesting particles
for (partsize_t ii = 0; ii < partsize_t(src.getLocalNbParticles()); ii++)
{
- partsize_t src_index = src.getParticlesIndexes()[ii];
+ partsize_t src_label = src.getParticlesIndexes()[ii];
for (partsize_t iii=0; iii < partsize_t(indices_to_copy.size()); iii++)
{
- if (src_index == indices_to_copy[iii])
+ if (src_label == indices_to_copy[iii])
{
- tmp_local_index[this->local_number_of_particles] = src_index;
- //tmp_local_index[this->local_number_of_particles] = iii;
+ tmp_local_indices[this->local_number_of_particles] = iii;
+ tmp_local_labels[this->local_number_of_particles] = src_label;
std::copy(src.getParticlesState()+state_size*ii,
src.getParticlesState()+state_size*(ii+1),
tmp_local_state + state_size*this->local_number_of_particles);
@@ -441,13 +461,207 @@ class particle_set: public abstract_particle_set
if (this->local_number_of_particles > 0)
{
this->local_state.reset(new particle_rnumber[state_size*this->local_number_of_particles]);
- this->local_index.reset(new partsize_t[this->local_number_of_particles]);
+ this->local_indices.reset(new partsize_t[this->local_number_of_particles]);
+ this->local_labels.reset(new partsize_t[this->local_number_of_particles]);
+ std::copy(tmp_local_state,
+ tmp_local_state + state_size*this->local_number_of_particles,
+ this->local_state.get());
+ std::copy(tmp_local_indices,
+ tmp_local_indices + this->local_number_of_particles,
+ this->local_indices.get());
+ std::copy(tmp_local_labels,
+ tmp_local_labels + this->local_number_of_particles,
+ this->local_labels.get());
+ }
+
+ particles_utils::partition_extra_z<partsize_t, state_size>(
+ &this->local_state[0],
+ this->local_number_of_particles,
+ partition_interval_size,
+ this->number_particles_per_partition.get(),
+ this->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_indices[idx1], this->local_indices[idx2]);
+ std::swap(this->local_labels[idx1], this->local_labels[idx2]);
+ });
+
+ delete[] tmp_local_state;
+ delete[] tmp_local_indices;
+ delete[] tmp_local_labels;
+ return EXIT_SUCCESS;
+ }
+
+ int operator=(
+ abstract_particle_set *src)
+ {
+ TIMEZONE("particle_set::operator=");
+ assert(src->getStateSize() == state_size);
+ particle_rnumber *tmp_local_state = new particle_rnumber[state_size * src->getLocalNumberOfParticles()];
+ partsize_t *tmp_local_indices = new partsize_t[src->getLocalNumberOfParticles()];
+ partsize_t *tmp_local_labels = new partsize_t[src->getLocalNumberOfParticles()];
+ this->local_number_of_particles = 0;
+
+ this->total_number_of_particles = src->getTotalNumberOfParticles();
+
+ for (partsize_t ii = 0; ii < partsize_t(src->getLocalNumberOfParticles()); ii++)
+ {
+ partsize_t src_label = src->getParticleLabels()[ii];
+ partsize_t src_index = src->getParticleIndices()[ii];
+
+ {
+ // set new index
+ tmp_local_indices[this->local_number_of_particles] = src_index;
+ tmp_local_labels[this->local_number_of_particles] = src_label;
+ std::copy(src->getParticleState()+state_size*ii,
+ src->getParticleState()+state_size*(ii+1),
+ tmp_local_state + state_size*this->local_number_of_particles);
+ this->local_number_of_particles++;
+ }
+ }
+
+ // now we actually put the data "here"
+ if (this->local_number_of_particles > 0)
+ {
+ this->local_state.reset(new particle_rnumber[state_size*this->local_number_of_particles]);
+ this->local_indices.reset(new partsize_t[this->local_number_of_particles]);
+ this->local_labels.reset(new partsize_t[this->local_number_of_particles]);
+ std::copy(tmp_local_state,
+ tmp_local_state + state_size*this->local_number_of_particles,
+ this->local_state.get());
+ std::copy(tmp_local_indices,
+ tmp_local_indices + this->local_number_of_particles,
+ this->local_indices.get());
+ std::copy(tmp_local_labels,
+ tmp_local_labels + this->local_number_of_particles,
+ this->local_labels.get());
+ }
+
+ particles_utils::partition_extra_z<partsize_t, state_size>(
+ &this->local_state[0],
+ this->local_number_of_particles,
+ partition_interval_size,
+ this->number_particles_per_partition.get(),
+ this->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_indices[idx1], this->local_indices[idx2]);
+ std::swap(this->local_labels[idx1], this->local_labels[idx2]);
+ });
+
+ delete[] tmp_local_state;
+ delete[] tmp_local_indices;
+ delete[] tmp_local_labels;
+ return EXIT_SUCCESS;
+ }
+
+ int reset_labels()
+ {
+ std::copy(
+ this->local_indices.get(),
+ this->local_indices.get() + this->getLocalNumberOfParticles(),
+ this->local_labels.get());
+ return EXIT_SUCCESS;
+ }
+
+ /** \brief Filter particle set versus small list of particle indices.
+ *
+ * Using this method, we can select particles from a `src` particle set into the current object.
+ * The code will work reasonably as long as a small list of "indices of interest" is used.
+ *
+ * \warning The indices of interest must be existing indices, otherwise the object will be broken.
+ * In principle we could have a general algorithm that could handle anything, but it would be a messy code,
+ * which would need multiple passes of the data to figure out the new indices.
+ * It's easier to just clean up the data before feeding it to this method.
+ *
+ */
+ int init_as_subset_of(
+ abstract_particle_set &src,
+ const std::vector<partsize_t> indices_of_interest,
+ bool use_set_complement)
+ {
+ TIMEZONE("particle_set::init_as_subset_of version 2");
+ assert(indices_of_interest.size() > 0);
+ assert(partsize_t(indices_of_interest.size()) < src.getTotalNumberOfParticles());
+ particle_rnumber *tmp_local_state = new particle_rnumber[state_size * src.getLocalNumberOfParticles()];
+ partsize_t *tmp_local_indices = new partsize_t[src.getLocalNumberOfParticles()];
+ partsize_t *tmp_local_labels = new partsize_t[src.getLocalNumberOfParticles()];
+ this->local_number_of_particles = 0;
+
+ // set new number of particles
+ if (use_set_complement)
+ this->total_number_of_particles = src.getTotalNumberOfParticles() - indices_of_interest.size();
+ else
+ this->total_number_of_particles = indices_of_interest.size();
+
+ // select interesting particles
+ for (partsize_t ii = 0; ii < partsize_t(src.getLocalNumberOfParticles()); ii++)
+ {
+ partsize_t src_label = src.getParticleLabels()[ii];
+ partsize_t src_index = src.getParticleIndices()[ii];
+
+ bool index_is_interesting = false;
+ // find out if index is interesting
+ partsize_t iii = 0;
+ for (; iii < partsize_t(indices_of_interest.size()); iii++)
+ {
+ index_is_interesting = (indices_of_interest[iii] == src_index);
+ if (indices_of_interest[iii] >= src_index)
+ break;
+ }
+ bool copy_data = false;
+
+ if (use_set_complement && (!index_is_interesting)) // we are effectively deleting particles present in "indices_of_interest"
+ copy_data = true;
+ if ((!use_set_complement) && index_is_interesting) // we are only keeping particles present in "indices_of_interest"
+ copy_data = true;
+ if (copy_data)
+ {
+ // set new index
+ tmp_local_indices[this->local_number_of_particles] = src_index;
+ if (use_set_complement)
+ tmp_local_indices[this->local_number_of_particles] -= iii;
+ else
+ tmp_local_indices[this->local_number_of_particles] -= src.getTotalNumberOfParticles() - iii;
+ // copy particle with label src_label
+ tmp_local_labels[this->local_number_of_particles] = src_label;
+ std::copy(src.getParticleState()+state_size*ii,
+ src.getParticleState()+state_size*(ii+1),
+ tmp_local_state + state_size*this->local_number_of_particles);
+ this->local_number_of_particles++;
+ }
+ }
+ DEBUG_MSG("particle_set::init_as_subset_of --- after particle selection local number of particles is %d\n",
+ this->local_number_of_particles);
+ for (partsize_t ii = 0; ii < this->local_number_of_particles; ii++)
+ DEBUG_MSG(" --- ii = %d, tmp_local_labels[ii] = %d, tmp_local_indices[ii] = %d\n",
+ ii,
+ tmp_local_labels[ii],
+ tmp_local_indices[ii]);
+
+ // now we actually put the data "here"
+ if (this->local_number_of_particles > 0)
+ {
+ this->local_state.reset(new particle_rnumber[state_size*this->local_number_of_particles]);
+ this->local_indices.reset(new partsize_t[this->local_number_of_particles]);
+ this->local_labels.reset(new partsize_t[this->local_number_of_particles]);
std::copy(tmp_local_state,
tmp_local_state + state_size*this->local_number_of_particles,
this->local_state.get());
- std::copy(tmp_local_index,
- tmp_local_index + this->local_number_of_particles,
- this->local_index.get());
+ std::copy(tmp_local_indices,
+ tmp_local_indices + this->local_number_of_particles,
+ this->local_indices.get());
+ std::copy(tmp_local_labels,
+ tmp_local_labels + this->local_number_of_particles,
+ this->local_labels.get());
}
particles_utils::partition_extra_z<partsize_t, state_size>(
@@ -462,11 +676,13 @@ class particle_set: public abstract_particle_set
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]);
+ std::swap(this->local_indices[idx1], this->local_indices[idx2]);
+ std::swap(this->local_labels[idx1], this->local_labels[idx2]);
});
delete[] tmp_local_state;
- delete[] tmp_local_index;
+ delete[] tmp_local_indices;
+ delete[] tmp_local_labels;
return EXIT_SUCCESS;
}
@@ -490,70 +706,13 @@ class particle_set: public abstract_particle_set
{
std::vector<partsize_t> local_indices;
for (partsize_t i=0; i<this->getLocalNumberOfParticles(); i++)
- for (partsize_t ii=0; ii<inGlobalIndicesToDelete.size(); ii++)
+ for (partsize_t ii=0; ii<partsize_t(inGlobalIndicesToDelete.size()); ii++)
{
- if (this->getParticleIndex()[i] == inGlobalIndicesToDelete[ii])
+ if (this->getParticleIndices()[i] == inGlobalIndicesToDelete[ii])
local_indices.push_back(inGlobalIndicesToDelete[ii]);
}
return local_indices;
}
-
- /** \brief Delete particles from arbitrary list
- *
- * \note List of particles given to this method may contain indices to
- * particles at arbitrary locations.
- */
- int deleteParticles(
- const std::vector<partsize_t> &inGlobalIndicesToDelete)
- {
- std::vector<partsize_t> inLocalIndicesToDelete =
- this->selectLocalParticleFromFewIndices(inGlobalIndicesToDelete);
- std::set<partsize_t> setIndicesToDelete(inLocalIndicesToDelete.begin(), inLocalIndicesToDelete.end());
- if (setIndicesToDelete.size() > 0)
- {
- partsize_t number_of_deleted_particles = 0;
-
- for(int idxPartition = 0 ; idxPartition < partition_interval_size ; ++idxPartition){
- const partsize_t nbDeletedParticlesBefore = number_of_deleted_particles;
- for (partsize_t idx = this->offset_particles_for_partition[idxPartition];
- idx < this->offset_particles_for_partition[idxPartition+1];
- ++idx)
- {
- if(setIndicesToDelete.find(this->local_index[idx]) != setIndicesToDelete.end())
- {
- number_of_deleted_particles += 1;
- }
- else if(number_of_deleted_particles)
- {
- this->local_index[idx-number_of_deleted_particles] = this->local_index[idx];
- for(int idx_pos = 0 ; idx_pos < state_size; ++idx_pos)
- {
- this->local_state[(idx-number_of_deleted_particles)*state_size+idx_pos] =
- this->local_state[idx*state_size+idx_pos];
- }
- }
- }
-
- this->offset_particles_for_partition[idxPartition] -= nbDeletedParticlesBefore;
- }
-
- this->offset_particles_for_partition[partition_interval_size] -= number_of_deleted_particles;
-
- this->local_number_of_particles -= number_of_deleted_particles;
- assert(this->local_number_of_particles == this->offset_particles_for_partition[partition_interval_size]);
- }
-
- AssertMpi(MPI_Allreduce(
- const_cast<partsize_t*>(&(this->local_number_of_particles)),
- &(this->total_number_of_particles),
- 1,
- particles_utils::GetMpiType(partsize_t()),
- MPI_SUM,
- this->mpi_comm));
- this->particle_file_layout.resize(1);
- this->particle_file_layout[0] = this->total_number_of_particles;
- return EXIT_SUCCESS;
- }
};
#endif//PARTICLE_SET_HPP
diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp
index 45390268..7534562c 100644
--- a/cpp/particles/particle_solver.cpp
+++ b/cpp/particles/particle_solver.cpp
@@ -174,7 +174,7 @@ int particle_solver::writeCheckpoint(
particles_output_writer->template save<state_size>(
this->pset->getParticleState(),
this->additional_states.data(),
- this->pset->getParticleIndex(),
+ this->pset->getParticleIndices(),
this->pset->getLocalNumberOfParticles(),
this->iteration);
return EXIT_SUCCESS;
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 3cb04ce3..97fad59d 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -54,6 +54,7 @@ class particles_input_hdf5 : public abstract_particles_input<partsize_t, real_nu
std::unique_ptr<real_number[]> my_particles_positions;
std::unique_ptr<partsize_t[]> my_particles_indexes;
+ std::unique_ptr<partsize_t[]> my_particles_labels;
std::vector<std::unique_ptr<real_number[]>> my_particles_rhs;
public:
@@ -318,6 +319,12 @@ public:
hdfret = H5Pclose(plist_id_par);
assert(hdfret >= 0);
}
+
+ // clone indices in labels:
+ my_particles_labels.reset(new partsize_t[this->getLocalNbParticles()]);
+ std::copy(my_particles_indexes.get(),
+ my_particles_indexes.get()+this->getLocalNbParticles(),
+ my_particles_labels.get());
}
~particles_input_hdf5() noexcept(false){
@@ -350,6 +357,11 @@ public:
return std::move(my_particles_indexes);
}
+ std::unique_ptr<partsize_t[]> getMyParticlesLabels() final {
+ assert(my_particles_labels != nullptr || nb_particles_for_me == 0);
+ return std::move(my_particles_labels);
+ }
+
std::vector<hsize_t> getParticleFileLayout(){
return std::vector<hsize_t>(this->particle_file_layout);
}
diff --git a/cpp/particles/particles_input_random.hpp b/cpp/particles/particles_input_random.hpp
index 1179b080..0b40ae57 100644
--- a/cpp/particles/particles_input_random.hpp
+++ b/cpp/particles/particles_input_random.hpp
@@ -44,6 +44,7 @@ class particles_input_random: public abstract_particles_input<partsize_t, partic
std::unique_ptr<particle_rnumber[]> local_particle_state;
std::unique_ptr<partsize_t[]> local_particle_index;
+ std::unique_ptr<partsize_t[]> local_particle_label;
std::vector<std::unique_ptr<particle_rnumber[]>> local_particle_rhs;
public:
particles_input_random(
@@ -219,6 +220,12 @@ class particles_input_random: public abstract_particles_input<partsize_t, partic
}
/*************************************************/
+
+ // clone indices in labels:
+ local_particle_label.reset(new partsize_t[this->getLocalNbParticles()]);
+ std::copy(local_particle_index.get(),
+ local_particle_index.get()+this->getLocalNbParticles(),
+ local_particle_label.get());
}
partsize_t getTotalNbParticles()
@@ -246,6 +253,12 @@ class particles_input_random: public abstract_particles_input<partsize_t, partic
return std::move(this->local_particle_index);
}
+ std::unique_ptr<partsize_t[]> getMyParticlesLabels()
+ {
+ assert(this->local_particle_label != nullptr || this->local_number_of_particles == 0);
+ return std::move(this->local_particle_label);
+ }
+
std::vector<std::unique_ptr<particle_rnumber[]>> getMyRhs()
{
return std::move(std::vector<std::unique_ptr<particle_rnumber[]>>());
--
GitLab
From 739f20bac871649dc7bc497b1318f48586cd9942 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 21:29:35 +0100
Subject: [PATCH 15/37] [wip] debugging particle delete functionality
---
TurTLE/test/particle_set/particle_deleter.cpp | 113 +++++++++++-------
TurTLE/test/test_particle_deleter.py | 12 +-
2 files changed, 76 insertions(+), 49 deletions(-)
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
index af7d93ca..fb20ca26 100644
--- a/TurTLE/test/particle_set/particle_deleter.cpp
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -159,12 +159,12 @@ int particle_deleter<rnumber>::do_stats()
return EXIT_SUCCESS;
// sample position
- this->pset->writeStateTriplet(
- 0,
- this->particles_sample_writer_mpi,
- "tracers0",
- "position",
- psolver->getIteration());
+ //this->pset->writeStateTriplet(
+ // 0,
+ // this->particles_sample_writer_mpi,
+ // "tracers0",
+ // "position",
+ // psolver->getIteration());
// sample velocity
// first ensure velocity field is computed, and is in real space
@@ -175,44 +175,68 @@ int particle_deleter<rnumber>::do_stats()
*this->tmp_vec_field = this->fs->cvelocity->get_cdata();
this->tmp_vec_field->ift();
}
- this->pset->writeSample(
- this->tmp_vec_field,
- this->particles_sample_writer_mpi,
- "tracers0",
- "velocity",
- psolver->getIteration());
-
- this->pset->writeParticleIndex(
- "particle_index.h5",
- "tracers0",
- "index",
- psolver->getIteration());
-
- std::vector<long long int> indices_to_delete(10);
- indices_to_delete[0] = 1;
- indices_to_delete[1] = 3;
- indices_to_delete[2] = 19;
- indices_to_delete[3] = 13;
- indices_to_delete[4] = 101;
- indices_to_delete[5] = 91;
- indices_to_delete[6] = 141;
- indices_to_delete[7] = 143;
- indices_to_delete[8] = 145;
- indices_to_delete[9] = 149;
-
- DEBUG_MSG("before delete particles");
- this->pset->_print_debug_info();
- this->pset->deleteParticles(indices_to_delete);
- DEBUG_MSG("after delete particles");
- this->pset->_print_debug_info();
-
- DEBUG_MSG("before writeParticleIndex");
- this->pset->writeParticleIndex(
- "particle_index.h5",
- "tracers0",
- "index",
- psolver->getIteration()+1);
- DEBUG_MSG("after writeParticleIndex");
+ //this->pset->writeSample(
+ // this->tmp_vec_field,
+ // this->particles_sample_writer_mpi,
+ // "tracers0",
+ // "velocity",
+ // psolver->getIteration());
+
+ if (this->pset->getTotalNumberOfParticles() > 3)
+ {
+ //this->pset->writeParticleLabels(
+ // "particle_index_before.h5",
+ // "tracers0",
+ // "index",
+ // psolver->getIteration());
+ std::vector<long long int> indices_to_delete(2);
+ indices_to_delete[0] = 1;
+ indices_to_delete[1] = (this->iteration*3)%this->pset->getTotalNumberOfParticles();
+ if (indices_to_delete[1] == indices_to_delete[0])
+ indices_to_delete[1] = 3;
+
+ DEBUG_MSG("before delete particles\n");
+ this->pset->_print_debug_info();
+ particle_set<3, 3, 2> *tmp_pset = new particle_set<3, 3, 2>(
+ this->fs->cvelocity->rlayout,
+ this->dkx,
+ this->dky,
+ this->dkz,
+ 0.5);
+ tmp_pset->init_as_subset_of(
+ *(this->pset),
+ indices_to_delete,
+ true);
+ *(this->pset) = tmp_pset;
+ delete tmp_pset;
+
+ DEBUG_MSG("after delete particles\n");
+ this->pset->_print_debug_info();
+
+ delete this->particles_output_writer_mpi;
+ delete this->particles_sample_writer_mpi;
+ this->particles_output_writer_mpi = new particles_output_hdf5<
+ long long int, double, 3>(
+ MPI_COMM_WORLD,
+ "tracers0",
+ this->pset->getTotalNumberOfParticles(),
+ tracers0_integration_steps);
+ this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+ long long int, double, double, 3>(
+ MPI_COMM_WORLD,
+ this->pset->getTotalNumberOfParticles(),
+ (this->simname + "_particles.h5"),
+ "tracers0",
+ "position/0");
+
+ DEBUG_MSG("before writeParticleIndex");
+ //this->pset->writeParticleLabels(
+ // "particle_index_after.h5",
+ // "tracers0",
+ // "index",
+ // psolver->getIteration());
+ DEBUG_MSG("after writeParticleIndex");
+ }
return EXIT_SUCCESS;
}
@@ -237,3 +261,4 @@ int particle_deleter<rnumber>::read_parameters(void)
template class particle_deleter<float>;
template class particle_deleter<double>;
+
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index 9e98271e..ce2b0ab7 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -77,9 +77,9 @@ class ADNS(TurTLE.DNS):
self.parameters['dkx'] = float(1.0)
self.parameters['dky'] = float(1.0)
self.parameters['dkz'] = float(1.0)
- self.parameters['niter_todo'] = int(200)
- self.parameters['niter_stat'] = int(20)
- self.parameters['niter_out'] = int(200)
+ self.parameters['niter_todo'] = int(20)
+ self.parameters['niter_stat'] = int(2)
+ self.parameters['niter_out'] = int(20)
self.parameters['checkpoints_per_file'] = int(1)
self.parameters['dt'] = float(0.0001)
self.parameters['histogram_bins'] = int(64)
@@ -88,7 +88,7 @@ class ADNS(TurTLE.DNS):
self.parameters['niter_part'] = int(1)
self.parameters['niter_part_fine_period'] = int(2)
self.parameters['niter_part_fine_duration'] = int(0)
- self.parameters['nparticles'] = int(1000)
+ self.parameters['nparticles'] = int(10)
self.parameters['tracers0_integration_steps'] = int(0)
self.parameters['tracers0_neighbours'] = int(3)
self.parameters['tracers0_smoothness'] = int(2)
@@ -308,7 +308,9 @@ class ADNS(TurTLE.DNS):
particle_file.create_group('tracers0/position')
particle_file.create_group('tracers0/velocity')
particle_file.create_group('collisions0')
- with h5py.File('particle_index.h5', 'w') as particle_file:
+ with h5py.File('particle_index_before.h5', 'w') as particle_file:
+ particle_file.create_group('tracers0/index')
+ with h5py.File('particle_index_after.h5', 'w') as particle_file:
particle_file.create_group('tracers0/index')
self.generate_tracer_state(integration_steps = self.parameters['tracers0_integration_steps'],
rseed = opt.particle_rand_seed,
--
GitLab
From a955fd34879be4da7a2d23870b45bb28c96aa4ba Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 1 Dec 2021 23:03:23 +0100
Subject: [PATCH 16/37] fixes memory leak
debug message hardcoded for more MPI processes than used during test.
for some reason, this didn't show up until now in the unit tests (but it
should have).
luckily, fsanitize=address catches it.
---
cpp/particles/p2p/p2p_distr_mpi.hpp | 15 ++++++++++-----
1 file changed, 10 insertions(+), 5 deletions(-)
diff --git a/cpp/particles/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index ecf6aa4d..06f93fbf 100644
--- a/cpp/particles/p2p/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p/p2p_distr_mpi.hpp
@@ -434,13 +434,18 @@ public:
DEBUG_MSG("my_down_z_cell_level = %d\n", int(my_down_z_cell_level));
DEBUG_MSG("first_cell_level_proc(0) = %d\n", int(first_cell_level_proc(0)));
DEBUG_MSG("first_cell_level_proc(1) = %d\n", int(first_cell_level_proc(1)));
- DEBUG_MSG("first_cell_level_proc(2) = %d\n", int(first_cell_level_proc(2)));
- DEBUG_MSG("first_cell_level_proc(3) = %d\n", int(first_cell_level_proc(3)));
+ if (nb_processes > 1)
+ DEBUG_MSG("first_cell_level_proc(2) = %d\n", int(first_cell_level_proc(2)));
+ if (nb_processes > 2)
+ DEBUG_MSG("first_cell_level_proc(3) = %d\n", int(first_cell_level_proc(3)));
DEBUG_MSG("nb_cell_levels[IDXC_Z] = %d\n", int(nb_cell_levels[IDXC_Z]));
DEBUG_MSG("last_cell_level_proc(0) = %d\n", int(last_cell_level_proc(0)));
- DEBUG_MSG("last_cell_level_proc(1) = %d\n", int(last_cell_level_proc(1)));
- DEBUG_MSG("last_cell_level_proc(2) = %d\n", int(last_cell_level_proc(2)));
- DEBUG_MSG("last_cell_level_proc(3) = %d\n", int(last_cell_level_proc(3)));
+ if (nb_processes > 1)
+ DEBUG_MSG("last_cell_level_proc(1) = %d\n", int(last_cell_level_proc(1)));
+ if (nb_processes > 2)
+ DEBUG_MSG("last_cell_level_proc(2) = %d\n", int(last_cell_level_proc(2)));
+ if (nb_processes > 3)
+ DEBUG_MSG("last_cell_level_proc(3) = %d\n", int(last_cell_level_proc(3)));
// Find process with at least one neighbor
{
--
GitLab
From 62fbe4a40ceede6fc279cb68e3290ab15bc5d68e Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Thu, 2 Dec 2021 09:41:29 +0100
Subject: [PATCH 17/37] [wip] takes out "fsanitize=address"
HDF5 still triggers an error, so I need this flag out to check the
pipeline.
---
CMakeLists.txt | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 3dc4417c..27308ff9 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -230,7 +230,7 @@ endif()
#####################################################################################
## Extra flags
-set(CMAKE_CXX_COMPILE_FLAGS "${CMAKE_CXX_COMPILE_FLAGS} $ENV{TURTLE_COMPILATION_FLAGS} -Wall -g -fsanitize=address")
+set(CMAKE_CXX_COMPILE_FLAGS "${CMAKE_CXX_COMPILE_FLAGS} $ENV{TURTLE_COMPILATION_FLAGS} -Wall -g")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CMAKE_CXX_COMPILE_FLAGS}")
#####################################################################################
--
GitLab
From 0b8614a7c0e61cf713e2525c6dab88107f84d87b Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 4 Dec 2021 18:17:48 +0100
Subject: [PATCH 18/37] gets rid of warning from gcc
apparently the const for an event pair thing should be in a different
place in order to prevent copying.
---
cpp/scope_timer.hpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/cpp/scope_timer.hpp b/cpp/scope_timer.hpp
index 51fc2d38..7d13cce3 100644
--- a/cpp/scope_timer.hpp
+++ b/cpp/scope_timer.hpp
@@ -415,7 +415,7 @@ public:
std::vector<SerializedEvent> myEvents;
myEvents.reserve(m_records.size());
- for(const std::pair<std::string, const CoreEvent*>& event : m_records){
+ for(const std::pair<const std::string, CoreEvent*>& event : m_records){
myEvents.emplace_back();
SerializedEvent& current_event = myEvents.back();
--
GitLab
From 9c110ac85303302099fe1f8b350d2085e97add3e Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 4 Dec 2021 19:16:07 +0100
Subject: [PATCH 19/37] fixes memory leak for NSVEparticle_set.cpp (test)
---
TurTLE/test/particle_set/NSVEparticle_set.cpp | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/TurTLE/test/particle_set/NSVEparticle_set.cpp b/TurTLE/test/particle_set/NSVEparticle_set.cpp
index e2196a85..c6e8ba81 100644
--- a/TurTLE/test/particle_set/NSVEparticle_set.cpp
+++ b/TurTLE/test/particle_set/NSVEparticle_set.cpp
@@ -129,12 +129,13 @@ template <typename rnumber>
int NSVEparticle_set<rnumber>::finalize(void)
{
TIMEZONE("NSVEparticle_set::finalize");
- delete this->particles_output_writer_mpi;
delete this->particles_sample_writer_mpi;
+ delete this->particles_output_writer_mpi;
delete this->psolver;
delete this->pset;
delete this->trhs;
delete this->fti;
+ delete this->previous_velocity;
this->NSVE<rnumber>::finalize();
return EXIT_SUCCESS;
}
--
GitLab
From 68a6841ee458747ef3b148fc23c12e35f6937dc0 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 4 Dec 2021 19:17:57 +0100
Subject: [PATCH 20/37] fixes destructor (mem leak copied from
NSVEparticle_set)
---
TurTLE/test/particle_set/particle_deleter.cpp | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
index fb20ca26..07d073f1 100644
--- a/TurTLE/test/particle_set/particle_deleter.cpp
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -129,12 +129,13 @@ template <typename rnumber>
int particle_deleter<rnumber>::finalize(void)
{
TIMEZONE("particle_deleter::finalize");
- delete this->particles_output_writer_mpi;
delete this->particles_sample_writer_mpi;
+ delete this->particles_output_writer_mpi;
delete this->psolver;
delete this->pset;
delete this->trhs;
delete this->fti;
+ delete this->previous_velocity;
this->NSVE<rnumber>::finalize();
return EXIT_SUCCESS;
}
--
GitLab
From 17f51dbb019a23403b7bacb39118753203562adb Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sun, 5 Dec 2021 11:17:24 +0100
Subject: [PATCH 21/37] [wip] fixes delete particle functionality
still a bunch of details to implement for the test.
---
TurTLE/test/particle_set/particle_deleter.cpp | 53 +++++++++++--------
TurTLE/test/test_particle_deleter.py | 33 +-----------
.../interpolation/abstract_particle_set.hpp | 2 +
cpp/particles/interpolation/particle_set.hpp | 50 ++++++++++++-----
4 files changed, 69 insertions(+), 69 deletions(-)
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
index 07d073f1..165698a9 100644
--- a/TurTLE/test/particle_set/particle_deleter.cpp
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -160,12 +160,12 @@ int particle_deleter<rnumber>::do_stats()
return EXIT_SUCCESS;
// sample position
- //this->pset->writeStateTriplet(
- // 0,
- // this->particles_sample_writer_mpi,
- // "tracers0",
- // "position",
- // psolver->getIteration());
+ this->pset->writeStateTriplet(
+ 0,
+ this->particles_sample_writer_mpi,
+ "tracers0",
+ "position",
+ psolver->getIteration());
// sample velocity
// first ensure velocity field is computed, and is in real space
@@ -176,27 +176,34 @@ int particle_deleter<rnumber>::do_stats()
*this->tmp_vec_field = this->fs->cvelocity->get_cdata();
this->tmp_vec_field->ift();
}
- //this->pset->writeSample(
- // this->tmp_vec_field,
- // this->particles_sample_writer_mpi,
- // "tracers0",
- // "velocity",
- // psolver->getIteration());
+ this->pset->writeSample(
+ this->tmp_vec_field,
+ this->particles_sample_writer_mpi,
+ "tracers0",
+ "velocity",
+ psolver->getIteration());
if (this->pset->getTotalNumberOfParticles() > 3)
{
- //this->pset->writeParticleLabels(
- // "particle_index_before.h5",
- // "tracers0",
- // "index",
- // psolver->getIteration());
+ this->pset->writeParticleLabels(
+ "particle_index_before.h5",
+ "tracers0",
+ "index",
+ psolver->getIteration());
std::vector<long long int> indices_to_delete(2);
indices_to_delete[0] = 1;
indices_to_delete[1] = (this->iteration*3)%this->pset->getTotalNumberOfParticles();
if (indices_to_delete[1] == indices_to_delete[0])
indices_to_delete[1] = 3;
+ std::sort(indices_to_delete.begin(), indices_to_delete.end());
- DEBUG_MSG("before delete particles\n");
+ DEBUG_MSG("computed indices_to_delete:\n");
+ DEBUG_MSG("indices_to_delete[0] = %d\n",
+ indices_to_delete[0]);
+ DEBUG_MSG("indices_to_delete[1] = %d\n",
+ indices_to_delete[1]);
+ DEBUG_MSG("before delete particles, iteration %d\n",
+ this->iteration);
this->pset->_print_debug_info();
particle_set<3, 3, 2> *tmp_pset = new particle_set<3, 3, 2>(
this->fs->cvelocity->rlayout,
@@ -231,11 +238,11 @@ int particle_deleter<rnumber>::do_stats()
"position/0");
DEBUG_MSG("before writeParticleIndex");
- //this->pset->writeParticleLabels(
- // "particle_index_after.h5",
- // "tracers0",
- // "index",
- // psolver->getIteration());
+ this->pset->writeParticleLabels(
+ "particle_index_after.h5",
+ "tracers0",
+ "index",
+ psolver->getIteration());
DEBUG_MSG("after writeParticleIndex");
}
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index ce2b0ab7..c9a0e483 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -328,39 +328,8 @@ class ADNS(TurTLE.DNS):
def main():
bla = ADNS()
bla.launch(sys.argv[1:])
- assert(bla.parameters['niter_part'] == 1)
- df = bla.get_particle_file()
- x = []
- v = []
- for ii in range(bla.parameters['niter_out']):
- x.append(df['tracers0/position/{0}'.format(ii)][()])
- v.append(df['tracers0/velocity/{0}'.format(ii)][()])
- df.close()
- x = np.array(x)
- v = np.array(v)
- d1x = (x[2:] - x[:-2]) / (2*bla.parameters['dt'])
- vv = v[1:-1]
- diff = vv - d1x
- relative_error = (np.sum(diff**2, axis = 2) / np.sum(vv**2, axis = 2))**0.5
- tindex, trajindex = np.unravel_index(
- np.argmax(relative_error),
- relative_error.shape)
- try:
- import matplotlib.pyplot as plt
- f = plt.figure()
- a = f.add_subplot(211)
- a.plot(vv[:, trajindex])
- a.plot(d1x[:, trajindex], dashes = (2, 2))
- a = f.add_subplot(212)
- a.plot(x[:, trajindex])
- a.set_title('maximum error for (tindex, trajindex) = ({0}, {1})'.format(tindex, trajindex))
- f.tight_layout()
- f.savefig('traj.pdf')
- except ImportError:
- print('couldn\'t import matplotlib, no figure')
- print(np.max(relative_error))
- assert(np.max(relative_error) < 1e-4)
return None
+
if __name__ == '__main__':
main()
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index f7188e62..f5481cd1 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -220,6 +220,8 @@ class abstract_particle_set
this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
+ // deallocate sample writer
+ delete particle_sample_writer;
// deallocate temporary array
delete[] pdata.release();
// deallocate position array
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 14ce618a..dff01348 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -161,9 +161,9 @@ class particle_set: public abstract_particle_set
assert(IDXC_Y == 1);
assert(IDXC_Z == 2);
- DEBUG_MSG("current partition interval %d %d\n",
- current_partition_interval.first,
- current_partition_interval.second);
+ //DEBUG_MSG("current partition interval %d %d\n",
+ // current_partition_interval.first,
+ // current_partition_interval.second);
this->number_particles_per_partition.reset(new partsize_t[partition_interval_size]);
this->offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
@@ -433,10 +433,10 @@ class particle_set: public abstract_particle_set
this->local_number_of_particles = 0;
this->total_number_of_particles = indices_to_copy.size();
- DEBUG_MSG("particle_set::init_as_subset_of, desired subset_size = %ld, src_local_number = %ld, src_total_number = %ld\n",
- this->total_number_of_particles,
- src.getLocalNbParticles(),
- src.getGlobalNbParticles());
+ //DEBUG_MSG("particle_set::init_as_subset_of, desired subset_size = %ld, src_local_number = %ld, src_total_number = %ld\n",
+ // this->total_number_of_particles,
+ // src.getLocalNbParticles(),
+ // src.getGlobalNbParticles());
// dumb selection of interesting particles
for (partsize_t ii = 0; ii < partsize_t(src.getLocalNbParticles()); ii++)
@@ -493,6 +493,13 @@ class particle_set: public abstract_particle_set
delete[] tmp_local_state;
delete[] tmp_local_indices;
delete[] tmp_local_labels;
+
+ std::vector<hsize_t> tmp_file_layout(1);
+ tmp_file_layout[0] = hsize_t(this->total_number_of_particles);
+
+ int file_layout_result = this->setParticleFileLayout(tmp_file_layout);
+ variable_used_only_in_assert(file_layout_result);
+ assert(file_layout_result == EXIT_SUCCESS);
return EXIT_SUCCESS;
}
@@ -560,6 +567,10 @@ class particle_set: public abstract_particle_set
delete[] tmp_local_state;
delete[] tmp_local_indices;
delete[] tmp_local_labels;
+
+ int file_layout_result = this->setParticleFileLayout(src->getParticleFileLayout());
+ variable_used_only_in_assert(file_layout_result);
+ assert(file_layout_result == EXIT_SUCCESS);
return EXIT_SUCCESS;
}
@@ -582,6 +593,8 @@ class particle_set: public abstract_particle_set
* which would need multiple passes of the data to figure out the new indices.
* It's easier to just clean up the data before feeding it to this method.
*
+ * \warning `indices_of_interest` must be in increasing order.
+ *
*/
int init_as_subset_of(
abstract_particle_set &src,
@@ -591,6 +604,7 @@ class particle_set: public abstract_particle_set
TIMEZONE("particle_set::init_as_subset_of version 2");
assert(indices_of_interest.size() > 0);
assert(partsize_t(indices_of_interest.size()) < src.getTotalNumberOfParticles());
+ // ensure that indices of interest are sorted
particle_rnumber *tmp_local_state = new particle_rnumber[state_size * src.getLocalNumberOfParticles()];
partsize_t *tmp_local_indices = new partsize_t[src.getLocalNumberOfParticles()];
partsize_t *tmp_local_labels = new partsize_t[src.getLocalNumberOfParticles()];
@@ -639,13 +653,13 @@ class particle_set: public abstract_particle_set
this->local_number_of_particles++;
}
}
- DEBUG_MSG("particle_set::init_as_subset_of --- after particle selection local number of particles is %d\n",
- this->local_number_of_particles);
- for (partsize_t ii = 0; ii < this->local_number_of_particles; ii++)
- DEBUG_MSG(" --- ii = %d, tmp_local_labels[ii] = %d, tmp_local_indices[ii] = %d\n",
- ii,
- tmp_local_labels[ii],
- tmp_local_indices[ii]);
+ //DEBUG_MSG("particle_set::init_as_subset_of --- after particle selection local number of particles is %d\n",
+ // this->local_number_of_particles);
+ //for (partsize_t ii = 0; ii < this->local_number_of_particles; ii++)
+ // DEBUG_MSG(" --- ii = %d, tmp_local_labels[ii] = %d, tmp_local_indices[ii] = %d\n",
+ // ii,
+ // tmp_local_labels[ii],
+ // tmp_local_indices[ii]);
// now we actually put the data "here"
if (this->local_number_of_particles > 0)
@@ -683,6 +697,14 @@ class particle_set: public abstract_particle_set
delete[] tmp_local_state;
delete[] tmp_local_indices;
delete[] tmp_local_labels;
+
+ // update particle file layout
+ std::vector<hsize_t> tmp_file_layout(1);
+ tmp_file_layout[0] = hsize_t(this->total_number_of_particles);
+
+ int file_layout_result = this->setParticleFileLayout(tmp_file_layout);
+ variable_used_only_in_assert(file_layout_result);
+ assert(file_layout_result == EXIT_SUCCESS);
return EXIT_SUCCESS;
}
--
GitLab
From b242e3d6d1ba965ebff8cde75ef2accc875f6466 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sun, 12 Dec 2021 17:53:01 +0100
Subject: [PATCH 22/37] adds comment
---
TurTLE/test/test_particle_deleter.py | 2 ++
1 file changed, 2 insertions(+)
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index c9a0e483..590aa6cb 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -328,6 +328,8 @@ class ADNS(TurTLE.DNS):
def main():
bla = ADNS()
bla.launch(sys.argv[1:])
+ # TODO
+ # here we need to replicate deletion functionality, confirm resulting particles are same as c++
return None
--
GitLab
From f7664fba707f3c5b6047be6ca21b6ef544ec0fac Mon Sep 17 00:00:00 2001
From: Tobias Baetge <tobias.baetge@ds.mpg.de>
Date: Mon, 13 Dec 2021 18:13:05 +0100
Subject: [PATCH 23/37] adding comparison to python deleting replication
---
TurTLE/test/test_particle_deleter.py | 33 +++++++++++++++++++++++++---
1 file changed, 30 insertions(+), 3 deletions(-)
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index 590aa6cb..ce58bb4e 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -326,12 +326,39 @@ class ADNS(TurTLE.DNS):
return None
def main():
+ #python routine replicating deleting pattern of the test in c++
+ def delete_index(before_py, iteration):
+ index2 = (iteration*3)%len(before_py)
+ index_to_delete = np.array([1, index2])
+ if index2==1:
+ index_to_delete = np.array([1, 3])
+ after_py = np.delete(before_py, index_to_delete)
+ return before_py, after_py
+
+ #running c++ test
bla = ADNS()
bla.launch(sys.argv[1:])
- # TODO
- # here we need to replicate deletion functionality, confirm resulting particles are same as c++
+ #read parameter file
+ f = h5py.File('tracer_set_particle_deleter.h5', 'r')
+ p = f['parameters']
+ #initialising list of indeces
+ before_py = np.arange(0,np.array(p['nparticles']))
+ #comparing list of particles for each iteration
+ for i in range(np.array(p['niter_todo'])):
+ #apply python routine
+ before_py, after_py = delete_index(before_py, i)
+ #load turtle indeces
+ before_turtle = h5py.File('particle_index_before.h5','r')
+ after_turtle = h5py.File('particle_index_after.h5','r')
+ index_before_turtle = np.array(before_turtle['tracers0/index/'+str(i)])
+ index_after_turtle = np.array(after_turtle['tracers0/index/'+str(i)])
+ index_before_turtle = index_before_turtle.flatten()
+ index_after_turtle = index_after_turtle.flatten()
+ assert(np.all(before_py==index_before_turtle) and np.all(after_py==index_after_turtle))
+ if(len(after_py)==2):
+ break
+ before_py = after_py
return None
-
if __name__ == '__main__':
main()
--
GitLab
From 9b9bfb90ee67c578ef84ae728854a864ca9649d5 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 09:45:23 +0100
Subject: [PATCH 24/37] updates ci config
---
.gitlab-ci.yml | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 8fbf29c0..925d1548 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -7,7 +7,7 @@ image: gitlab-registry.mpcdf.mpg.de/mpcdf/module-image
.load_modules: &load_modules |
module purge
- module load cmake $COMPILER $MPI gsl hdf5-mpi fftw-mpi anaconda/3
+ module load cmake $COMPILER $MPI gsl hdf5-mpi/1.12.1 fftw-mpi/3.3.10 anaconda/3/2021.05
export FFTW_DIR=$FFTW_HOME
.build: &build |
@@ -18,8 +18,8 @@ image: gitlab-registry.mpcdf.mpg.de/mpcdf/module-image
.run_tests: &run_tests |
cd build
- mkdir -p /usr/local/lib/python3.7/site-packages
- export PYTHONPATH=/usr/local/lib/python3.7/site-packages:${PYTHONPATH}
+ mkdir -p /usr/local/lib/python3.8/site-packages
+ export PYTHONPATH=/usr/local/lib/python3.8/site-packages:${PYTHONPATH}
make VERBOSE=1 install
export CMAKE_PREFIX_PATH=/usr/local/lib:${CMAKE_PREFIX_PATH}
env CTEST_OUTPUT_ON_FAILURE=1 make test
@@ -30,7 +30,7 @@ build-gcc:
- *load_modules
- *build
variables:
- COMPILER: "gcc"
+ COMPILER: "gcc/10"
MPI: "impi"
MPICXX: "mpigxx"
tags:
@@ -66,7 +66,7 @@ test-gcc:
tags:
- docker
variables:
- COMPILER: "gcc"
+ COMPILER: "gcc/10"
MPI: "impi"
MPICXX: "mpigxx"
needs:
--
GitLab
From 7ec5babb32b6037211a8af7309f0dc7c546583b8 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 09:54:13 +0100
Subject: [PATCH 25/37] reverts a couple of ci changes
---
.gitlab-ci.yml | 6 +++---
1 file changed, 3 insertions(+), 3 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 925d1548..21733373 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -7,7 +7,7 @@ image: gitlab-registry.mpcdf.mpg.de/mpcdf/module-image
.load_modules: &load_modules |
module purge
- module load cmake $COMPILER $MPI gsl hdf5-mpi/1.12.1 fftw-mpi/3.3.10 anaconda/3/2021.05
+ module load cmake $COMPILER $MPI gsl hdf5-mpi/1.12.1 fftw-mpi anaconda
export FFTW_DIR=$FFTW_HOME
.build: &build |
@@ -18,8 +18,8 @@ image: gitlab-registry.mpcdf.mpg.de/mpcdf/module-image
.run_tests: &run_tests |
cd build
- mkdir -p /usr/local/lib/python3.8/site-packages
- export PYTHONPATH=/usr/local/lib/python3.8/site-packages:${PYTHONPATH}
+ mkdir -p /usr/local/lib/python3.7/site-packages
+ export PYTHONPATH=/usr/local/lib/python3.7/site-packages:${PYTHONPATH}
make VERBOSE=1 install
export CMAKE_PREFIX_PATH=/usr/local/lib:${CMAKE_PREFIX_PATH}
env CTEST_OUTPUT_ON_FAILURE=1 make test
--
GitLab
From e1b85e03a7c759d902693b6488a86e424ccfd674 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 10:15:21 +0100
Subject: [PATCH 26/37] bumps intel version for CI
---
.gitlab-ci.yml | 8 ++++----
1 file changed, 4 insertions(+), 4 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 21733373..2a40109e 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -47,8 +47,8 @@ build-intel:
- *load_modules
- *build
variables:
- COMPILER: "intel"
- MPI: "impi"
+ COMPILER: "intel/21.3.0"
+ MPI: "impi/2021.3"
MPICXX: "mpiicpc"
tags:
- docker
@@ -81,8 +81,8 @@ test-intel:
tags:
- docker
variables:
- COMPILER: "intel"
- MPI: "impi"
+ COMPILER: "intel/21.3.0"
+ MPI: "impi/2021.3"
MPICXX: "mpiicpc"
needs:
- job: build-intel
--
GitLab
From 351e2c9b8753604580ba376511a434b81002dacd Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 10:52:23 +0100
Subject: [PATCH 27/37] reverts intel CI changes
---
.gitlab-ci.yml | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 2a40109e..b23350e2 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -7,7 +7,7 @@ image: gitlab-registry.mpcdf.mpg.de/mpcdf/module-image
.load_modules: &load_modules |
module purge
- module load cmake $COMPILER $MPI gsl hdf5-mpi/1.12.1 fftw-mpi anaconda
+ module load cmake ${COMPILER} ${MPI} gsl hdf5-mpi/1.12.1 fftw-mpi anaconda
export FFTW_DIR=$FFTW_HOME
.build: &build |
@@ -47,8 +47,8 @@ build-intel:
- *load_modules
- *build
variables:
- COMPILER: "intel/21.3.0"
- MPI: "impi/2021.3"
+ COMPILER: "intel"
+ MPI: "impi"
MPICXX: "mpiicpc"
tags:
- docker
@@ -81,8 +81,8 @@ test-intel:
tags:
- docker
variables:
- COMPILER: "intel/21.3.0"
- MPI: "impi/2021.3"
+ COMPILER: "intel"
+ MPI: "impi"
MPICXX: "mpiicpc"
needs:
- job: build-intel
--
GitLab
From b0c876b3ae02caab082c555003a311d408ed3fa0 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 15:38:29 +0100
Subject: [PATCH 28/37] [wip] switches intel module for CI
---
.gitlab-ci.yml | 16 ++++++++++++----
1 file changed, 12 insertions(+), 4 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index abc11a69..ba80e058 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -37,6 +37,8 @@ build-gcc:
script:
- *load_modules
- *export_GCC_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
COMPILER: "gcc/10"
@@ -55,10 +57,12 @@ build-intel:
script:
- *load_modules
- *export_INTEL_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
- COMPILER: "intel"
- MPI: "impi"
+ COMPILER: "intel/21.3.0"
+ MPI: "impi/2021.3"
MPICXX: "mpiicpc"
tags:
- docker
@@ -73,6 +77,8 @@ test-gcc:
script:
- *load_modules
- *export_GCC_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *run_tests
tags:
- docker
@@ -89,12 +95,14 @@ test-intel:
script:
- *load_modules
- *export_INTEL_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *run_tests
tags:
- docker
variables:
- COMPILER: "intel"
- MPI: "impi"
+ COMPILER: "intel/21.3.0"
+ MPI: "impi/2021.3"
MPICXX: "mpiicpc"
needs:
- job: build-intel
--
GitLab
From 8f4503d2673f71c05e69bed09df702b56bf7a4da Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 17:01:29 +0100
Subject: [PATCH 29/37] reverts intel CI changes
---
.gitlab-ci.yml | 12 ++++++------
1 file changed, 6 insertions(+), 6 deletions(-)
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index ba80e058..bf0e4b0b 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -41,7 +41,7 @@ build-gcc:
export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
- COMPILER: "gcc/10"
+ COMPILER: "gcc"
MPI: "impi"
MPICXX: "mpigxx"
tags:
@@ -61,8 +61,8 @@ build-intel:
export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
- COMPILER: "intel/21.3.0"
- MPI: "impi/2021.3"
+ COMPILER: "intel"
+ MPI: "impi"
MPICXX: "mpiicpc"
tags:
- docker
@@ -83,7 +83,7 @@ test-gcc:
tags:
- docker
variables:
- COMPILER: "gcc/10"
+ COMPILER: "gcc"
MPI: "impi"
MPICXX: "mpigxx"
needs:
@@ -101,8 +101,8 @@ test-intel:
tags:
- docker
variables:
- COMPILER: "intel/21.3.0"
- MPI: "impi/2021.3"
+ COMPILER: "intel"
+ MPI: "impi"
MPICXX: "mpiicpc"
needs:
- job: build-intel
--
GitLab
From a6a1d3d91cdef39a1853a30a0526f0c9f6aa8531 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 17:53:53 +0100
Subject: [PATCH 30/37] adds Barrier after H5Fclose
---
cpp/full_code/joint_acc_vel_stats.cpp | 4 ++++
1 file changed, 4 insertions(+)
diff --git a/cpp/full_code/joint_acc_vel_stats.cpp b/cpp/full_code/joint_acc_vel_stats.cpp
index f0f500dd..ae4ddda0 100644
--- a/cpp/full_code/joint_acc_vel_stats.cpp
+++ b/cpp/full_code/joint_acc_vel_stats.cpp
@@ -62,6 +62,8 @@ int joint_acc_vel_stats<rnumber>::initialize(void)
else
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -76,6 +78,8 @@ int joint_acc_vel_stats<rnumber>::initialize(void)
H5Dread(dset, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &this->max_velocity_estimate);
H5Dclose(dset);
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
if (this->myrank == 0)
{
// set caching parameters
--
GitLab
From 582749fa5d8b73dd43b1315c39579b9dab3fc18b Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Wed, 5 Jan 2022 18:13:15 +0100
Subject: [PATCH 31/37] adds calls to MPI_barrier after H5Fclose
---
cpp/full_code/bandpass_stats.cpp | 4 ++++
cpp/full_code/field_single_to_double.cpp | 4 ++++
cpp/full_code/get_3D_correlations.cpp | 4 ++++
cpp/full_code/get_rfields.cpp | 4 ++++
cpp/full_code/get_velocity.cpp | 4 ++++
cpp/full_code/pressure_stats.cpp | 4 ++++
cpp/full_code/resize.cpp | 4 ++++
cpp/full_code/write_rpressure.cpp | 2 ++
8 files changed, 30 insertions(+)
diff --git a/cpp/full_code/bandpass_stats.cpp b/cpp/full_code/bandpass_stats.cpp
index b176f892..55bd5301 100644
--- a/cpp/full_code/bandpass_stats.cpp
+++ b/cpp/full_code/bandpass_stats.cpp
@@ -47,6 +47,8 @@ int bandpass_stats<rnumber>::initialize(void)
if (this->checkpoints_per_file == INT_MAX) // value returned if dataset does not exist
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -67,6 +69,8 @@ int bandpass_stats<rnumber>::initialize(void)
parameter_file,
"/bandpass_stats/parameters/filter_type");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
assert(this->k0list.size() == this->k1list.size());
if (this->myrank == 0)
{
diff --git a/cpp/full_code/field_single_to_double.cpp b/cpp/full_code/field_single_to_double.cpp
index 37eb482f..f494216e 100644
--- a/cpp/full_code/field_single_to_double.cpp
+++ b/cpp/full_code/field_single_to_double.cpp
@@ -58,6 +58,8 @@ int field_single_to_double<rnumber>::initialize(void)
else
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -67,6 +69,8 @@ int field_single_to_double<rnumber>::initialize(void)
parameter_file,
"/field_single_to_double/parameters/iteration_list");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/get_3D_correlations.cpp b/cpp/full_code/get_3D_correlations.cpp
index 9a672f1c..ce15bfb2 100644
--- a/cpp/full_code/get_3D_correlations.cpp
+++ b/cpp/full_code/get_3D_correlations.cpp
@@ -62,6 +62,8 @@ int get_3D_correlations<rnumber>::initialize(void)
if (this->checkpoints_per_file == INT_MAX) // value returned if dataset does not exist
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -70,6 +72,8 @@ int get_3D_correlations<rnumber>::initialize(void)
parameter_file,
"/get_3D_correlations/parameters/iteration_list");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/get_rfields.cpp b/cpp/full_code/get_rfields.cpp
index c3cf6c9b..55f85e6e 100644
--- a/cpp/full_code/get_rfields.cpp
+++ b/cpp/full_code/get_rfields.cpp
@@ -65,6 +65,8 @@ int get_rfields<rnumber>::initialize(void)
if (this->checkpoints_per_file == INT_MAX) // value returned if dataset does not exist
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -76,6 +78,8 @@ int get_rfields<rnumber>::initialize(void)
parameter_file,
"/get_rfields/parameters/TrS2_on");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/get_velocity.cpp b/cpp/full_code/get_velocity.cpp
index 1509854d..3f4ae193 100644
--- a/cpp/full_code/get_velocity.cpp
+++ b/cpp/full_code/get_velocity.cpp
@@ -53,6 +53,8 @@ int get_velocity<rnumber>::initialize(void)
if (this->checkpoints_per_file == INT_MAX) // value returned if dataset does not exist
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -61,6 +63,8 @@ int get_velocity<rnumber>::initialize(void)
parameter_file,
"/get_velocity/parameters/iteration_list");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/pressure_stats.cpp b/cpp/full_code/pressure_stats.cpp
index fc060b4d..2abb333e 100644
--- a/cpp/full_code/pressure_stats.cpp
+++ b/cpp/full_code/pressure_stats.cpp
@@ -83,6 +83,8 @@ int pressure_stats<rnumber>::initialize(void)
if (this->checkpoints_per_file == INT_MAX) // value returned if dataset does not exist
this->checkpoints_per_file = 1;
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -94,6 +96,8 @@ int pressure_stats<rnumber>::initialize(void)
parameter_file,
"/pressure_stats/parameters/max_pressure_estimate");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
if (this->myrank == 0)
{
// set caching parameters
diff --git a/cpp/full_code/resize.cpp b/cpp/full_code/resize.cpp
index f074a46c..e81601ef 100644
--- a/cpp/full_code/resize.cpp
+++ b/cpp/full_code/resize.cpp
@@ -43,6 +43,8 @@ int resize<rnumber>::initialize(void)
this->niter_out = hdf5_tools::read_value<int>(
parameter_file, "/parameters/niter_out");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
parameter_file = H5Fopen(
(this->simname + std::string("_post.h5")).c_str(),
H5F_ACC_RDONLY,
@@ -61,6 +63,8 @@ int resize<rnumber>::initialize(void)
this->new_simname = hdf5_tools::read_string(
parameter_file, "/resize/parameters/new_simname");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
this->new_field = new field<rnumber, FFTW, THREE>(
this->new_nx, this->new_ny, this->new_nz,
diff --git a/cpp/full_code/write_rpressure.cpp b/cpp/full_code/write_rpressure.cpp
index 5835ee33..f0972a7f 100644
--- a/cpp/full_code/write_rpressure.cpp
+++ b/cpp/full_code/write_rpressure.cpp
@@ -20,6 +20,8 @@ int write_rpressure<rnumber>::initialize(void)
parameter_file,
"/write_rpressure/parameters/iteration_list");
H5Fclose(parameter_file);
+ // the following ensures the file is free for rank 0 to open in read/write mode
+ MPI_Barrier(this->comm);
if (this->myrank==0) DEBUG_MSG("Iteration list[0] = %d \n", this->iteration_list[0]);
//get the pressure from here
--
GitLab
From 2934dbea144ed60ea89ba944b2e846e6501e33ca Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Thu, 6 Jan 2022 16:13:50 +0100
Subject: [PATCH 32/37] fixes label redistribution
---
TurTLE/test/particle_set/NSVEparticle_set.cpp | 7 +-
TurTLE/test/test_particle_deleter.py | 25 ++++--
cpp/field.cpp | 2 +
.../interpolation/abstract_particle_set.hpp | 3 +-
cpp/particles/interpolation/particle_set.hpp | 3 +-
cpp/particles/p2p/p2p_distr_mpi.hpp | 57 ++++++++++++-
cpp/particles/particles_distr_mpi.hpp | 85 ++++++++++++++++++-
7 files changed, 167 insertions(+), 15 deletions(-)
diff --git a/TurTLE/test/particle_set/NSVEparticle_set.cpp b/TurTLE/test/particle_set/NSVEparticle_set.cpp
index c6e8ba81..ee89c2e3 100644
--- a/TurTLE/test/particle_set/NSVEparticle_set.cpp
+++ b/TurTLE/test/particle_set/NSVEparticle_set.cpp
@@ -167,15 +167,20 @@ int NSVEparticle_set<rnumber>::do_stats()
"position",
psolver->getIteration());
+ MPI_Barrier(this->comm);
+
// sample velocity
// first ensure velocity field is computed, and is in real space
+ DEBUG_MSG("test 00\n");
if (!(this->iteration % this->niter_stat == 0))
{
// we need to compute velocity field manually, because it didn't happen in NSVE::do_stats()
this->fs->compute_velocity(this->fs->cvorticity);
- *this->tmp_vec_field = this->fs->cvelocity->get_cdata();
+ *this->tmp_vec_field = *(this->fs->cvelocity);
this->tmp_vec_field->ift();
}
+ DEBUG_MSG("test 01\n");
+ MPI_Barrier(this->comm);
this->pset->writeSample(
this->tmp_vec_field,
this->particles_sample_writer_mpi,
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index ce58bb4e..15acd35b 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -88,7 +88,7 @@ class ADNS(TurTLE.DNS):
self.parameters['niter_part'] = int(1)
self.parameters['niter_part_fine_period'] = int(2)
self.parameters['niter_part_fine_duration'] = int(0)
- self.parameters['nparticles'] = int(10)
+ self.parameters['nparticles'] = int(100)
self.parameters['tracers0_integration_steps'] = int(0)
self.parameters['tracers0_neighbours'] = int(3)
self.parameters['tracers0_smoothness'] = int(2)
@@ -326,7 +326,14 @@ class ADNS(TurTLE.DNS):
return None
def main():
- #python routine replicating deleting pattern of the test in c++
+ #####################################################################
+ ## run c++ test
+ bla = ADNS()
+ bla.launch(sys.argv[1:])
+ #####################################################################
+
+ #####################################################################
+ ## python routine replicating deleting pattern of the test in c++
def delete_index(before_py, iteration):
index2 = (iteration*3)%len(before_py)
index_to_delete = np.array([1, index2])
@@ -334,11 +341,10 @@ def main():
index_to_delete = np.array([1, 3])
after_py = np.delete(before_py, index_to_delete)
return before_py, after_py
-
- #running c++ test
- bla = ADNS()
- bla.launch(sys.argv[1:])
- #read parameter file
+ #####################################################################
+
+ #####################################################################
+ ## perform sanity check
f = h5py.File('tracer_set_particle_deleter.h5', 'r')
p = f['parameters']
#initialising list of indeces
@@ -354,10 +360,13 @@ def main():
index_after_turtle = np.array(after_turtle['tracers0/index/'+str(i)])
index_before_turtle = index_before_turtle.flatten()
index_after_turtle = index_after_turtle.flatten()
- assert(np.all(before_py==index_before_turtle) and np.all(after_py==index_after_turtle))
+ print('testing labels iteration {0}'.format(i))
+ assert(np.all(before_py==index_before_turtle) and
+ np.all(after_py==index_after_turtle))
if(len(after_py)==2):
break
before_py = after_py
+ #####################################################################
return None
if __name__ == '__main__':
diff --git a/cpp/field.cpp b/cpp/field.cpp
index cead918e..6e950ad4 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -1604,6 +1604,7 @@ template <typename rnumber,
void field<rnumber, be, fc>::symmetrize_alternate()
{
TIMEZONE("field::symmetrize");
+ MPI_Barrier(this->clayout->comm); // TODO: figure out if this can be taken out by careful MPI tag generation
start_mpi_profiling_zone(turtle_mpi_pcontrol::FIELD);
assert(!this->real_space_representation);
typename fftw_interface<rnumber>::complex *cdata = this->get_cdata();
@@ -1834,6 +1835,7 @@ void field<rnumber, be, fc>::symmetrize_alternate()
}
}
finish_mpi_profiling_zone(turtle_mpi_pcontrol::FIELD);
+ MPI_Barrier(this->clayout->comm); // TODO: figure out if this can be taken out by careful MPI tag generation
}
/** \brief Enforce Hermitian symmetry (fast and reasonable)
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index f5481cd1..1227f066 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -185,7 +185,8 @@ class abstract_particle_set
this->getParticleState(),
additional_data.data(),
int(additional_data.size()),
- this->getParticleIndices());
+ this->getParticleIndices(),
+ this->getParticleLabels());
return EXIT_SUCCESS;
}
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index dff01348..a530cc14 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -383,7 +383,8 @@ class particle_set: public abstract_particle_set
&this->local_state,
additional_data.data(),
int(additional_data.size()),
- &this->local_indices);
+ &this->local_indices,
+ &this->local_labels);
return EXIT_SUCCESS;
}
diff --git a/cpp/particles/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index afd0c195..3051bcee 100644
--- a/cpp/particles/p2p/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p/p2p_distr_mpi.hpp
@@ -50,6 +50,7 @@ protected:
TAG_NB_PARTICLES,
TAG_POSITION_PARTICLES,
TAG_INDEX_PARTICLES,
+ TAG_LABEL_PARTICLES,
TAG_RESULT_PARTICLES,
};
@@ -64,6 +65,7 @@ protected:
std::unique_ptr<real_number[]> toCompute;
std::unique_ptr<real_number[]> results;
std::unique_ptr<partsize_t[]> indexes;
+ std::unique_ptr<partsize_t[]> labels;
};
enum Action{
@@ -288,7 +290,8 @@ public:
const partsize_t current_my_nb_particles_per_partition[],
real_number particles_positions[],
std::unique_ptr<real_number[]> particles_current_rhs[], const int nb_rhs,
- partsize_t inout_index_particles[]){
+ partsize_t inout_index_particles[],
+ partsize_t inout_label_particles[] = nullptr){
TIMEZONE("p2p_distr_mpi::compute_distr");
// Some processes might not be involved
@@ -297,6 +300,8 @@ public:
return;
}
+ const bool labels_present = (inout_label_particles != nullptr);
+
const long int my_top_z_cell_level = last_cell_level_proc(my_rank);
assert(my_top_z_cell_level < nb_cell_levels[IDXC_Z]);
const long int my_down_z_cell_level = first_cell_level_proc(my_rank);
@@ -369,6 +374,13 @@ public:
part_to_sort.data(),
inout_index_particles,
&buffer);
+ if (labels_present)
+ permute_copy<partsize_t, 1>(
+ current_offset_particles_for_partition[idxPartition],
+ current_my_nb_particles_per_partition[idxPartition],
+ part_to_sort.data(),
+ inout_label_particles,
+ &buffer);
permute_copy<long int, 1>(
current_offset_particles_for_partition[idxPartition],
current_my_nb_particles_per_partition[idxPartition],
@@ -554,6 +566,21 @@ public:
current_com,
&mpiRequests.back()));
+ if (labels_present)
+ {
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(descriptor.nbParticlesToExchange*size_particle_positions < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Isend(
+ const_cast<partsize_t*>(&inout_label_particles[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]]),
+ int(descriptor.nbParticlesToExchange),
+ particles_utils::GetMpiType(partsize_t()),
+ descriptor.destProc,
+ TAG_LABEL_PARTICLES,
+ current_com,
+ &mpiRequests.back()));
+ }
+
assert(descriptor.toRecvAndMerge == nullptr);
descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
@@ -680,6 +707,8 @@ public:
assert(NbParticlesToReceive != -1);
assert(descriptor.toCompute == nullptr);
assert(descriptor.indexes == nullptr);
+ if (labels_present)
+ assert(descriptor.labels == nullptr);
if(NbParticlesToReceive){
descriptor.toCompute.reset(new real_number[NbParticlesToReceive*size_particle_positions]);
@@ -707,6 +736,22 @@ public:
TAG_INDEX_PARTICLES,
current_com,
&mpiRequests.back()));
+
+ if (labels_present)
+ {
+ descriptor.labels.reset(new partsize_t[NbParticlesToReceive]);
+ whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
+ mpiRequests.emplace_back();
+ assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(
+ descriptor.labels.get(),
+ int(NbParticlesToReceive),
+ particles_utils::GetMpiType(partsize_t()),
+ destProc,
+ TAG_LABEL_PARTICLES,
+ current_com,
+ &mpiRequests.back()));
+ }
}
}
@@ -716,9 +761,12 @@ public:
if(releasedAction.first == COMPUTE_PARTICLES){
NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
descriptor.nbReceived += 1;
- assert(descriptor.nbReceived <= 2);
- if(descriptor.nbReceived == 2){
+ int nbReceived_desired_value = 2;
+ if (labels_present)
+ nbReceived_desired_value = 3;
+
+ if(descriptor.nbReceived == nbReceived_desired_value){
TIMEZONE("compute-particles");
NeighborDescriptor& descriptor = neigDescriptors[releasedAction.second];
assert(descriptor.isRecv);
@@ -726,6 +774,8 @@ public:
assert(descriptor.toCompute != nullptr);
assert(descriptor.indexes != nullptr);
+ if (labels_present)
+ assert(descriptor.labels != nullptr);
// allocate rhs buffer
descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
// clean up rhs buffer
@@ -823,6 +873,7 @@ public:
&mpiRequests.back()));
delete[] descriptor.toCompute.release();
delete[] descriptor.indexes.release();
+ delete[] descriptor.labels.release();
}
}
//////////////////////////////////////////////////////////////////////
diff --git a/cpp/particles/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index 1681ea3d..c4dac2ad 100644
--- a/cpp/particles/particles_distr_mpi.hpp
+++ b/cpp/particles/particles_distr_mpi.hpp
@@ -60,6 +60,9 @@ protected:
TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
+ TAG_LOW_UP_MOVED_PARTICLES_LABELS,
+ TAG_UP_LOW_MOVED_PARTICLES_LABELS,
+
TAG_LOW_UP_MOVED_PARTICLES_RHS,
TAG_LOW_UP_MOVED_PARTICLES_RHS_MAX = TAG_LOW_UP_MOVED_PARTICLES_RHS+MaxNbRhs,
@@ -557,7 +560,8 @@ public:
partsize_t* nb_particles,
std::unique_ptr<real_number[]>* inout_positions_particles,
std::unique_ptr<real_number[]> inout_rhs_particles[], const int in_nb_rhs,
- std::unique_ptr<partsize_t[]>* inout_index_particles){
+ std::unique_ptr<partsize_t[]>* inout_index_particles,
+ std::unique_ptr<partsize_t[]>* inout_label_particles = nullptr){
TIMEZONE("redistribute");
// Some latest processes might not be involved
@@ -565,6 +569,8 @@ public:
return;
}
+ const bool labels_present = (inout_label_particles != nullptr);
+
{// TODO remove
partsize_t partOffset = 0;
for(int idxPartition = 0 ; idxPartition < current_partition_size ; ++idxPartition){
@@ -603,6 +609,15 @@ public:
(*inout_index_particles)[size_particle_index*idx2+idx_val]);
}
+ if (labels_present)
+ {
+ for(int idx_val = 0 ; idx_val < size_particle_index ; ++idx_val)
+ {
+ std::swap((*inout_label_particles)[size_particle_index*idx1+idx_val],
+ (*inout_label_particles)[size_particle_index*idx2+idx_val]);
+ }
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
std::swap(inout_rhs_particles[idx_rhs][idx1*size_particle_rhs + idx_val],
@@ -629,6 +644,15 @@ public:
(*inout_index_particles)[size_particle_index*idx2+idx_val]);
}
+ if (labels_present)
+ {
+ for(int idx_val = 0 ; idx_val < size_particle_index ; ++idx_val)
+ {
+ std::swap((*inout_label_particles)[size_particle_index*idx1+idx_val],
+ (*inout_label_particles)[size_particle_index*idx2+idx_val]);
+ }
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
std::swap(inout_rhs_particles[idx_rhs][idx1*size_particle_rhs + idx_val],
@@ -645,6 +669,8 @@ public:
std::unique_ptr<real_number[]> newParticlesUp;
std::unique_ptr<partsize_t[]> newParticlesLowIndexes;
std::unique_ptr<partsize_t[]> newParticlesUpIndexes;
+ std::unique_ptr<partsize_t[]> newParticlesLowLabels;
+ std::unique_ptr<partsize_t[]> newParticlesUpLabels;
std::vector<std::unique_ptr<real_number[]>> newParticlesLowRhs(in_nb_rhs);
std::vector<std::unique_ptr<real_number[]>> newParticlesUpRhs(in_nb_rhs);
@@ -679,6 +705,13 @@ public:
(my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if (labels_present)
+ {
+ AssertMpi(MPI_Isend(&(*inout_label_particles)[0], int(nbOutLower*size_particle_index), particles_utils::GetMpiType(partsize_t()),
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_LABELS,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
@@ -715,6 +748,12 @@ public:
AssertMpi(MPI_Isend(&(*inout_index_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_index], int(nbOutUpper*size_particle_index),
particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if (labels_present)
+ {
+ AssertMpi(MPI_Isend(&(*inout_label_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_index], int(nbOutUpper*size_particle_index),
+ particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_LABELS,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -758,6 +797,18 @@ public:
(my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if (labels_present)
+ {
+ newParticlesLowLabels.reset(new partsize_t[nbNewFromLow*size_particle_index]);
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(nbNewFromLow*size_particle_index < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(&newParticlesLowLabels[0], int(nbNewFromLow*size_particle_index),
+ particles_utils::GetMpiType(partsize_t()),
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_LABELS,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
newParticlesLowRhs[idx_rhs].reset(new real_number[nbNewFromLow*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -788,6 +839,18 @@ public:
(my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
+ if (labels_present)
+ {
+ newParticlesUpLabels.reset(new partsize_t[nbNewFromUp*size_particle_index]);
+ whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
+ mpiRequests.emplace_back();
+ assert(nbNewFromUp*size_particle_index < std::numeric_limits<int>::max());
+ AssertMpi(MPI_Irecv(&newParticlesUpLabels[0], int(nbNewFromUp*size_particle_index),
+ particles_utils::GetMpiType(partsize_t()),
+ (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_LABELS,
+ MPI_COMM_WORLD, &mpiRequests.back()));
+ }
+
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
newParticlesUpRhs[idx_rhs].reset(new real_number[nbNewFromUp*size_particle_rhs]);
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -818,6 +881,7 @@ public:
std::unique_ptr<real_number[]> newArray(new real_number[myTotalNewNbParticles*size_particle_positions]);
std::unique_ptr<partsize_t[]> newArrayIndexes(new partsize_t[myTotalNewNbParticles*size_particle_index]);
+ std::unique_ptr<partsize_t[]> newArrayLabels(new partsize_t[myTotalNewNbParticles*size_particle_index]);
std::vector<std::unique_ptr<real_number[]>> newArrayRhs(in_nb_rhs);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
newArrayRhs[idx_rhs].reset(new real_number[myTotalNewNbParticles*size_particle_rhs]);
@@ -828,6 +892,8 @@ public:
const particles_utils::fixed_copy fcp(0, 0, nbNewFromLow);
fcp.copy(newArray, newParticlesLow, size_particle_positions);
fcp.copy(newArrayIndexes, newParticlesLowIndexes, size_particle_index);
+ if (labels_present)
+ fcp.copy(newArrayLabels, newParticlesLowLabels, size_particle_index);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
fcp.copy(newArrayRhs[idx_rhs], newParticlesLowRhs[idx_rhs], size_particle_rhs);
}
@@ -838,6 +904,10 @@ public:
const particles_utils::fixed_copy fcp(nbNewFromLow, nbOutLower, nbOldParticlesInside);
fcp.copy(newArray, (*inout_positions_particles), size_particle_positions);
fcp.copy(newArrayIndexes, (*inout_index_particles), size_particle_index);
+ if (labels_present)
+ {
+ fcp.copy(newArrayLabels, (*inout_label_particles), size_particle_index);
+ }
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
fcp.copy(newArrayRhs[idx_rhs], inout_rhs_particles[idx_rhs], size_particle_rhs);
}
@@ -848,6 +918,8 @@ public:
const particles_utils::fixed_copy fcp(nbNewFromLow+nbOldParticlesInside, 0, nbNewFromUp);
fcp.copy(newArray, newParticlesUp, size_particle_positions);
fcp.copy(newArrayIndexes, newParticlesUpIndexes, size_particle_index);
+ if (labels_present)
+ fcp.copy(newArrayLabels, newParticlesUpLabels, size_particle_index);
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
fcp.copy(newArrayRhs[idx_rhs], newParticlesUpRhs[idx_rhs], size_particle_rhs);
}
@@ -855,6 +927,10 @@ public:
(*inout_positions_particles) = std::move(newArray);
(*inout_index_particles) = std::move(newArrayIndexes);
+ if (labels_present)
+ {
+ (*inout_label_particles) = std::move(newArrayLabels);
+ }
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
inout_rhs_particles[idx_rhs] = std::move(newArrayRhs[idx_rhs]);
}
@@ -878,6 +954,13 @@ public:
std::swap((*inout_index_particles)[size_particle_index*idx1 + idx_val],
(*inout_index_particles)[size_particle_index*idx2 + idx_val]);
}
+ if (labels_present)
+ {
+ for(int idx_val = 0 ; idx_val < size_particle_index ; ++idx_val){
+ std::swap((*inout_label_particles)[size_particle_index*idx1 + idx_val],
+ (*inout_label_particles)[size_particle_index*idx2 + idx_val]);
+ }
+ }
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
--
GitLab
From 6b298e1447474e79270efdbcc479ea039ec7d117 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 15 Jan 2022 18:10:17 +0100
Subject: [PATCH 33/37] partial header clean up
I want to ensure that we're not making a very complicated header
inclusion graph, so I untangled "base.hpp" from "hdf5_tools.hpp".
---
cpp/base.hpp | 11 +++++------
cpp/fftw_tools.hpp | 8 ++++----
cpp/field.hpp | 8 ++++----
cpp/field_binary_IO.hpp | 8 ++++----
cpp/field_layout.hpp | 7 ++++---
cpp/full_code/Gauss_field_test.cpp | 1 +
cpp/full_code/NSE.cpp | 1 +
cpp/full_code/NSVE.cpp | 1 +
cpp/full_code/NSVE_field_stats.cpp | 1 +
cpp/full_code/bandpass_stats.cpp | 1 +
cpp/full_code/code_base.cpp | 1 +
cpp/full_code/dealias_test.cpp | 1 +
cpp/full_code/direct_numerical_simulation.hpp | 1 +
cpp/full_code/field_single_to_double.cpp | 1 +
cpp/full_code/field_test.cpp | 1 +
cpp/full_code/filter_test.cpp | 1 +
cpp/full_code/get_3D_correlations.cpp | 2 +-
cpp/full_code/get_rfields.cpp | 1 +
cpp/full_code/get_velocity.cpp | 1 +
cpp/full_code/joint_acc_vel_stats.cpp | 1 +
cpp/full_code/native_binary_to_hdf5.cpp | 1 +
cpp/full_code/pressure_stats.cpp | 1 +
cpp/full_code/resize.cpp | 1 +
cpp/full_code/symmetrize_test.cpp | 1 +
cpp/full_code/write_rpressure.cpp | 3 ++-
cpp/kspace.hpp | 8 ++++----
cpp/particles/particles_input_hdf5.hpp | 1 +
cpp/particles/particles_output_sampling_hdf5.hpp | 2 +-
cpp/spectrum_function.hpp | 6 ++++++
cpp/vorticity_equation.hpp | 9 +++++----
30 files changed, 59 insertions(+), 32 deletions(-)
diff --git a/cpp/base.hpp b/cpp/base.hpp
index be05de04..55562ebc 100644
--- a/cpp/base.hpp
+++ b/cpp/base.hpp
@@ -23,17 +23,16 @@
**********************************************************************/
+#ifndef BASE_HPP
+
+#define BASE_HPP
+
#include <cassert>
#include <mpi.h>
#include <stdarg.h>
#include <iostream>
#include <typeinfo>
-#include "hdf5_tools.hpp"
-
-#ifndef BASE
-
-#define BASE
/////////////////////////////////////////////////////////////
@@ -210,5 +209,5 @@ inline void DEBUG_MSG_WAIT(MPI_Comm communicator, const char * format, ...)
#define variable_used_only_in_assert(x) ((void)(x))
-#endif//BASE
+#endif//BASE_HPP
diff --git a/cpp/fftw_tools.hpp b/cpp/fftw_tools.hpp
index 11c1cbf9..5fbcdf81 100644
--- a/cpp/fftw_tools.hpp
+++ b/cpp/fftw_tools.hpp
@@ -24,16 +24,16 @@
+#ifndef FFTW_TOOLS
+#define FFTW_TOOLS
+
+
#include <mpi.h>
#include <fftw3-mpi.h>
#include <string>
#include <map>
#include <cmath>
-#ifndef FFTW_TOOLS
-
-#define FFTW_TOOLS
-
// helper function
inline void complex_number_phase_shifter(double &xx, double &yy, const double cos_val, const double sin_val)
diff --git a/cpp/field.hpp b/cpp/field.hpp
index 8da33703..f6c69a7d 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -24,6 +24,10 @@
+#ifndef FIELD_HPP
+
+#define FIELD_HPP
+
#include <hdf5.h>
#include <unordered_map>
#include <vector>
@@ -31,10 +35,6 @@
#include "kspace.hpp"
#include "omputils.hpp"
-#ifndef FIELD_HPP
-
-#define FIELD_HPP
-
/** \class field
* \brief Holds field data, performs FFTs and HDF5 I/O operations.
*
diff --git a/cpp/field_binary_IO.hpp b/cpp/field_binary_IO.hpp
index 4983ffe8..c615f9e2 100644
--- a/cpp/field_binary_IO.hpp
+++ b/cpp/field_binary_IO.hpp
@@ -24,6 +24,10 @@
+#ifndef FIELD_BINARY_IO_HPP
+
+#define FIELD_BINARY_IO_HPP
+
#include <vector>
#include <string>
#include "base.hpp"
@@ -31,10 +35,6 @@
#include "field_layout.hpp"
#include "field.hpp"
-#ifndef FIELD_BINARY_IO_HPP
-
-#define FIELD_BINARY_IO_HPP
-
/* could this be a boolean somehow?*/
enum field_representation: bool {
REAL = true,
diff --git a/cpp/field_layout.hpp b/cpp/field_layout.hpp
index 27d0ea9a..e7ba578c 100644
--- a/cpp/field_layout.hpp
+++ b/cpp/field_layout.hpp
@@ -24,13 +24,14 @@
-#include <vector>
-#include "base.hpp"
-
#ifndef FIELD_LAYOUT_HPP
#define FIELD_LAYOUT_HPP
+#include <vector>
+#include <hdf5.h>
+#include "base.hpp"
+
enum field_components {ONE, THREE, THREExTHREE};
constexpr unsigned int ncomp(
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index a9c604ee..6f7cdfd1 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -29,6 +29,7 @@
#include <random>
#include "Gauss_field_test.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
diff --git a/cpp/full_code/NSE.cpp b/cpp/full_code/NSE.cpp
index 7c1c058d..25a4c395 100644
--- a/cpp/full_code/NSE.cpp
+++ b/cpp/full_code/NSE.cpp
@@ -28,6 +28,7 @@
#include "NSE.hpp"
#include "scope_timer.hpp"
#include "fftw_tools.hpp"
+#include "hdf5_tools.hpp"
#include "shared_array.hpp"
diff --git a/cpp/full_code/NSVE.cpp b/cpp/full_code/NSVE.cpp
index 9abb40fa..09ca7f28 100644
--- a/cpp/full_code/NSVE.cpp
+++ b/cpp/full_code/NSVE.cpp
@@ -28,6 +28,7 @@
#include "NSVE.hpp"
#include "scope_timer.hpp"
#include "fftw_tools.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index 3063943d..4ee109b8 100644
--- a/cpp/full_code/NSVE_field_stats.cpp
+++ b/cpp/full_code/NSVE_field_stats.cpp
@@ -28,6 +28,7 @@
#include "NSVE_field_stats.hpp"
#include "fftw_tools.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/bandpass_stats.cpp b/cpp/full_code/bandpass_stats.cpp
index 55bd5301..7985a26a 100644
--- a/cpp/full_code/bandpass_stats.cpp
+++ b/cpp/full_code/bandpass_stats.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "bandpass_stats.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/code_base.cpp b/cpp/full_code/code_base.cpp
index 1015aed9..b62f5015 100644
--- a/cpp/full_code/code_base.cpp
+++ b/cpp/full_code/code_base.cpp
@@ -26,6 +26,7 @@
#include "code_base.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
code_base::code_base(
const MPI_Comm COMMUNICATOR,
diff --git a/cpp/full_code/dealias_test.cpp b/cpp/full_code/dealias_test.cpp
index e3afa2ac..56980b54 100644
--- a/cpp/full_code/dealias_test.cpp
+++ b/cpp/full_code/dealias_test.cpp
@@ -29,6 +29,7 @@
#include <random>
#include "dealias_test.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
diff --git a/cpp/full_code/direct_numerical_simulation.hpp b/cpp/full_code/direct_numerical_simulation.hpp
index 87f19078..be5d4988 100644
--- a/cpp/full_code/direct_numerical_simulation.hpp
+++ b/cpp/full_code/direct_numerical_simulation.hpp
@@ -31,6 +31,7 @@
#include <sys/types.h>
#include <sys/stat.h>
#include "base.hpp"
+#include "hdf5_tools.hpp"
#include "full_code/code_base.hpp"
class direct_numerical_simulation: public code_base
diff --git a/cpp/full_code/field_single_to_double.cpp b/cpp/full_code/field_single_to_double.cpp
index f494216e..9f2082e6 100644
--- a/cpp/full_code/field_single_to_double.cpp
+++ b/cpp/full_code/field_single_to_double.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "field_single_to_double.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/field_test.cpp b/cpp/full_code/field_test.cpp
index 9b4c4460..0aa2b940 100644
--- a/cpp/full_code/field_test.cpp
+++ b/cpp/full_code/field_test.cpp
@@ -28,6 +28,7 @@
#include <random>
#include "field_test.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/filter_test.cpp b/cpp/full_code/filter_test.cpp
index 8d0ae00f..e70cc61c 100644
--- a/cpp/full_code/filter_test.cpp
+++ b/cpp/full_code/filter_test.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "filter_test.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/get_3D_correlations.cpp b/cpp/full_code/get_3D_correlations.cpp
index ce15bfb2..6970e7af 100644
--- a/cpp/full_code/get_3D_correlations.cpp
+++ b/cpp/full_code/get_3D_correlations.cpp
@@ -27,7 +27,7 @@
#include <cmath>
#include "get_3D_correlations.hpp"
#include "scope_timer.hpp"
-
+#include "hdf5_tools.hpp"
template <typename rnumber>
int get_3D_correlations<rnumber>::initialize(void)
diff --git a/cpp/full_code/get_rfields.cpp b/cpp/full_code/get_rfields.cpp
index 55f85e6e..85864ee4 100644
--- a/cpp/full_code/get_rfields.cpp
+++ b/cpp/full_code/get_rfields.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "get_rfields.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/get_velocity.cpp b/cpp/full_code/get_velocity.cpp
index 3f4ae193..a29cc660 100644
--- a/cpp/full_code/get_velocity.cpp
+++ b/cpp/full_code/get_velocity.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "get_velocity.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/joint_acc_vel_stats.cpp b/cpp/full_code/joint_acc_vel_stats.cpp
index ae4ddda0..8fe55d20 100644
--- a/cpp/full_code/joint_acc_vel_stats.cpp
+++ b/cpp/full_code/joint_acc_vel_stats.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "joint_acc_vel_stats.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/native_binary_to_hdf5.cpp b/cpp/full_code/native_binary_to_hdf5.cpp
index a1a3277c..ac0973d6 100644
--- a/cpp/full_code/native_binary_to_hdf5.cpp
+++ b/cpp/full_code/native_binary_to_hdf5.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "native_binary_to_hdf5.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/pressure_stats.cpp b/cpp/full_code/pressure_stats.cpp
index 2abb333e..abeccf5c 100644
--- a/cpp/full_code/pressure_stats.cpp
+++ b/cpp/full_code/pressure_stats.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "pressure_stats.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/resize.cpp b/cpp/full_code/resize.cpp
index e81601ef..12f125e0 100644
--- a/cpp/full_code/resize.cpp
+++ b/cpp/full_code/resize.cpp
@@ -27,6 +27,7 @@
#include <cmath>
#include "resize.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/symmetrize_test.cpp b/cpp/full_code/symmetrize_test.cpp
index da2bef88..904b57ae 100644
--- a/cpp/full_code/symmetrize_test.cpp
+++ b/cpp/full_code/symmetrize_test.cpp
@@ -28,6 +28,7 @@
#include "symmetrize_test.hpp"
#include "fftw_tools.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <typename rnumber>
diff --git a/cpp/full_code/write_rpressure.cpp b/cpp/full_code/write_rpressure.cpp
index f0972a7f..a9e5e4ff 100644
--- a/cpp/full_code/write_rpressure.cpp
+++ b/cpp/full_code/write_rpressure.cpp
@@ -1,8 +1,9 @@
#include <string>
#include <cmath>
+#include <hdf5.h>
#include "write_rpressure.hpp"
#include "scope_timer.hpp"
-
+#include "hdf5_tools.hpp"
template <typename rnumber>
int write_rpressure<rnumber>::initialize(void)
diff --git a/cpp/kspace.hpp b/cpp/kspace.hpp
index ab68a2b0..6457c8bd 100644
--- a/cpp/kspace.hpp
+++ b/cpp/kspace.hpp
@@ -24,6 +24,10 @@
+#ifndef KSPACE_HPP
+
+#define KSPACE_HPP
+
#include <hdf5.h>
#include <vector>
#include <string>
@@ -31,10 +35,6 @@
#include "fftw_interface.hpp"
#include "field_layout.hpp"
-#ifndef KSPACE_HPP
-
-#define KSPACE_HPP
-
enum field_backend {FFTW};
enum kspace_dealias_type {ONE_HALF, TWO_THIRDS, SMOOTH};
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 97fad59d..b1ea662e 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -37,6 +37,7 @@
#include "alltoall_exchanger.hpp"
#include "particles/particles_utils.hpp"
#include "scope_timer.hpp"
+#include "hdf5_tools.hpp"
template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
diff --git a/cpp/particles/particles_output_sampling_hdf5.hpp b/cpp/particles/particles_output_sampling_hdf5.hpp
index 73bc253f..e3cc2e02 100644
--- a/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -27,8 +27,8 @@
#define PARTICLES_OUTPUT_SAMPLING_HDF5_HPP
#include "abstract_particles_output.hpp"
+#include "hdf5_tools.hpp"
-#include <hdf5.h>
template <class partsize_t,
class position_type,
diff --git a/cpp/spectrum_function.hpp b/cpp/spectrum_function.hpp
index 2df2ecab..7be4474c 100644
--- a/cpp/spectrum_function.hpp
+++ b/cpp/spectrum_function.hpp
@@ -1,3 +1,6 @@
+#ifndef SPECTRUM_FUNCTION_CPP
+#define SPECTRUM_FUNCTION_CPP
+
#include<cmath>
#include<vector>
@@ -30,3 +33,6 @@ class spectrum_function
return this->values[index];
}
};
+
+#endif//SPECTRUM_FUNCTION_CPP
+
diff --git a/cpp/vorticity_equation.hpp b/cpp/vorticity_equation.hpp
index 2ba11293..32c6547b 100644
--- a/cpp/vorticity_equation.hpp
+++ b/cpp/vorticity_equation.hpp
@@ -22,6 +22,11 @@
* *
**********************************************************************/
+
+
+#ifndef VORTICITY_EQUATION
+#define VORTICITY_EQUATION
+
#include <sys/stat.h>
#include <stdio.h>
#include <stdlib.h>
@@ -29,10 +34,6 @@
#include "field.hpp"
-#ifndef VORTICITY_EQUATION
-
-#define VORTICITY_EQUATION
-
extern int myrank, nprocs;
--
GitLab
From 804fdd290bb4ea5889972270159eb62359f53479 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 15 Jan 2022 18:20:53 +0100
Subject: [PATCH 34/37] adds optimization to debug mode
---
CMakeLists.txt | 1 +
1 file changed, 1 insertion(+)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b3638881..6c104cca 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -491,6 +491,7 @@ set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}
list(APPEND TURTLE_LIBS "${OpenMP_CXX_LIB_NAMES}")
target_link_libraries(TurTLE ${TURTLE_LIBS})
+target_compile_options(TurTLE PRIVATE $<$<CONFIG:DEBUG>:-O1>)
install(TARGETS TurTLE EXPORT TURTLE_EXPORT DESTINATION lib/ )
install(DIRECTORY ${PROJECT_SOURCE_DIR}/cpp/ DESTINATION include/TurTLE/ FILES_MATCHING PATTERN "*.h*")
--
GitLab
From dbc9a9534334dfaf261b420cb1d5635a5ce3c314 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 15 Jan 2022 18:59:52 +0100
Subject: [PATCH 35/37] adds call to "extract subset" functionality
---
TurTLE/test/particle_set/particle_deleter.cpp | 27 +++++++++++++++++--
TurTLE/test/test_particle_deleter.py | 5 +++-
cpp/hdf5_tools.cpp | 6 +++++
cpp/particles/interpolation/particle_set.hpp | 2 +-
4 files changed, 36 insertions(+), 4 deletions(-)
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
index 165698a9..806d9961 100644
--- a/TurTLE/test/particle_set/particle_deleter.cpp
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -237,13 +237,36 @@ int particle_deleter<rnumber>::do_stats()
"tracers0",
"position/0");
- DEBUG_MSG("before writeParticleIndex");
+ DEBUG_MSG("before writeParticleIndex\n");
this->pset->writeParticleLabels(
"particle_index_after.h5",
"tracers0",
"index",
psolver->getIteration());
- DEBUG_MSG("after writeParticleIndex");
+ DEBUG_MSG("after writeParticleIndex\n");
+ }
+
+ // test particle subset extraction
+ if (this->pset->getTotalNumberOfParticles() > 4)
+ {
+ particle_set<3, 1, 0> temp_pset(
+ this->fs->cvelocity->rlayout,
+ this->dkx,
+ this->dky,
+ this->dkz);
+ std::vector<long long int> indices_to_extract(2);
+ // TODO: possibly use a more complicated pattern, but remember to keep
+ // indices smaller than current total number of particles
+ indices_to_extract[0] = 2;
+ indices_to_extract[1] = 4;
+ temp_pset.init_as_subset_of(*(this->pset), indices_to_extract, false);
+ DEBUG_MSG("before writeParticleLabels for extracted subset\n");
+ temp_pset.writeParticleLabels(
+ "particle_subset_index.h5",
+ "tracers_subset",
+ "index",
+ psolver->getIteration());
+ DEBUG_MSG("after writeParticleLabels for extracted subset\n");
}
return EXIT_SUCCESS;
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index 15acd35b..057abe76 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -312,6 +312,8 @@ class ADNS(TurTLE.DNS):
particle_file.create_group('tracers0/index')
with h5py.File('particle_index_after.h5', 'w') as particle_file:
particle_file.create_group('tracers0/index')
+ with h5py.File('particle_subset_index.h5', 'w') as particle_file:
+ particle_file.create_group('tracers_subset/index')
self.generate_tracer_state(integration_steps = self.parameters['tracers0_integration_steps'],
rseed = opt.particle_rand_seed,
ncomponents = 3)
@@ -345,7 +347,7 @@ def main():
#####################################################################
## perform sanity check
- f = h5py.File('tracer_set_particle_deleter.h5', 'r')
+ f = h5py.File(bla.simname + '.h5', 'r')
p = f['parameters']
#initialising list of indeces
before_py = np.arange(0,np.array(p['nparticles']))
@@ -366,6 +368,7 @@ def main():
if(len(after_py)==2):
break
before_py = after_py
+ ## TODO: check subset extraction as well
#####################################################################
return None
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 920db174..0d0854f9 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -213,6 +213,8 @@ std::vector<number> hdf5_tools::read_vector(
mem_dtype = H5Tcopy(H5T_NATIVE_INT);
else if (typeid(number) == typeid(double))
mem_dtype = H5Tcopy(H5T_NATIVE_DOUBLE);
+ else
+ return result;
dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
dspace = H5Dget_space(dset);
assert(H5Sget_simple_extent_ndims(dspace) == 1);
@@ -240,6 +242,8 @@ number hdf5_tools::read_value(
mem_dtype = H5Tcopy(H5T_NATIVE_LLONG);
else if (typeid(number) == typeid(double))
mem_dtype = H5Tcopy(H5T_NATIVE_DOUBLE);
+ else
+ return result;
if (H5Lexists(group, dset_name.c_str(), H5P_DEFAULT))
{
dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
@@ -276,6 +280,8 @@ int hdf5_tools::write_value_with_single_rank(
mem_dtype = H5Tcopy(H5T_NATIVE_INT);
else if (typeid(number) == typeid(double))
mem_dtype = H5Tcopy(H5T_NATIVE_DOUBLE);
+ else
+ return EXIT_FAILURE;
if (H5Lexists(group, dset_name.c_str(), H5P_DEFAULT))
{
dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index a530cc14..071614ad 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -645,7 +645,7 @@ class particle_set: public abstract_particle_set
if (use_set_complement)
tmp_local_indices[this->local_number_of_particles] -= iii;
else
- tmp_local_indices[this->local_number_of_particles] -= src.getTotalNumberOfParticles() - iii;
+ tmp_local_indices[this->local_number_of_particles] = iii;
// copy particle with label src_label
tmp_local_labels[this->local_number_of_particles] = src_label;
std::copy(src.getParticleState()+state_size*ii,
--
GitLab
From 5c4bf88237908adcf8058320ed889318d7819a86 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Sat, 15 Jan 2022 19:37:36 +0100
Subject: [PATCH 36/37] adds call to Barrier
---
TurTLE/test/particle_set/particle_deleter.cpp | 1 +
1 file changed, 1 insertion(+)
diff --git a/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
index 806d9961..884517cc 100644
--- a/TurTLE/test/particle_set/particle_deleter.cpp
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -286,6 +286,7 @@ int particle_deleter<rnumber>::read_parameters(void)
this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
H5Fclose(parameter_file);
+ MPI_Barrier(this->comm);
return EXIT_SUCCESS;
}
--
GitLab
From 968635b8837e675a26c435e86c6a62204b1e8f3f Mon Sep 17 00:00:00 2001
From: Tobias Baetge <tobias.baetge@ds.mpg.de>
Date: Tue, 18 Jan 2022 10:51:22 +0100
Subject: [PATCH 37/37] add subset test
---
TurTLE/test/test_particle_deleter.py | 7 +++++--
1 file changed, 5 insertions(+), 2 deletions(-)
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
index 057abe76..a7bbe788 100644
--- a/TurTLE/test/test_particle_deleter.py
+++ b/TurTLE/test/test_particle_deleter.py
@@ -147,7 +147,6 @@ class ADNS(TurTLE.DNS):
bla /= np.sum(bla**2, axis = 1)[:, None]**.5
return bla
while nn > 0:
- #TODO introduce initial vel parameter as in established Stokes class, now set to zero
if nn > batch_size:
dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
if dset.shape[1] == 6:
@@ -358,17 +357,21 @@ def main():
#load turtle indeces
before_turtle = h5py.File('particle_index_before.h5','r')
after_turtle = h5py.File('particle_index_after.h5','r')
+ subset_turtle = h5py.File('particle_subset_index.h5', 'r')
index_before_turtle = np.array(before_turtle['tracers0/index/'+str(i)])
index_after_turtle = np.array(after_turtle['tracers0/index/'+str(i)])
+ index_subset = np.array(subset_turtle['tracers_subset/index/'+str(i)])
index_before_turtle = index_before_turtle.flatten()
index_after_turtle = index_after_turtle.flatten()
+ index_subset = index_subset.flatten()
print('testing labels iteration {0}'.format(i))
assert(np.all(before_py==index_before_turtle) and
np.all(after_py==index_after_turtle))
+ print('testing subset labels iteration {0}'.format(i))
+ assert(np.all(after_py[[2, 4]]==index_subset))
if(len(after_py)==2):
break
before_py = after_py
- ## TODO: check subset extraction as well
#####################################################################
return None
--
GitLab