From c27efd2b833fe19087b1f69d08b304eb28fc3287 Mon Sep 17 00:00:00 2001
From: Berenger Bramas <bbramas@mpcdf.mpg.de>
Date: Fri, 19 May 2017 11:47:55 +0200
Subject: [PATCH] Add a sampling method -- update to be generic against the
field -- need to update NSVEparticles::do_stats
---
bfps/cpp/full_code/NSVEparticles.cpp | 13 +-
.../particles/abstract_particles_output.hpp | 8 +-
.../particles/abstract_particles_system.hpp | 21 +++
...cles_distr.hpp => particles_distr_mpi.hpp} | 70 ++++------
.../particles/particles_field_computer.hpp | 22 +---
.../particles_output_sampling_hdf5.hpp | 123 ++++++++++++++++++
bfps/cpp/particles/particles_sampling.hpp | 40 ++++++
bfps/cpp/particles/particles_system.hpp | 70 ++++++++--
.../particles/particles_system_builder.hpp | 27 ++--
setup.py | 6 +-
10 files changed, 308 insertions(+), 92 deletions(-)
rename bfps/cpp/particles/{abstract_particles_distr.hpp => particles_distr_mpi.hpp} (93%)
create mode 100644 bfps/cpp/particles/particles_output_sampling_hdf5.hpp
create mode 100644 bfps/cpp/particles/particles_sampling.hpp
diff --git a/bfps/cpp/full_code/NSVEparticles.cpp b/bfps/cpp/full_code/NSVEparticles.cpp
index 9c0603b1..35cbf47d 100644
--- a/bfps/cpp/full_code/NSVEparticles.cpp
+++ b/bfps/cpp/full_code/NSVEparticles.cpp
@@ -2,7 +2,7 @@
#include <cmath>
#include "NSVEparticles.hpp"
#include "scope_timer.hpp"
-
+#include "particles/particles_sampling.hpp"
template <typename rnumber>
int NSVEparticles<rnumber>::initialize(void)
@@ -69,9 +69,14 @@ int NSVEparticles<rnumber>::do_stats()
// fluid stats go here
this->NSVE<rnumber>::do_stats();
// particle sampling should go here
- //if (this->iteration % this->niter_part == 0)
- //{
- //}
+// if (this->iteration % this->niter_part == 0)
+// {
+// sample_from_particles_system(*this->fs->cvelocity,// field to save TODO
+// this->ps,
+// 0, // hdf5 datagroup TODO
+// "TODO", // dataset basename TODO
+// this->particles_output_writer_mpi->getTotalNbParticles());
+// }
return EXIT_SUCCESS;
}
diff --git a/bfps/cpp/particles/abstract_particles_output.hpp b/bfps/cpp/particles/abstract_particles_output.hpp
index dfc8ec37..a6eccaea 100644
--- a/bfps/cpp/particles/abstract_particles_output.hpp
+++ b/bfps/cpp/particles/abstract_particles_output.hpp
@@ -45,10 +45,6 @@ protected:
return mpi_com_writer;
}
- partsize_t getTotalNbParticles() const {
- return total_nb_particles;
- }
-
int getNbRhs() const {
return nb_rhs;
}
@@ -126,6 +122,10 @@ public:
if(current_is_involved){
AssertMpi( MPI_Comm_free(&mpi_com_writer) );
}
+ }
+
+ partsize_t getTotalNbParticles() const {
+ return total_nb_particles;
}
void releaseMemory(){
diff --git a/bfps/cpp/particles/abstract_particles_system.hpp b/bfps/cpp/particles/abstract_particles_system.hpp
index 33a94513..6961eabc 100644
--- a/bfps/cpp/particles/abstract_particles_system.hpp
+++ b/bfps/cpp/particles/abstract_particles_system.hpp
@@ -3,6 +3,12 @@
#include <memory>
+//- Not generic to enable sampling begin
+#include "field.hpp"
+#include "kspace.hpp"
+//- Not generic to enable sampling end
+
+
template <class partsize_t, class real_number>
class abstract_particles_system {
public:
@@ -27,6 +33,21 @@ public:
virtual partsize_t getLocalNbParticles() const = 0;
virtual int getNbRhs() const = 0;
+
+ //- Not generic to enable sampling begin
+ virtual void sample_compute_field(const field<float, FFTW, ONE>& sample_field,
+ real_number sample_rhs[]) = 0;
+ virtual void sample_compute_field(const field<float, FFTW, THREE>& sample_field,
+ real_number sample_rhs[]) = 0;
+ virtual void sample_compute_field(const field<float, FFTW, THREExTHREE>& sample_field,
+ real_number sample_rhs[]) = 0;
+ virtual void sample_compute_field(const field<double, FFTW, ONE>& sample_field,
+ real_number sample_rhs[]) = 0;
+ virtual void sample_compute_field(const field<double, FFTW, THREE>& sample_field,
+ real_number sample_rhs[]) = 0;
+ virtual void sample_compute_field(const field<double, FFTW, THREExTHREE>& sample_field,
+ real_number sample_rhs[]) = 0;
+ //- Not generic to enable sampling end
};
#endif
diff --git a/bfps/cpp/particles/abstract_particles_distr.hpp b/bfps/cpp/particles/particles_distr_mpi.hpp
similarity index 93%
rename from bfps/cpp/particles/abstract_particles_distr.hpp
rename to bfps/cpp/particles/particles_distr_mpi.hpp
index b012f824..48559518 100644
--- a/bfps/cpp/particles/abstract_particles_distr.hpp
+++ b/bfps/cpp/particles/particles_distr_mpi.hpp
@@ -1,5 +1,5 @@
-#ifndef ABSTRACT_PARTICLES_DISTR_HPP
-#define ABSTRACT_PARTICLES_DISTR_HPP
+#ifndef PARTICLES_DISTR_MPI_HPP
+#define PARTICLES_DISTR_MPI_HPP
#include <mpi.h>
@@ -14,8 +14,8 @@
#include "particles_utils.hpp"
-template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs, int size_particle_index>
-class abstract_particles_distr {
+template <class partsize_t, class real_number>
+class particles_distr_mpi {
protected:
static const int MaxNbRhs = 100;
@@ -92,7 +92,7 @@ protected:
public:
////////////////////////////////////////////////////////////////////////////
- abstract_particles_distr(MPI_Comm in_current_com,
+ particles_distr_mpi(MPI_Comm in_current_com,
const std::pair<int,int>& in_current_partitions,
const std::array<size_t,3>& in_field_grid_dim)
: current_com(in_current_com),
@@ -130,11 +130,14 @@ public:
assert(int(field_grid_dim[IDX_Z]) == partition_interval_offset_per_proc[nb_processes_involved]);
}
- virtual ~abstract_particles_distr(){}
+ virtual ~particles_distr_mpi(){}
////////////////////////////////////////////////////////////////////////////
- void compute_distr(const partsize_t current_my_nb_particles_per_partition[],
+ template <class computer_class, class field_class, int size_particle_positions, int size_particle_rhs>
+ void compute_distr(computer_class& in_computer,
+ field_class& in_field,
+ const partsize_t current_my_nb_particles_per_partition[],
const real_number particles_positions[],
real_number particles_current_rhs[],
const int interpolation_size){
@@ -377,10 +380,10 @@ public:
assert(descriptor.toCompute != nullptr);
descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
- init_result_array(descriptor.results.get(), NbParticlesToReceive);
+ in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
if(more_than_one_thread == false){
- apply_computation(descriptor.toCompute.get(), descriptor.results.get(), NbParticlesToReceive);
+ in_computer.template apply_computation<field_class, size_particle_rhs>(in_field, descriptor.toCompute.get(), descriptor.results.get(), NbParticlesToReceive);
}
else{
TIMEZONE_OMP_INIT_PRETASK(timeZoneTaskKey)
@@ -392,8 +395,8 @@ public:
#pragma omp task default(shared) firstprivate(ptr_descriptor, idxPart, sizeToDo) priority(10) \
TIMEZONE_OMP_PRAGMA_TASK_KEY(timeZoneTaskKey)
{
- TIMEZONE_OMP_TASK("apply_computation", timeZoneTaskKey);
- apply_computation(&ptr_descriptor->toCompute[idxPart*size_particle_positions],
+ TIMEZONE_OMP_TASK("in_computer.apply_computation", timeZoneTaskKey);
+ in_computer.template apply_computation<field_class, size_particle_rhs>(in_field, &ptr_descriptor->toCompute[idxPart*size_particle_positions],
&ptr_descriptor->results[idxPart*size_particle_rhs], sizeToDo);
}
}
@@ -425,13 +428,13 @@ public:
if(descriptor.isLower){
TIMEZONE("reduce");
assert(descriptor.toRecvAndMerge != nullptr);
- reduce_particles_rhs(&particles_current_rhs[0], descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToSend);
+ in_computer.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[0], descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToSend);
descriptor.toRecvAndMerge.release();
}
else {
TIMEZONE("reduce");
assert(descriptor.toRecvAndMerge != nullptr);
- reduce_particles_rhs(&particles_current_rhs[(current_offset_particles_for_partition[current_partition_size]-descriptor.nbParticlesToSend)*size_particle_rhs],
+ in_computer.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[(current_offset_particles_for_partition[current_partition_size]-descriptor.nbParticlesToSend)*size_particle_rhs],
descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToSend);
descriptor.toRecvAndMerge.release();
}
@@ -452,8 +455,8 @@ public:
// Low priority to help master thread when possible
#pragma omp task default(shared) firstprivate(idxPart, sizeToDo) priority(0) TIMEZONE_OMP_PRAGMA_TASK_KEY(timeZoneTaskKey)
{
- TIMEZONE_OMP_TASK("apply_computation", timeZoneTaskKey);
- apply_computation(&particles_positions[idxPart*size_particle_positions],
+ TIMEZONE_OMP_TASK("in_computer.apply_computation", timeZoneTaskKey);
+ in_computer.template apply_computation<field_class, size_particle_rhs>(in_field, &particles_positions[idxPart*size_particle_positions],
&particles_current_rhs[idxPart*size_particle_rhs],
sizeToDo);
}
@@ -470,13 +473,13 @@ public:
if(descriptor.isLower){
TIMEZONE("reduce_later");
assert(descriptor.toRecvAndMerge != nullptr);
- reduce_particles_rhs(&particles_current_rhs[0], descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToSend);
+ in_computer.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[0], descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToSend);
descriptor.toRecvAndMerge.release();
}
else {
TIMEZONE("reduce_later");
assert(descriptor.toRecvAndMerge != nullptr);
- reduce_particles_rhs(&particles_current_rhs[(current_offset_particles_for_partition[current_partition_size]-descriptor.nbParticlesToSend)*size_particle_rhs],
+ in_computer.template reduce_particles_rhs<size_particle_rhs>(&particles_current_rhs[(current_offset_particles_for_partition[current_partition_size]-descriptor.nbParticlesToSend)*size_particle_rhs],
descriptor.toRecvAndMerge.get(), descriptor.nbParticlesToSend);
descriptor.toRecvAndMerge.release();
}
@@ -489,7 +492,7 @@ public:
TIMEZONE("compute-my_compute");
// Compute my particles
if(myTotalNbParticles){
- apply_computation(particles_positions, particles_current_rhs, myTotalNbParticles);
+ in_computer.template apply_computation<field_class, size_particle_rhs>(in_field, particles_positions, particles_current_rhs, myTotalNbParticles);
}
}
@@ -497,20 +500,12 @@ public:
assert(mpiRequests.size() == 0);
}
- ////////////////////////////////////////////////////////////////////////////
-
- virtual void init_result_array(real_number particles_current_rhs[],
- const partsize_t nb_particles) const = 0;
- virtual void apply_computation(const real_number particles_positions[],
- real_number particles_current_rhs[],
- const partsize_t nb_particles) const = 0;
- virtual void reduce_particles_rhs(real_number particles_current_rhs[],
- const real_number extra_particles_current_rhs[],
- const partsize_t nb_particles) const = 0;
////////////////////////////////////////////////////////////////////////////
- void redistribute(partsize_t current_my_nb_particles_per_partition[],
+ template <class computer_class, int size_particle_positions, int size_particle_rhs, int size_particle_index>
+ void redistribute(computer_class& in_computer,
+ partsize_t current_my_nb_particles_per_partition[],
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,
@@ -533,7 +528,7 @@ public:
// Find particles outside my interval
const partsize_t nbOutLower = particles_utils::partition_extra<partsize_t, size_particle_positions>(&(*inout_positions_particles)[0], current_my_nb_particles_per_partition[0],
[&](const real_number val[]){
- const int partition_level = pbc_field_layer(val[IDX_Z], IDX_Z);
+ const int partition_level = in_computer.pbc_field_layer(val[IDX_Z], IDX_Z);
assert(partition_level == current_partition_interval.first
|| partition_level == (current_partition_interval.first-1+int(field_grid_dim[IDX_Z]))%int(field_grid_dim[IDX_Z])
|| partition_level == (current_partition_interval.first+1)%int(field_grid_dim[IDX_Z]));
@@ -558,7 +553,7 @@ public:
&(*inout_positions_particles)[(current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow)*size_particle_positions],
myTotalNbParticles - (current_offset_particles_for_partition[current_partition_size-1]+offesetOutLow),
[&](const real_number val[]){
- const int partition_level = pbc_field_layer(val[IDX_Z], IDX_Z);
+ const int partition_level = in_computer.pbc_field_layer(val[IDX_Z], IDX_Z);
assert(partition_level == (current_partition_interval.second-1)
|| partition_level == ((current_partition_interval.second-1)-1+int(field_grid_dim[IDX_Z]))%int(field_grid_dim[IDX_Z])
|| partition_level == ((current_partition_interval.second-1)+1)%int(field_grid_dim[IDX_Z]));
@@ -807,7 +802,7 @@ public:
myTotalNbParticles,current_partition_size,
current_my_nb_particles_per_partition, current_offset_particles_for_partition.get(),
[&](const real_number& z_pos){
- const int partition_level = pbc_field_layer(z_pos, IDX_Z);
+ const int partition_level = in_computer.pbc_field_layer(z_pos, IDX_Z);
assert(current_partition_interval.first <= partition_level && partition_level < current_partition_interval.second);
return partition_level - current_partition_interval.first;
},
@@ -829,7 +824,7 @@ public:
assert(current_my_nb_particles_per_partition[idxPartition] ==
current_offset_particles_for_partition[idxPartition+1] - current_offset_particles_for_partition[idxPartition]);
for(partsize_t idx = current_offset_particles_for_partition[idxPartition] ; idx < current_offset_particles_for_partition[idxPartition+1] ; ++idx){
- assert(pbc_field_layer((*inout_positions_particles)[idx*3+IDX_Z], IDX_Z)-current_partition_interval.first == idxPartition);
+ assert(in_computer.pbc_field_layer((*inout_positions_particles)[idx*3+IDX_Z], IDX_Z)-current_partition_interval.first == idxPartition);
}
}
}
@@ -838,15 +833,6 @@ public:
assert(mpiRequests.size() == 0);
}
-
- virtual int pbc_field_layer(const real_number& a_z_pos, const int idx_dim) const = 0;
-
- ////////////////////////////////////////////////////////////////////////////
-
- virtual void move_particles(real_number particles_positions[],
- const partsize_t nb_particles,
- const std::unique_ptr<real_number[]> particles_current_rhs[],
- const int nb_rhs, const real_number dt) const = 0;
};
#endif
diff --git a/bfps/cpp/particles/particles_field_computer.hpp b/bfps/cpp/particles/particles_field_computer.hpp
index 1c7c7e77..f68f2fc0 100644
--- a/bfps/cpp/particles/particles_field_computer.hpp
+++ b/bfps/cpp/particles/particles_field_computer.hpp
@@ -4,15 +4,13 @@
#include <array>
#include <utility>
-#include "abstract_particles_distr.hpp"
#include "scope_timer.hpp"
#include "particles_utils.hpp"
template <class partsize_t,
class real_number,
class interpolator_class,
- int interp_neighbours,
- class positions_updater_class>
+ int interp_neighbours>
class particles_field_computer {
const std::array<int,3> field_grid_dim;
@@ -20,8 +18,6 @@ class particles_field_computer {
const interpolator_class& interpolator;
- const positions_updater_class positions_updater;
-
const std::array<real_number,3> spatial_box_width;
const std::array<real_number,3> spatial_box_offset;
const std::array<real_number,3> box_step_width;
@@ -33,11 +29,10 @@ public:
particles_field_computer(const std::array<size_t,3>& in_field_grid_dim,
const std::pair<int,int>& in_current_partitions,
const interpolator_class& in_interpolator,
- const field_class& in_field,
const std::array<real_number,3>& in_spatial_box_width, const std::array<real_number,3>& in_spatial_box_offset,
const std::array<real_number,3>& in_box_step_width)
: field_grid_dim({{int(in_field_grid_dim[0]),int(in_field_grid_dim[1]),int(in_field_grid_dim[2])}}), current_partition_interval(in_current_partitions),
- interpolator(in_interpolator), field(in_field), positions_updater(),
+ interpolator(in_interpolator),
spatial_box_width(in_spatial_box_width), spatial_box_offset(in_spatial_box_offset), box_step_width(in_box_step_width){
deriv[IDX_X] = 0;
deriv[IDX_Y] = 0;
@@ -174,19 +169,6 @@ public:
}
}
- ////////////////////////////////////////////////////////////////////////
- /// Update position
- ////////////////////////////////////////////////////////////////////////
-
- void move_particles(real_number particles_positions[],
- const partsize_t nb_particles,
- const std::unique_ptr<real_number[]> particles_current_rhs[],
- const int nb_rhs, const real_number dt) const {
- TIMEZONE("particles_field_computer::move_particles");
- positions_updater.move_particles(particles_positions, nb_particles,
- particles_current_rhs, nb_rhs, dt);
- }
-
////////////////////////////////////////////////////////////////////////
/// Re-distribution related
////////////////////////////////////////////////////////////////////////
diff --git a/bfps/cpp/particles/particles_output_sampling_hdf5.hpp b/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
new file mode 100644
index 00000000..7e9a54dd
--- /dev/null
+++ b/bfps/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -0,0 +1,123 @@
+#ifndef PARTICLES_OUTPUT_SAMPLING_HDF5_HPP
+#define PARTICLES_OUTPUT_SAMPLING_HDF5_HPP
+
+#include "abstract_particles_output.hpp"
+
+#include <hdf5.h>
+
+template <class partsize_t,
+ class real_number,
+ int size_particle_positions,
+ int size_particle_rhs>
+class particles_output_sampling_hdf5 : public abstract_particles_output<partsize_t,
+ real_number,
+ size_particle_positions,
+ size_particle_rhs>{
+ using Parent = abstract_particles_output<partsize_t,
+ real_number,
+ size_particle_positions,
+ size_particle_rhs>;
+
+ hid_t parent_group;
+ const std::string dataset_name;
+ const bool use_collective_io;
+
+public:
+ particles_output_sampling_hdf5(MPI_Comm in_mpi_com,
+ const partsize_t inTotalNbParticles,
+ const hid_t in_parent_group,
+ const std::string& in_dataset_name,
+ const bool in_use_collective_io = false)
+ : Parent(in_mpi_com, inTotalNbParticles, 1),
+ parent_group(in_parent_group),
+ dataset_name(in_dataset_name),
+ use_collective_io(in_use_collective_io){}
+
+ void write(
+ const int idx_time_step,
+ const real_number* /*particles_positions*/,
+ const std::unique_ptr<real_number[]>* particles_rhs,
+ const partsize_t nb_particles,
+ const partsize_t particles_idx_offset) final{
+ assert(Parent::isInvolved());
+
+ TIMEZONE("particles_output_hdf5::write");
+
+ assert(particles_idx_offset < Parent::getTotalNbParticles() || (particles_idx_offset == Parent::getTotalNbParticles() && 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);
+
+ hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
+ assert(plist_id >= 0);
+ {
+ int rethdf = H5Pset_dxpl_mpio(plist_id, use_collective_io ? H5FD_MPIO_COLLECTIVE : H5FD_MPIO_INDEPENDENT);
+ assert(rethdf >= 0);
+ }
+ {
+ assert(size_particle_rhs >= 0);
+ const hsize_t datacount[3] = {hsize_t(Parent::getNbRhs()),
+ hsize_t(Parent::total_nb_particles),
+ hsize_t(size_particle_rhs)};
+ hid_t dataspace = H5Screate_simple(3, datacount, NULL);
+ assert(dataspace >= 0);
+
+ hid_t dataset_id = H5Dcreate( parent_group,
+ dataset_name.c_str(),
+ type_id,
+ dataspace,
+ H5P_DEFAULT,
+ H5P_DEFAULT,
+ H5P_DEFAULT);
+ assert(dataset_id >= 0);
+
+ assert(particles_idx_offset >= 0);
+ const hsize_t count[3] = {
+ 1,
+ hsize_t(nb_particles),
+ hsize_t(size_particle_rhs)};
+ const hsize_t offset[3] = {
+ 0,
+ hsize_t(particles_idx_offset),
+ 0};
+ hid_t memspace = H5Screate_simple(3, count, NULL);
+ assert(memspace >= 0);
+
+ hid_t filespace = H5Dget_space(dataset_id);
+ assert(filespace >= 0);
+ int rethdf = H5Sselect_hyperslab(
+ filespace,
+ H5S_SELECT_SET,
+ offset,
+ NULL,
+ count,
+ NULL);
+ assert(rethdf >= 0);
+
+ herr_t status = H5Dwrite(
+ dataset_id,
+ type_id,
+ memspace,
+ filespace,
+ plist_id,
+ particles_rhs[0].get());
+ assert(status >= 0);
+ rethdf = H5Sclose(filespace);
+ assert(rethdf >= 0);
+ rethdf = H5Sclose(memspace);
+ assert(rethdf >= 0);
+ rethdf = H5Dclose(dataset_id);
+ assert(rethdf >= 0);
+ }
+
+ {
+ int rethdf = H5Pclose(plist_id);
+ assert(rethdf >= 0);
+ }
+ }
+};
+
+#endif
diff --git a/bfps/cpp/particles/particles_sampling.hpp b/bfps/cpp/particles/particles_sampling.hpp
new file mode 100644
index 00000000..cc987f69
--- /dev/null
+++ b/bfps/cpp/particles/particles_sampling.hpp
@@ -0,0 +1,40 @@
+#ifndef PARTICLES_SAMPLING_HPP
+#define PARTICLES_SAMPLING_HPP
+
+#include <memory>
+#include <string>
+
+#include "abstract_particles_system.hpp"
+#include "particles_output_sampling_hdf5.hpp"
+
+#include "field.hpp"
+#include "kspace.hpp"
+
+
+template <class partsize_t, class particles_rnumber, class rnumber, field_backend be, field_components fc>
+void sample_from_particles_system(const field<rnumber, be, fc>& in_field, // a pointer to a field<rnumber, FFTW, fc>
+ std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>& ps, // a pointer to an particles_system<double>
+ hid_t gid, // an hid_t identifying an HDF5 group
+ const std::string& fname,
+ const partsize_t totalNbParticles){
+// const int size_particle_rhs = ncomp(fc);
+// const partsize_t nb_particles = ps->getLocalNbParticles();
+// std::unique_ptr<particles_rnumber[]> sample_rhs(new particles_rnumber[size_particle_rhs*nb_particles]);
+
+// ps->sample_compute_field(in_field, sample_rhs.get());
+
+// const std::string datasetname = fname + std::string("/") + std::to_string(ps->step_idx);
+
+// particles_output_sampling_hdf5<partsize_t, particles_rnumber, 3, size_particle_rhs> outputclass(MPI_COMM_WORLD,
+// totalNbParticles,
+// gid,
+// datasetname);
+// outputclass.save(ps->getParticlesPositions(),
+// &sample_rhs,
+// ps->getParticlesIndexes(),
+// ps->getLocalNbParticles(),
+// in_field->iteration);
+}
+
+#endif
+
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 444f0cc7..1e5ec2dd 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -4,6 +4,7 @@
#include <array>
#include "abstract_particles_system.hpp"
+#include "particles_distr_mpi.hpp"
#include "particles_output_hdf5.hpp"
#include "particles_output_mpiio.hpp"
#include "particles_field_computer.hpp"
@@ -12,7 +13,7 @@
#include "scope_timer.hpp"
template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours,
- class size_particle_rhs>
+ int size_particle_rhs>
class particles_system : public abstract_particles_system<partsize_t, real_number> {
MPI_Comm mpi_com;
@@ -21,7 +22,14 @@ class particles_system : public abstract_particles_system<partsize_t, real_numbe
interpolator_class interpolator;
- particles_field_computer<partsize_t, real_number, interpolator_class, field_class, interp_neighbours, particles_adams_bashforth<partsize_t, real_number, 3, size_particle_rhs>, size_particle_rhs> computer;
+ particles_distr_mpi<partsize_t, real_number> particles_distr;
+
+ particles_adams_bashforth<partsize_t, real_number, 3, size_particle_rhs> positions_updater;
+
+ using computer_class = particles_field_computer<partsize_t, real_number, interpolator_class, interp_neighbours>;
+ computer_class computer;
+
+ field_class default_field;
std::unique_ptr<partsize_t[]> current_my_nb_particles_per_partition;
std::unique_ptr<partsize_t[]> current_offset_particles_for_partition;
@@ -52,8 +60,11 @@ public:
current_partition_interval({in_local_field_offset[IDX_Z], in_local_field_offset[IDX_Z] + in_local_field_dims[IDX_Z]}),
partition_interval_size(current_partition_interval.second - current_partition_interval.first),
interpolator(),
- computer(in_mpi_com, field_grid_dim, current_partition_interval,
- interpolator, in_field, in_spatial_box_width, in_spatial_box_offset, in_spatial_partition_width),
+ particles_distr(in_mpi_com, current_partition_interval,field_grid_dim),
+ positions_updater(),
+ computer(field_grid_dim, current_partition_interval,
+ interpolator, in_spatial_box_width, in_spatial_box_offset, in_spatial_partition_width),
+ default_field(in_field),
spatial_box_width(in_spatial_box_width), spatial_partition_width(in_spatial_partition_width),
my_spatial_low_limit(in_my_spatial_low_limit), my_spatial_up_limit(in_my_spatial_up_limit),
my_nb_particles(0), step_idx(in_current_iteration){
@@ -65,7 +76,7 @@ public:
~particles_system(){
}
- void init(abstract_particles_input<partsize_t, real_number>& particles_input){
+ void init(abstract_particles_input<partsize_t, real_number>& particles_input) {
TIMEZONE("particles_system::init");
my_particles_positions = particles_input.getMyParticles();
@@ -110,22 +121,65 @@ public:
void compute() final {
TIMEZONE("particles_system::compute");
- computer.compute_distr(current_my_nb_particles_per_partition.get(),
+ particles_distr.template compute_distr<computer_class, field_class, 3, size_particle_rhs>(
+ computer, default_field,
+ current_my_nb_particles_per_partition.get(),
my_particles_positions.get(),
my_particles_rhs.front().get(),
interp_neighbours);
}
+ template <class sample_field_class, int sample_size_particle_rhs>
+ void sample_compute(const sample_field_class& sample_field,
+ real_number sample_rhs[]) {
+ TIMEZONE("particles_system::compute");
+ particles_distr.template compute_distr<computer_class, sample_field_class, 3, sample_size_particle_rhs>(
+ computer, sample_field,
+ current_my_nb_particles_per_partition.get(),
+ my_particles_positions.get(),
+ sample_rhs,
+ interp_neighbours);
+ }
+
+ //- Not generic to enable sampling begin
+ void sample_compute_field(const field<float, FFTW, ONE>& sample_field,
+ real_number sample_rhs[]) final {
+ // sample_compute<decltype(sample_field), 1>(sample_field, sample_rhs);
+ }
+ void sample_compute_field(const field<float, FFTW, THREE>& sample_field,
+ real_number sample_rhs[]) final {
+ sample_compute<decltype(sample_field), 3>(sample_field, sample_rhs);
+ }
+ void sample_compute_field(const field<float, FFTW, THREExTHREE>& sample_field,
+ real_number sample_rhs[]) final {
+ sample_compute<decltype(sample_field), 9>(sample_field, sample_rhs);
+ }
+ void sample_compute_field(const field<double, FFTW, ONE>& sample_field,
+ real_number sample_rhs[]) final {
+ sample_compute<decltype(sample_field), 1>(sample_field, sample_rhs);
+ }
+ void sample_compute_field(const field<double, FFTW, THREE>& sample_field,
+ real_number sample_rhs[]) final {
+ sample_compute<decltype(sample_field), 3>(sample_field, sample_rhs);
+ }
+ void sample_compute_field(const field<double, FFTW, THREExTHREE>& sample_field,
+ real_number sample_rhs[]) final {
+ sample_compute<decltype(sample_field), 9>(sample_field, sample_rhs);
+ }
+ //- Not generic to enable sampling end
+
void move(const real_number dt) final {
TIMEZONE("particles_system::move");
- computer.move_particles(my_particles_positions.get(), my_nb_particles,
+ positions_updater.move_particles(my_particles_positions.get(), my_nb_particles,
my_particles_rhs.data(), std::min(step_idx,int(my_particles_rhs.size())),
dt);
}
void redistribute() final {
TIMEZONE("particles_system::redistribute");
- computer.redistribute(current_my_nb_particles_per_partition.get(),
+ particles_distr.template redistribute<computer_class, 3, size_particle_rhs, 1>(
+ computer,
+ current_my_nb_particles_per_partition.get(),
&my_nb_particles,
&my_particles_positions,
my_particles_rhs.data(), int(my_particles_rhs.size()),
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index a3c0d094..321f6d20 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -192,18 +192,21 @@ struct particles_system_build_container {
const particles_rnumber my_spatial_up_limit_z = particles_rnumber(local_field_offset[IDX_Z]+local_field_dims[IDX_Z])*spatial_partition_width[IDX_Z];
// Create the particles system
- particles_system<partsize_t, particles_rnumber, field_rnumber, field<field_rnumber, be, fc>, particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>, interpolation_size, ncomp(fc)>>* part_sys
- = new particles_system<partsize_t, particles_rnumber, field_rnumber, field<field_rnumber, be, fc>, particles_generic_interp<particles_rnumber, interpolation_size,spline_mode, interpolation_size, ncomp(fc)>>(field_grid_dim,
- spatial_box_width,
- spatial_box_offset,
- spatial_partition_width,
- my_spatial_low_limit_z,
- my_spatial_up_limit_z,
- local_field_dims,
- local_field_offset,
- (*fs_field),
- mpi_comm,
- in_current_iteration);
+ using particles_system_type = particles_system<partsize_t, particles_rnumber, field_rnumber,
+ field<field_rnumber, be, fc>,
+ particles_generic_interp<particles_rnumber, interpolation_size,spline_mode>,
+ interpolation_size, ncomp(fc)>;
+ particles_system_type* part_sys = new particles_system_type(field_grid_dim,
+ spatial_box_width,
+ spatial_box_offset,
+ spatial_partition_width,
+ my_spatial_low_limit_z,
+ my_spatial_up_limit_z,
+ local_field_dims,
+ local_field_offset,
+ (*fs_field),
+ mpi_comm,
+ in_current_iteration);
// Load particles from hdf5
particles_input_hdf5<partsize_t, particles_rnumber, 3,3> generator(mpi_comm, fname_input,
diff --git a/setup.py b/setup.py
index 139a4411..c531263a 100644
--- a/setup.py
+++ b/setup.py
@@ -121,7 +121,7 @@ src_file_list = ['full_code/direct_numerical_simulation',
'scope_timer']
particle_headers = [
- 'cpp/particles/abstract_particles_distr.hpp',
+ 'cpp/particles/particles_distr_mpi.hpp',
'cpp/particles/abstract_particles_input.hpp',
'cpp/particles/abstract_particles_output.hpp',
'cpp/particles/abstract_particles_system.hpp',
@@ -135,7 +135,9 @@ particle_headers = [
'cpp/particles/particles_system_builder.hpp',
'cpp/particles/particles_system.hpp',
'cpp/particles/particles_utils.hpp',
- 'cpp/particles/env_utils.hpp' ]
+ 'cpp/particles/particles_output_sampling_hdf5.hpp',
+ 'cpp/particles/particles_sampling.hpp',
+ 'cpp/particles/env_utils.hpp']
full_code_headers = ['cpp/full_code/main_code.hpp',
'cpp/full_code/codes_with_no_output.hpp',
--
GitLab