change particle I/O file interaction

the "simname_particles.h5" file is now only opened once, and then
samples are written the same as before.
some sanity checks of the velocity and acceleration data should still be
performed.

bfps/cpp/full_code/NSVEparticles.cpp

@@ -2,7 +2,6 @@
 #include <cmath>
 #include "NSVEparticles.hpp"
 #include "scope_timer.hpp"
-#include "particles/particles_sampling.hpp"
 
 template <typename rnumber>
 int NSVEparticles<rnumber>::initialize(void)
@@ -27,6 +26,13 @@ int NSVEparticles<rnumber>::initialize(void)
                 "tracers0",
                 nparticles,
                 tracers0_integration_steps);
+    this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+        long long int, double, 3, 3>(
+                MPI_COMM_WORLD,
+                this->ps->getGlobalNbParticles(),
+                (this->simname + "_particles.h5"),
+                "tracers0",
+                "position/0");
     return EXIT_SUCCESS;
 }
 
@@ -60,6 +66,7 @@ int NSVEparticles<rnumber>::finalize(void)
 {
     this->ps.release();
     delete this->particles_output_writer_mpi;
+    delete this->particles_sample_writer_mpi;
     this->NSVE<rnumber>::finalize();
     return EXIT_SUCCESS;
 }
@@ -77,30 +84,65 @@ int NSVEparticles<rnumber>::do_stats()
     if (!(this->iteration % this->niter_part == 0))
         return EXIT_SUCCESS;
 
+    // allocate temporary data array
+    std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
+
+    // copy position data
+    std::copy(this->ps->getParticlesPositions(),
+              this->ps->getParticlesPositions()+this->ps->getLocalNbParticles(),
+              pdata.get());
+
     /// sample position
-    sample_particles_system_position(
-                                 this->ps,
-                                 (this->simname + "_particles.h5"), // filename
-                                 "tracers0",                        // hdf5 parent group
-                                 "position"                         // dataset basename TODO
-                                 );
+    this->particles_sample_writer_mpi->save_dataset(
+            "tracers0",
+            "position",
+            this->ps->getParticlesPositions(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx());
+
+    ////sample_particles_system_position(
+    ////                             this->ps,
+    ////                             (this->simname + "_particles.h5"), // filename
+    ////                             "tracers0",                        // hdf5 parent group
+    ////                             "position"                         // dataset basename TODO
+    ////                             );
 
     /// sample velocity
-    sample_from_particles_system(*this->tmp_vec_field,              // field to save
-                                 this->ps,
-                                 (this->simname + "_particles.h5"), // filename
-                                 "tracers0",                        // hdf5 parent group
-                                 "velocity"                         // dataset basename TODO
-                                 );
+    this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
+    this->particles_sample_writer_mpi->save_dataset(
+            "tracers0",
+            "velocity",
+            this->ps->getParticlesPositions(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx());
+    //sample_from_particles_system(*this->tmp_vec_field,              // field to save
+    //                             this->ps,
+    //                             (this->simname + "_particles.h5"), // filename
+    //                             "tracers0",                        // hdf5 parent group
+    //                             "velocity"                         // dataset basename TODO
+    //                             );
 
     /// compute acceleration and sample it
     this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field);
     this->tmp_vec_field->ift();
-    sample_from_particles_system(*this->tmp_vec_field,
-                                 this->ps,
-                                 (this->simname + "_particles.h5"),
-                                 "tracers0",
-                                 "acceleration");
+    this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
+    this->particles_sample_writer_mpi->save_dataset(
+            "tracers0",
+            "acceleration",
+            this->ps->getParticlesPositions(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx());
+    //sample_from_particles_system(*this->tmp_vec_field,
+    //                             this->ps,
+    //                             (this->simname + "_particles.h5"),
+    //                             "tracers0",
+    //                             "acceleration");
 
     return EXIT_SUCCESS;
 }

bfps/cpp/full_code/NSVEparticles.hpp

@@ -35,6 +35,7 @@
 #include "full_code/NSVE.hpp"
 #include "particles/particles_system_builder.hpp"
 #include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
 
 /** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
  *
@@ -58,6 +59,7 @@ class NSVEparticles: public NSVE<rnumber>
         /* other stuff */
         std::unique_ptr<abstract_particles_system<long long int, double>> ps;
         particles_output_hdf5<long long int, double,3,3> *particles_output_writer_mpi;
+        particles_output_sampling_hdf5<long long int, double, 3, 3> *particles_sample_writer_mpi;
 
 
         NSVEparticles(

bfps/cpp/full_code/field_test.cpp

@@ -24,7 +24,7 @@ int field_test<rnumber>::read_parameters()
     this->test::read_parameters();
     // in case any parameters are needed, this is where they should be read
     hid_t parameter_file;
-    hid_t dset, memtype, space;
+    hid_t dset;
     parameter_file = H5Fopen(
             (this->simname + std::string(".h5")).c_str(),
             H5F_ACC_RDONLY,

bfps/cpp/full_code/filter_test.cpp

@@ -40,7 +40,7 @@ int filter_test<rnumber>::read_parameters()
 {
     this->test::read_parameters();
     hid_t parameter_file;
-    hid_t dset, memtype, space;
+    hid_t dset;
     parameter_file = H5Fopen(
             (this->simname + std::string(".h5")).c_str(),
             H5F_ACC_RDONLY,

bfps/cpp/full_code/test.cpp

@@ -22,9 +22,8 @@ int test::main_loop(void)
 int test::read_parameters()
 {
     hid_t parameter_file;
-    hid_t dset, memtype, space;
+    hid_t dset;
     char fname[256];
-    char *string_data;
     sprintf(fname, "%s.h5", this->simname.c_str());
     parameter_file = H5Fopen(fname, H5F_ACC_RDONLY, H5P_DEFAULT);
     dset = H5Dopen(parameter_file, "/parameters/dealias_type", H5P_DEFAULT);

bfps/cpp/particles/abstract_particles_output.hpp

@@ -41,6 +41,10 @@ class abstract_particles_output {
     partsize_t particles_chunk_current_offset;
 
 protected:
+    MPI_Comm& getCom(){
+        return mpi_com;
+    }
+
     MPI_Comm& getComWriter(){
         return mpi_com_writer;
     }

bfps/cpp/particles/particles_output_sampling_hdf5.hpp

@@ -9,10 +9,11 @@ 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>{
+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,
@@ -20,7 +21,7 @@ class particles_output_sampling_hdf5 : public abstract_particles_output<partsize
 
     hid_t file_id, pgroup_id;
 
-    const std::string dataset_name;
+    std::string dataset_name;
     const bool use_collective_io;
 
 public:
@@ -34,7 +35,6 @@ public:
         int dataset_exists = -1;
 
         if(my_rank == 0){
-            // Parallel HDF5 write
             hid_t file_id = H5Fopen(
                     in_filename.c_str(),
                     H5F_ACC_RDWR | H5F_ACC_DEBUG,
@@ -54,16 +54,18 @@ public:
         return dataset_exists;
     }
 
-    particles_output_sampling_hdf5(MPI_Comm in_mpi_com,
-                          const partsize_t inTotalNbParticles,
-                                   const std::string& in_filename,
-                                   const std::string& in_groupname,
-                          const std::string& in_dataset_name,
-                          const bool in_use_collective_io = false)
+    particles_output_sampling_hdf5(
+            MPI_Comm in_mpi_com,
+            const partsize_t inTotalNbParticles,
+            const std::string& in_filename,
+            const std::string& in_groupname,
+            const std::string& in_dataset_name,
+            const bool in_use_collective_io = false)
             : Parent(in_mpi_com, inTotalNbParticles, 1),
               dataset_name(in_dataset_name),
               use_collective_io(in_use_collective_io){
         if(Parent::isInvolved()){
+            // prepare parallel MPI access property list
             hid_t plist_id_par = H5Pcreate(H5P_FILE_ACCESS);
             assert(plist_id_par >= 0);
             int retTest = H5Pset_fapl_mpio(
@@ -72,7 +74,7 @@ public:
                     MPI_INFO_NULL);
             assert(retTest >= 0);
 
-            // Parallel HDF5 write
+            // open file for parallel HDF5 access
             file_id = H5Fopen(
                     in_filename.c_str(),
                     H5F_ACC_RDWR | H5F_ACC_DEBUG,
@@ -81,6 +83,7 @@ public:
             retTest = H5Pclose(plist_id_par);
             assert(retTest >= 0);
 
+            // open group
             pgroup_id = H5Gopen(
                     file_id,
                     in_groupname.c_str(),
@@ -91,13 +94,65 @@ public:
 
     ~particles_output_sampling_hdf5(){
         if(Parent::isInvolved()){
+            // close group
             int retTest = H5Gclose(pgroup_id);
             assert(retTest >= 0);
+            // close file
             retTest = H5Fclose(file_id);
             assert(retTest >= 0);
         }
     }
 
+    int switch_to_group(
+            const std::string &in_groupname)
+    {
+        if(Parent::isInvolved()){
+            // close old group
+            int retTest = H5Gclose(pgroup_id);
+            assert(retTest >= 0);
+
+            // open new group
+            pgroup_id = H5Gopen(
+                    file_id,
+                    in_groupname.c_str(),
+                    H5P_DEFAULT);
+            assert(pgroup_id >= 0);
+        }
+        return EXIT_SUCCESS;
+    }
+
+    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 partsize_t index_particles[],
+            const partsize_t nb_particles,
+            const int idx_time_step)
+    {
+        // update group
+        int retTest = this->switch_to_group(
+                in_groupname);
+        assert(retTest == EXIT_SUCCESS);
+        // update dataset name
+        dataset_name = in_dataset_name + "/" + std::to_string(idx_time_step);
+        int dataset_exists;
+        if (this->getMyRank() == 0)
+            dataset_exists = H5Lexists(
+                pgroup_id,
+                dataset_name.c_str(),
+                H5P_DEFAULT);
+        AssertMpi(MPI_Bcast(&dataset_exists, 1, MPI_INT, 0, this->getCom()));
+        if (dataset_exists == 0)
+            this->save(
+                input_particles_positions,
+                input_particles_rhs,
+                index_particles,
+                nb_particles,
+                idx_time_step);
+        return EXIT_SUCCESS;
+    }
+
     void write(
             const int /*idx_time_step*/,
             const real_number* /*particles_positions*/,
@@ -108,18 +163,26 @@ public:
 
         TIMEZONE("particles_output_hdf5::write");
 
-        assert(particles_idx_offset < Parent::getTotalNbParticles() || (particles_idx_offset == Parent::getTotalNbParticles() && nb_particles == 0));
+        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);
+        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);
+            int rethdf = H5Pset_dxpl_mpio(
+                    plist_id,
+                    (use_collective_io ?
+                     H5FD_MPIO_COLLECTIVE :
+                     H5FD_MPIO_INDEPENDENT));
             assert(rethdf >= 0);
         }
         {

setup.py

@@ -88,7 +88,8 @@ print('This is bfps version ' + VERSION)
 
 
 ### lists of files and MANIFEST.in
-src_file_list = ['full_code/joint_acc_vel_stats',
+src_file_list = ['full_code/NSVEparticles',
+                 'full_code/joint_acc_vel_stats',
                  'full_code/test',
                  'full_code/filter_test',
                  'full_code/field_test',
@@ -127,8 +128,7 @@ src_file_list = ['full_code/joint_acc_vel_stats',
                  'spline_n9',
                  'spline_n10',
                  'Lagrange_polys',
-                 'scope_timer',
-                 'full_code/NSVEparticles']
+                 'scope_timer']
 
 particle_headers = [
         'cpp/particles/particles_distr_mpi.hpp',