adds basic particle tracking for static field

checkpoint is somewhat improper because checkpointing functionality was
implemented in vorticity_equation instead of
direct_numerical_simulation, so that would be the next step.

preliminary test works with no obvious problems.

cpp/full_code/direct_numerical_simulation.hpp

@@ -62,6 +62,15 @@ class direct_numerical_simulation: public code_base
         int read_iteration(void);
         int write_iteration(void);
         int grow_file_datasets(void);
+
+        const inline std::string get_current_fname()
+        {
+            return (
+                    std::string(this->simname) +
+                    std::string("_checkpoint_") +
+                    std::to_string(this->checkpoint) +
+                    std::string(".h5"));
+        }
 };
 
 #endif//DIRECT_NUMERICAL_SIMULATION_HPP

cpp/full_code/static_field.cpp

@@ -74,7 +74,7 @@ int static_field<rnumber>::initialize(void)
     this->velocity->real_space_representation = false;
     //read vorticity field
     this->vorticity->io(
-	    this->simname + std::string("_checkpoint_0.h5"),
+	    this->get_current_fname(),
 	    "vorticity",
 	    this->iteration,
 	    true);
@@ -90,12 +90,44 @@ int static_field<rnumber>::initialize(void)
     this->velocity->ift();
     this->vorticity->ift();
 
+    // initialize particles
+    this->ps = particles_system_builder(
+                this->velocity,              // (field object)
+                this->kk,                     // (kspace object, contains dkx, dky, dkz)
+                tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
+                (long long int)nparticles,  // to check coherency between parameters and hdf input file
+	            this->get_current_fname(),    // particles input filename
+                std::string("/tracers0/state/") + std::to_string(this->iteration), // dataset name for initial input
+                std::string("/tracers0/rhs/")  + std::to_string(this->iteration),  // dataset name for initial input
+                tracers0_neighbours,        // parameter (interpolation no neighbours)
+                tracers0_smoothness,        // parameter
+                this->comm,
+                this->iteration+1);
+    // initialize output objects
+    this->particles_output_writer_mpi = new particles_output_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                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>(
+                MPI_COMM_WORLD,
+                this->ps->getGlobalNbParticles(),
+                (this->simname + "_particles.h5"),
+                "tracers0",
+                "position/0");
+    this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
+
     return EXIT_SUCCESS;
 }
 
 template <typename rnumber>
 int static_field<rnumber>::step(void)
 {
+    TIMEZONE("static_file::step");
+    this->ps->completeLoop(this->dt);
     this->iteration++;
     return EXIT_SUCCESS;
 }
@@ -103,6 +135,15 @@ int static_field<rnumber>::step(void)
 template <typename rnumber>
 int static_field<rnumber>::write_checkpoint(void)
 {
+    TIMEZONE("static_file::write_checkpoint");
+    this->particles_output_writer_mpi->open_file(this->get_current_fname());
+    this->particles_output_writer_mpi->template save<3>(
+            this->ps->getParticlesState(),
+            this->ps->getParticlesRhs(),
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->iteration);
+    this->particles_output_writer_mpi->close_file();
     return EXIT_SUCCESS;
 }
 
@@ -113,6 +154,9 @@ int static_field<rnumber>::finalize(void)
     if (this->myrank == 0)
         H5Fclose(this->stat_file);
     // good practice rule number n+1: always delete in inverse order of allocation
+    delete this->ps.release();
+    delete this->particles_output_writer_mpi;
+    delete this->particles_sample_writer_mpi;
     delete this->kk;
     delete this->velocity;
     delete this->vorticity;
@@ -126,6 +170,48 @@ int static_field<rnumber>::finalize(void)
 template <typename rnumber>
 int static_field<rnumber>::do_stats()
 {
+    TIMEZONE("static_field::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;
+
+    // allocate temporary data array
+    std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
+
+    /// copy position data
+
+    /// sample position
+    std::copy(this->ps->getParticlesState(),
+              this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
+              pdata.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "position",
+            this->ps->getParticlesState(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+
+    /// sample velocity
+    std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
+    this->ps->sample_compute_field(*this->velocity, pdata.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "velocity",
+            this->ps->getParticlesState(),
+            &pdata,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+
+    // deallocate temporary data array
+    delete[] pdata.release();
+
     return EXIT_SUCCESS;
 }
 
@@ -136,6 +222,14 @@ int static_field<rnumber>::read_parameters(void)
     this->direct_numerical_simulation::read_parameters();
     hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
     this->fftw_plan_rigor = hdf5_tools::read_string(parameter_file, "parameters/fftw_plan_rigor");
+    this->dt = hdf5_tools::read_value<double>(parameter_file, "parameters/dt");
+    this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+    this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+    this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+    this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+    this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
+    this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
+    this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
     H5Fclose(parameter_file);
     return EXIT_SUCCESS;
 }

cpp/full_code/static_field.hpp

@@ -34,6 +34,9 @@
 #include "kspace.hpp"
 #include "field.hpp"
 #include "full_code/direct_numerical_simulation.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
 
 template <typename rnumber>
 class static_field: public direct_numerical_simulation
@@ -41,13 +44,29 @@ class static_field: public direct_numerical_simulation
     public:
 
         /* parameters that are read in read_parameters */
+        double dt;
         std::string fftw_plan_rigor;
 
+        /* particle 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 */
         kspace<FFTW, SMOOTH> *kk;
         field<rnumber, FFTW, THREE> *vorticity;
         field<rnumber, FFTW, THREE> *velocity;
 
+        /* other stuff */
+        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;
+
 
         static_field(
                 const MPI_Comm COMMUNICATOR,