diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 60a3fe604b2e26c55ddfe55e0cf2a8098712927f..bf0e4b0b29a98e95df25ee13be948ce60eeecb18 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/3.22 $COMPILER $MPI gsl hdf5-mpi fftw-mpi anaconda/3
+ module load cmake/3.22 ${COMPILER} ${MPI} gsl hdf5-mpi/1.12.1 fftw-mpi anaconda
export FFTW_DIR=$FFTW_HOME
.build: &build |
@@ -37,6 +37,8 @@ build-gcc:
script:
- *load_modules
- *export_GCC_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
COMPILER: "gcc"
@@ -55,6 +57,8 @@ build-intel:
script:
- *load_modules
- *export_INTEL_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *build
variables:
COMPILER: "intel"
@@ -73,6 +77,8 @@ test-gcc:
script:
- *load_modules
- *export_GCC_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *run_tests
tags:
- docker
@@ -89,6 +95,8 @@ test-intel:
script:
- *load_modules
- *export_INTEL_compilers
+ - >
+ export MPI_HOME=${I_MPI_ROOT}
- *run_tests
tags:
- docker
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d28cce28d684237e08dcbcc6b00c6c0ded1e6049..da1d4bd2924c017dd82018265cd1f04d04a8d97e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -182,6 +182,7 @@ if(TIMING_OUTPUT)
endif()
#####################################################################################
+
set(TURTLE_LIBS "")
#####################################################################################
## HDF5
@@ -488,7 +489,15 @@ 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})
+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*")
@@ -585,6 +594,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/TurTLE/test/particle_set/NSVEparticle_set.cpp b/TurTLE/test/particle_set/NSVEparticle_set.cpp
index 8a32ebc0172282929ca70f6bef1733d0b2f78b0b..ee89c2e36df943fd38dc04372dd2b459d5e21632 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 *
* *
@@ -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"),
@@ -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;
}
@@ -166,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/particle_set/NSVEparticle_set.hpp b/TurTLE/test/particle_set/NSVEparticle_set.hpp
index 64a5beebbcc3feefb06cb7fc1cc6df8dcd578dec..57d7a584fe34d6c9a271a0754712d8cb97657344 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 *
* *
@@ -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/TurTLE/test/particle_set/particle_deleter.cpp b/TurTLE/test/particle_set/particle_deleter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..884517cc8ff891fbc19cf9785153245dd019bc5a
--- /dev/null
+++ b/TurTLE/test/particle_set/particle_deleter.cpp
@@ -0,0 +1,296 @@
+/**********************************************************************
+* *
+* 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_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;
+}
+
+/** \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());
+
+ 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;
+ std::sort(indices_to_delete.begin(), indices_to_delete.end());
+
+ 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,
+ 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\n");
+ this->pset->writeParticleLabels(
+ "particle_index_after.h5",
+ "tracers0",
+ "index",
+ psolver->getIteration());
+ 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;
+}
+
+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);
+ MPI_Barrier(this->comm);
+ 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 0000000000000000000000000000000000000000..0561b56d9cb9c2915ac9a8e133449e50ebd239b3
--- /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/profiler/Kraichnan_scalar_v1.cpp b/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
index 1b33f36fd29111844c09cea1639fbc8b8b809e7e..7e9a647da18ccf034b4a1ff050aedf575ea9abe4 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 cb50290650294cf795fc74011d64fe9925a3fab7..97aa9cdb6e37c49dc6ae446dc3ea88984894eb31 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 fa52286ea563588cc25260be9dff989e0b8e187c..5e112ea31ecc86aabbfd67f38f326ec5a0627d37 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 48eca0284006ed3344469441fbe23bd8709a3809..0eea813183e52d1f489e502f4205b00fcef48833 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 *
* *
diff --git a/TurTLE/test/test_particle_deleter.py b/TurTLE/test/test_particle_deleter.py
new file mode 100644
index 0000000000000000000000000000000000000000..a7bbe788433c201c01e16bf6f27f8ad9cf706b92
--- /dev/null
+++ b/TurTLE/test/test_particle_deleter.py
@@ -0,0 +1,379 @@
+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(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)
+ 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(100)
+ 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:
+ 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_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')
+ 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)
+ 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():
+ #####################################################################
+ ## 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])
+ if index2==1:
+ index_to_delete = np.array([1, 3])
+ after_py = np.delete(before_py, index_to_delete)
+ return before_py, after_py
+ #####################################################################
+
+ #####################################################################
+ ## perform sanity check
+ f = h5py.File(bla.simname + '.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')
+ 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
+ #####################################################################
+ return None
+
+if __name__ == '__main__':
+ main()
diff --git a/cpp/base.hpp b/cpp/base.hpp
index be05de0490194d34148f7afc888281a918169d8f..55562ebcefd600fbeb6087f60552bede030ec64f 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 11c1cbf9e32377f6b59caa90f1350c1a73606d84..5fbcdf81bb3e3e82def8e897ed80a8c803584282 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.cpp b/cpp/field.cpp
index 3efc8ee29007f7a2305a9d53d9d09d3d5280826f..6e950ad450553dbfccfa3ba33f0ae0bab168c58b 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,
@@ -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();
@@ -1812,8 +1813,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]);
@@ -1836,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)
@@ -3085,7 +3085,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;
diff --git a/cpp/field.hpp b/cpp/field.hpp
index a42d480111bd62339bf7d547d3c049685291de23..f6c69a7d1e10d6751e0be0055518b1fa4c2d2e6c 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.
*
@@ -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 0742a2cb0408ea5aee2022aca02d62104fdfea13..c615f9e266fa7de5816c0cd98b5a830dc81d6e48 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,
@@ -72,7 +72,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 edc094a953e16b2817d1d68bb4248c8b4cbc7d13..e7ba578cbeb5f7e990e2c21f95e0cb9cc1054985 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(
@@ -52,7 +53,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 +86,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/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index a9c604ee46352acf85c505d05cfa7c6a72f0279c..6f7cdfd16877f3d207dd8f1a30c5058bf2a03dfb 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 7c1c058d27fda089418c343145478e530f264e3c..25a4c3956f98143bba4be477490499fa67533921 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 9abb40fa2d6beeb96887af79641165bbebb5a2e8..09ca7f286be9d1d17ba3bce5125187bf88fb8024 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.hpp b/cpp/full_code/NSVE.hpp
index df8ffc5daf2aea10dd8891334ff9fd6c1aa62182..2f29f91fef5c990f1676dc547e30458c0ff97a0b 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_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp
index 83a1a7eff29b9b02f017df3476546e5d4d37715f..8da3bcd0ac5dd8091d05b7c938856858dd8ffef8 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 ecc2d135064f245a36d905714e05770e3574d074..fb8497c99074e52c246ac83ba9a071e53d3023c6 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/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index 3063943d4baf4111b06e98ea1acdf5efce50f240..4ee109b88463436a9448a421dad31b5833af76b0 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/NSVE_field_stats.hpp b/cpp/full_code/NSVE_field_stats.hpp
index de2078eedf23a576a58b7faca1f3b34bb314bf16..a04215389a7773f8ecbdf052d4a2a6fd4119b1b0 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 d912aec40fe32e5423300495251f4cd0f55dc40c..d739745cff1df3fe98ad555b0c1c919834cd64da 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/NSVEcomplex_particles.cpp b/cpp/full_code/NSVEcomplex_particles.cpp
index f42518f77c90eda1efb73a6a769c06bbd3123f75..0c71cf723ca18e982a9b684a76b9f9c9f8834ae9 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 a63995f986d9793307b3c57230866ddde2b25e84..3f6adaac45e79594e1e15a1e1fde737e4a1ea27f 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 e9bdfaf760565ba477f8634d34b98b5719143ff2..f2c53a4eadb5770ffecb4bb440fbc5e899b662a6 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();
@@ -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 bbb8ca9f508f3a4e4b1e0b6279b68cb7ffa847a6..dfa5598a911c3a6f7071312a8c22a43e9d0095e7 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/bandpass_stats.cpp b/cpp/full_code/bandpass_stats.cpp
index b176f892abdc5338dee3cd9a6127148c968bed44..7985a26af9be291915a0fcc088d4ccc084a6280f 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>
@@ -47,6 +48,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 +70,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/code_base.cpp b/cpp/full_code/code_base.cpp
index 1015aed91151bbcc94b5516a30a60ff6cf96f516..b62f50150fe815fc554c5bf63eab8568fa0218e6 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/code_base.hpp b/cpp/full_code/code_base.hpp
index d42a7f99f10f5f71a6131e8d8bb72192315623e5..017785d4a3e52cd92c06c2fa356a6c1c76763eaa 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/dealias_test.cpp b/cpp/full_code/dealias_test.cpp
index e3afa2ac17fe66b79b4354310203e700c9601cf7..56980b54a9bbb8c7e35574e7c26908cd948a831f 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/dealias_test.hpp b/cpp/full_code/dealias_test.hpp
index 1757892819d68ce89b7c1296b8bbeba4fa9ae1ba..ba765c45086760910cd8bbc7270b55189e58b48e 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/direct_numerical_simulation.hpp b/cpp/full_code/direct_numerical_simulation.hpp
index 668fa53f50d6a86dd1c281001af7ca1e64feb4e9..be5d498866372c2a7843fbac0768add479dcc271 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
@@ -50,7 +51,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/field_single_to_double.cpp b/cpp/full_code/field_single_to_double.cpp
index 37eb482f511abd495f3999f029acbb0ad90bfcb3..9f2082e63613e7c2d64943c01d8c782f57e4a1cf 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>
@@ -58,6 +59,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 +70,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/field_single_to_double.hpp b/cpp/full_code/field_single_to_double.hpp
index 04f56357c34b5365e1aa60cbb21ad5c7c60a30b9..de98e63ed5fda3fa9e6ffadfbb076c38436d414f 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/field_test.cpp b/cpp/full_code/field_test.cpp
index 9b4c44605584efef0033c4930a7f48f831eb9603..0aa2b940fe69a5919d7ce6a7c38bc2d1ffb07f93 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 8d0ae00fd0ec0fa5c3fd2085163953d6d3303f90..e70cc61c7ce8d40b606262a451e99f3829958832 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 9a672f1c9cdd08730a4936e719e4a5e115d75f8e..6970e7af4a10e2f97d0a9fca542a0ee65989d800 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)
@@ -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 c3cf6c9b1d7a802f831c4faac939ed07d5552f3c..85864ee4d8f4fda468539a94cd43512a6aea5573 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>
@@ -65,6 +66,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 +79,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 1509854d8d329e98fe2028257315941e7ce4cafc..a29cc660c72fd6e07cfcd41ad9f03482be7f1c1b 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>
@@ -53,6 +54,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 +64,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/joint_acc_vel_stats.cpp b/cpp/full_code/joint_acc_vel_stats.cpp
index f0f500dd3a561486bca0643570040aabc5c3c53f..8fe55d2025e1cf7c0c1b91144db2aab56ba6b4bc 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>
@@ -62,6 +63,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 +79,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
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 32502b47d17d75dc0082036db37667474b39fa43..bb63f5def0bef1b470801aff19bbd2d916e16235 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 0a9d4a3f833e9cd85bc5c896181b53b43a09dd1f..020cd0cdc7b2c5b8e35114bc7d78106e9a8c1fb2 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/native_binary_to_hdf5.cpp b/cpp/full_code/native_binary_to_hdf5.cpp
index a1a3277c0a0235ca97751accf55af2211eea25c2..ac0973d61144e0e7bdd5bb08533ece0779d7d757 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/native_binary_to_hdf5.hpp b/cpp/full_code/native_binary_to_hdf5.hpp
index d8ded73b4f24d6f73ae03c3b148ff740d72d8a88..3baaf609bfa4bba093dd097dce3be3d3aa590a4e 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 8e805300cd38d2283e8cb5eb8fb1f276d35f1fa7..e98f76cd568239a88f15da18c4fa72bc255fae9b 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 c232a6b3d2395515f94ce3ad43178f0efa825da6..da22c91e9f7ca187cf7b3038994c0852e255924b 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.cpp b/cpp/full_code/pressure_stats.cpp
index fc060b4dc152ec5a0a578879e82329a0a214b3a6..abeccf5cacca445032bd1e1f8898522f280c745f 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>
@@ -83,6 +84,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 +97,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/pressure_stats.hpp b/cpp/full_code/pressure_stats.hpp
index d8a54494660198adab52a0d12277b2fda407a317..f20ba08ee767fcd44bf6017670f3a3a41d7ba96e 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/resize.cpp b/cpp/full_code/resize.cpp
index f074a46c1d734b1873b9b09a1eda7d10aaed18cc..12f125e026e9faf864b5d510d59371390d5874e9 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>
@@ -43,6 +44,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 +64,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/static_field.cpp b/cpp/full_code/static_field.cpp
index 25760cf322debb4bb2438bb5242ef8fb5c01b96e..1fa52b646deb48f53b2066440221e10b01c3ab55 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 6393ae21151f2b6475972fdab991b808804a091a..90b563d1cd777768a63591e78ef121c9a594fc27 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/symmetrize_test.cpp b/cpp/full_code/symmetrize_test.cpp
index da2bef880c13dd76221477e2088dafdff6453892..904b57ae2e32004c97f0ca1046e450822f853809 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/symmetrize_test.hpp b/cpp/full_code/symmetrize_test.hpp
index 52fbac28c57efb85518b69384f3e0319c60dfa26..69443d1a8dbeee28d035c5c2364f2d79dce8f8af 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/test_interpolation.cpp b/cpp/full_code/test_interpolation.cpp
index 7325b894bac280bcf384da46a1734a52845d9727..9b9f2f50b53d05205eef6ece4287f2e8f2b136a2 100644
--- a/cpp/full_code/test_interpolation.cpp
+++ b/cpp/full_code/test_interpolation.cpp
@@ -101,7 +101,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 59ff49b81ea146b465b4e18c907e1f36ffbfa2bd..7c9e0fff1f1d286865fd1f29b62005063f141e75 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_interpolation_methods.cpp b/cpp/full_code/test_interpolation_methods.cpp
index 5d4a63760b8ea2b1b730eba3db137610e79eebd1..17cdfaf0ad75a02f78859bdedf61edf8536d9929 100644
--- a/cpp/full_code/test_interpolation_methods.cpp
+++ b/cpp/full_code/test_interpolation_methods.cpp
@@ -123,7 +123,7 @@ int test_interpolation_methods<rnumber>::do_work()
// initialize writer
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
this->comm,
pset.getTotalNumberOfParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/test_interpolation_methods.hpp b/cpp/full_code/test_interpolation_methods.hpp
index d0c43b9fbabcf3be9f52548b7d0f1eabc3fb12a7..62f661729cfbec51a25e6bb73819c8daac6dec66 100644
--- a/cpp/full_code/test_interpolation_methods.hpp
+++ b/cpp/full_code/test_interpolation_methods.hpp
@@ -47,7 +47,7 @@ class test_interpolation_methods: public test
int nzparticles;
int random_seed;
- 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, ONE> *phi;
@@ -59,10 +59,10 @@ class test_interpolation_methods: public test
test(
COMMUNICATOR,
simulation_name),
+ random_seed(0),
particles_sample_writer_mpi(nullptr),
phi(nullptr),
- kk(nullptr),
- random_seed(0){}
+ kk(nullptr){}
~test_interpolation_methods(){}
int initialize(void);
diff --git a/cpp/full_code/test_particle_integration.cpp b/cpp/full_code/test_particle_integration.cpp
index 92d5eff4582b20f6aa9949606b8924d705767002..a1d3cbaf957ffb6dab9070746ac70d3d22f57695 100644
--- a/cpp/full_code/test_particle_integration.cpp
+++ b/cpp/full_code/test_particle_integration.cpp
@@ -196,7 +196,7 @@ int test_particle_integration<rnumber>::do_work()
// initialize writer
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
- long long int, double, 3>(
+ long long int, double, double, 3>(
this->comm,
pset.getTotalNumberOfParticles(),
(this->simname + "_particles.h5"),
diff --git a/cpp/full_code/test_particle_integration.hpp b/cpp/full_code/test_particle_integration.hpp
index c303e080d21b04600720bd43207214af24a14a4a..88a32dc19216887d619832fddcf6c5dc6ec18118 100644
--- a/cpp/full_code/test_particle_integration.hpp
+++ b/cpp/full_code/test_particle_integration.hpp
@@ -53,7 +53,7 @@ class test_particle_integration: public test
field_tinterpolator<rnumber, FFTW, THREE, NONE> *velocity_front;
- 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_back;
@@ -65,11 +65,11 @@ class test_particle_integration: public test
test(
COMMUNICATOR,
simulation_name),
- particles_sample_writer_mpi(nullptr),
+ random_seed(0),
velocity_front(nullptr),
+ particles_sample_writer_mpi(nullptr),
velocity_back(nullptr),
- kk(nullptr),
- random_seed(0){}
+ kk(nullptr){}
~test_particle_integration(){}
int initialize(void);
diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
index d457a090b06723c08ffc523fd483f47471de629d..ec8893577686cab0410ec49b396133535868b586 100644
--- a/cpp/full_code/test_tracer_set.cpp
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -109,8 +109,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/full_code/write_rpressure.cpp b/cpp/full_code/write_rpressure.cpp
index 5835ee33625b6ac0240bb3c46e5aa9932b3a514a..a9e5e4ffa6972030c597befdb7e1020464af3f73 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)
@@ -20,6 +21,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
diff --git a/cpp/full_code/write_rpressure.hpp b/cpp/full_code/write_rpressure.hpp
index 64fbcfc75a5feb17871fba515e089fd20f6a22a3..2a14dcb6d58188cb5ad7b6d9333cff76ef963282 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);
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 920db174b8e665ec8e9c1101fd98b6d15ff8c9c2..0d0854f92e023db0ca8d05301bf53d38b3e51dc3 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/kspace.hpp b/cpp/kspace.hpp
index 656b8ca290b32a9cf87bbdb0305ac7e8a2b7fe43..6457c8bd0a65d229d154984c05341f69a7728f7b 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};
@@ -78,7 +78,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/particles/abstract_particle_rhs.hpp b/cpp/particles/abstract_particle_rhs.hpp
index d1a47bdcf164f24d1617e9f3f9f7b972fe06c2e5..947e152430cb14bc439021d83e48860411c2c175 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 ab078a423834b075e05d864b25d4b07b994a6a49..20884283491a8f316277ef78be20026080a34e9a 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;
@@ -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 2674946190a9de678c3c839c79f8704715426cbf..af690dcbd72ca06cf756d1d639a6efaebfd6b677 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;
@@ -89,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),
@@ -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]);
}
}
}
@@ -215,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;
@@ -250,23 +262,25 @@ 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){
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 +289,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,18 +304,29 @@ 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]);
}
}
{
+ // 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);
assert(dst_idx < particles_chunk_current_size);
@@ -322,7 +347,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 abcb8e684c5f269f10f47f337e36260486fba931..4edf20c7ca2d34fb023c6a0ea1e370c29e288df0 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 085956f4d377ef9ba24f42995a844d4f10da1a6f..1227f0667107429c0516aac78b67e904a6b9ee18 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,11 +57,12 @@ 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;
- 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;
@@ -183,7 +185,48 @@ class abstract_particle_set
this->getParticleState(),
additional_data.data(),
int(additional_data.size()),
- this->getParticleIndex());
+ this->getParticleIndices(),
+ this->getParticleLabels());
+ return EXIT_SUCCESS;
+ }
+
+ int writeParticleLabels(
+ const std::string file_name,
+ const std::string species_name,
+ const std::string field_name,
+ const int iteration)
+ {
+ 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(),
+ 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->getParticleLabels(),
+ this->getParticleLabels() + this->getLocalNumberOfParticles(),
+ pdata.get());
+ particle_sample_writer->template save_dataset<1>(
+ species_name,
+ field_name,
+ xx.get(),
+ &pdata,
+ this->getParticleIndices(),
+ this->getLocalNumberOfParticles(),
+ iteration);
+ // deallocate sample writer
+ delete particle_sample_writer;
+ // deallocate temporary array
+ delete[] pdata.release();
+ // deallocate position array
+ delete[] xx.release();
return EXIT_SUCCESS;
}
@@ -192,7 +235,10 @@ class abstract_particle_set
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)
@@ -212,7 +258,7 @@ class abstract_particle_set
field_name,
xx.get(),
&pdata,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
// deallocate temporary array
@@ -257,7 +303,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)
@@ -279,7 +328,7 @@ class abstract_particle_set
field_name,
xx.get(),
&yy,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
break;
@@ -289,7 +338,7 @@ class abstract_particle_set
field_name,
xx.get(),
&yy,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
break;
@@ -299,7 +348,7 @@ class abstract_particle_set
field_name,
xx.get(),
&yy,
- this->getParticleIndex(),
+ this->getParticleIndices(),
this->getLocalNumberOfParticles(),
iteration);
break;
@@ -316,7 +365,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 +384,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/field_tinterpolator.hpp b/cpp/particles/interpolation/field_tinterpolator.hpp
index 60175a288cf55a246e82a7a8baac71431b02eaaa..2dc27b56a69d5462b131347b23c8fce1903864af 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/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 206e9da519432a1226c441517f97d26c5f03c123..071614ad52369d596565a50fa2e8ac34d5c22e99 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"
@@ -73,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;
@@ -159,15 +161,77 @@ 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]);
}
- ~particle_set(){}
+ 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<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");
+ }
+ MPI_Barrier(this->mpi_comm);
+ }
+ return EXIT_SUCCESS;
+ }
+
+ ~particle_set() noexcept(false){}
particle_rnumber* getParticleState() const
{
@@ -192,9 +256,14 @@ class particle_set: public abstract_particle_set
return EXIT_SUCCESS;
}
- partsize_t* getParticleIndex() const
+ partsize_t* getParticleIndices() const
+ {
+ return this->local_indices.get();
+ }
+
+ partsize_t* getParticleLabels() const
{
- return this->local_index.get();
+ return this->local_labels.get();
}
partsize_t getLocalNumberOfParticles() const
@@ -314,7 +383,8 @@ class particle_set: public abstract_particle_set
&this->local_state,
additional_data.data(),
int(additional_data.size()),
- &this->local_index);
+ &this->local_indices,
+ &this->local_labels);
return EXIT_SUCCESS;
}
@@ -323,10 +393,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,
@@ -339,7 +411,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());
@@ -356,25 +429,26 @@ 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();
- 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++)
{
- 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);
@@ -388,13 +462,221 @@ 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;
+
+ 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;
+ }
+
+ 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;
+
+ 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;
+ }
+
+ 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.
+ *
+ * \warning `indices_of_interest` must be in increasing order.
+ *
+ */
+ 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());
+ // 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()];
+ 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] = 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>(
@@ -409,11 +691,21 @@ 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;
+
+ // 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;
}
@@ -431,6 +723,19 @@ 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<partsize_t(inGlobalIndicesToDelete.size()); ii++)
+ {
+ if (this->getParticleIndices()[i] == inGlobalIndicesToDelete[ii])
+ local_indices.push_back(inGlobalIndicesToDelete[ii]);
+ }
+ return local_indices;
+ }
};
#endif//PARTICLE_SET_HPP
diff --git a/cpp/particles/lock_free_bool_array.hpp b/cpp/particles/lock_free_bool_array.hpp
index bfcd48880cfc6308f63afe235f0528b1a927f68b..9c0ab5de0c13ab49e88b3d05403631e49656599b 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/p2p/p2p_distr_mpi.hpp b/cpp/particles/p2p/p2p_distr_mpi.hpp
index ecf6aa4ddf7c3bd56a745f4ebec01009325cf725..3051bcee05b5d6078abdbac571cc4500ff6af475 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{
@@ -197,7 +199,7 @@ public:
nb_cell_levels[IDXC_Z] = nb_cells_factor;
}
- virtual ~p2p_distr_mpi(){}
+ virtual ~p2p_distr_mpi() noexcept(false){}
////////////////////////////////////////////////////////////////////////////
@@ -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],
@@ -434,13 +446,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
{
@@ -549,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});
@@ -675,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]);
@@ -702,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()));
+ }
}
}
@@ -711,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);
@@ -721,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
@@ -818,6 +873,7 @@ public:
&mpiRequests.back()));
delete[] descriptor.toCompute.release();
delete[] descriptor.indexes.release();
+ delete[] descriptor.labels.release();
}
}
//////////////////////////////////////////////////////////////////////
diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp
index f53b24bdf731f70dd7437966baecac1571583f53..e22290e0a815f9e072b5fdc4a524aa31ef089d9c 100644
--- a/cpp/particles/particle_solver.cpp
+++ b/cpp/particles/particle_solver.cpp
@@ -315,7 +315,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/particle_solver.hpp b/cpp/particles/particle_solver.hpp
index 004581c8961364714fdbfd775ad683be43748117..54bbc3c846022c00f220d9f427e121241fd10223 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/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index 7a4d646c2e663aa31a972a8c488e9c6aca6f6e18..0a468253dc4d38d362d2dec189445de55df3ec7e 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,
@@ -160,7 +163,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){}
////////////////////////////////////////////////////////////////////////////
@@ -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);
+
{// this fails when particles move more than one cell during a single time-step
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){
diff --git a/cpp/particles/particles_input_grid.hpp b/cpp/particles/particles_input_grid.hpp
index e132095ee20c527039d23824a5243816957b32df..784c08e0623b602dc6e81a7b4bd13ed3c8a8dfe0 100644
--- a/cpp/particles/particles_input_grid.hpp
+++ b/cpp/particles/particles_input_grid.hpp
@@ -45,8 +45,10 @@ class particles_input_grid: public abstract_particles_input<partsize_t, particle
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_grid() noexcept(false){}
particles_input_grid(
const MPI_Comm in_mpi_comm,
const partsize_t NXPARTICLES,
@@ -179,6 +181,12 @@ class particles_input_grid: public abstract_particles_input<partsize_t, particle
}
/*************************************************/
+
+ // 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()
@@ -206,6 +214,12 @@ class particles_input_grid: public abstract_particles_input<partsize_t, particle
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[]>>());
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 680daff93c55ad1e6fe6951dfc3705b1229d17de..b1ea662ed14eaf4e725c532429eb631ade4ef8b7 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>
@@ -54,6 +55,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:
@@ -166,7 +168,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;
@@ -318,9 +320,15 @@ 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(){
+ ~particles_input_hdf5() noexcept(false){
}
partsize_t getTotalNbParticles() final{
@@ -350,6 +358,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 1179b080f4079515b42d16e64a0de9daa7cb8af2..0b40ae57f70729b82984c4c5b4591568e4078485 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[]>>());
diff --git a/cpp/particles/particles_output_hdf5.hpp b/cpp/particles/particles_output_hdf5.hpp
index 236504e15d00642723eb9cf7b46db042028b2a9c..0ef809e63114948509c96f5cbccc5cdd615b3599 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 b05579e2e3cf013689ef0b2e30b763d285952e86..ee9bd9cee42d3973e1e9e1bb93eb67a6329836d3 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 11a7dcc13ad2b292a9f8215b669bdcf4c97addb9..e3cc2e021b7062471f384d1eac5c288afc882939 100644
--- a/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -27,18 +27,21 @@
#define PARTICLES_OUTPUT_SAMPLING_HDF5_HPP
#include "abstract_particles_output.hpp"
+#include "hdf5_tools.hpp"
-#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,12 +200,11 @@ 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 ?
- H5T_NATIVE_DOUBLE :
- H5T_NATIVE_FLOAT);
+ static_assert(std::is_same<output_type, double>::value ||
+ 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);
diff --git a/cpp/particles/particles_sampling.hpp b/cpp/particles/particles_sampling.hpp
index 1606b0228b334e15355644610e02410ec4253653..8e83a5a53690b6cb3abaddaafe2245b1f448d164 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 b47763028aa25d040c51ffcde11a3dd3cfe7a87d..9589e4dcc2738c554f7ff364235b7bc92d8a155a 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) {
diff --git a/cpp/particles/rhs/deformation_tensor_rhs.hpp b/cpp/particles/rhs/deformation_tensor_rhs.hpp
index cbf43cbce49e039cc8ebb794dd541caed6889853..0835f9869645d14eacedea5990562dffbd0203e4 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 03bfab35d4803ff039f910a4d92305f455e26711..680b46fdbc330d3bd892a1343094dd066c1f83a5 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 9ba41c3092fab1f5de5605e050e9a7facb670511..5aa15213d02702eb96c198fab0e333584efb1598 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
{
diff --git a/cpp/scope_timer.hpp b/cpp/scope_timer.hpp
index 91365e01695b9a4b577d8863302b9fbed1002a12..211cccd49cabd63b3f874150406b285417a8f72d 100644
--- a/cpp/scope_timer.hpp
+++ b/cpp/scope_timer.hpp
@@ -244,7 +244,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);
@@ -416,7 +416,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();
diff --git a/cpp/spectrum_function.hpp b/cpp/spectrum_function.hpp
index fcc729ff4763f2a417b8da20a4e5dc0aab31110d..7be4474c0a6151e4eed81248140ad95fc2cdec51 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>
@@ -20,7 +23,7 @@ class spectrum_function
{
assert(this->values.size() == this->kk->nshells);
}
- ~spectrum_function(){}
+ ~spectrum_function() noexcept(false){}
double operator()(double kvalue)
{
@@ -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 8bf2ea8ef61d8f022b2e70280bac8916e514dee9..32c6547b9781d23bf180d2b28a13d5cfdcd4e69a 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;
@@ -87,7 +88,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);
diff --git a/setup.py.in b/setup.py.in
index 28eb48704653c4ec6ba90792f86a2c570bc4562a..1dd8a24eac4b11b30300391ee88363a883b4dce8 100644
--- a/setup.py.in
+++ b/setup.py.in
@@ -48,7 +48,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,
########################################################################