Add first version of the P2P collision counter

cpp/full_code/static_field_with_ghost_collisions.cpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
+*                                                                             *
+*  This file is part of bfps.                                                 *
+*                                                                             *
+*  bfps 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,                    *
+*  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/>               *
+*                                                                             *
+* Contact: Cristian.Lalescu@ds.mpg.de                                         *
+*                                                                             *
+******************************************************************************/
+
+
+#include <string>
+#include <cmath>
+#include "static_field_with_ghost_collisions.hpp"
+#include "scope_timer.hpp"
+#include "fftw_tools.hpp"
+
+
+template <typename rnumber>
+int static_field<rnumber>::initialize(void)
+{
+    TIMEZONE("static_file::initialize");
+    this->read_iteration();
+    this->read_parameters();
+    if (this->myrank == 0)
+    {
+        // set caching parameters
+        hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
+        herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
+        variable_used_only_in_assert(cache_err);
+        DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
+        this->stat_file = H5Fopen(
+                (this->simname + ".h5").c_str(),
+                H5F_ACC_RDWR,
+                fapl);
+    }
+    int data_file_problem;
+    if (this->myrank == 0)
+        data_file_problem = this->grow_file_datasets();
+    MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
+    if (data_file_problem > 0)
+    {
+        std::cerr <<
+            data_file_problem <<
+            " problems growing file datasets.\ntrying to exit now." <<
+            std::endl;
+        return EXIT_FAILURE;
+    }
+
+    this->vorticity = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+            fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+    this->vorticity->real_space_representation = false;
+
+    this->velocity = new field<rnumber, FFTW, THREE>(
+            nx, ny, nz,
+            this->comm,
+            fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+    this->velocity->real_space_representation = false;
+    //read vorticity field
+    this->vorticity->io(
+            this->get_current_fname(),
+            "vorticity",
+            this->iteration,
+            true);
+    this->kk = new kspace<FFTW, SMOOTH>(
+            this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+
+    // compute the velocity field and store
+    invert_curl(
+        this->kk,
+        this->vorticity,
+        velocity);
+    // transform velocity and vorticity fields to real space
+    this->velocity->ift();
+    this->vorticity->ift();
+
+    p2p_ghost_collisions<double, long long int> p2p_ghost_collisions;
+    p2p_ghost_collisions.setEnable(this->enable_p2p);
+
+    this->ps = Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<long long int, double>>,
+                       int, 1, 11, 1, // interpolation_size
+                       int, 0, 3, 1, // spline_mode
+                       particles_system_build_container<long long int, rnumber,FFTW,THREE,double,
+                                                        collision_detector<double, long long int>,
+                                                        particles_inner_computer_empty<double,long long int>,
+                                                        3,3>>(
+                           this->tracers0_neighbours, // template iterator 1
+                           this->tracers0_smoothness, // template iterator 2
+                           this->velocity,
+                           this->kk,
+                           tracers0_integration_steps,
+                           (long long int)nparticles,
+                           this->get_current_fname(),
+                           std::string("/tracers0/state/") + std::to_string(this->iteration),
+                           std::string("/tracers0/rhs/") + std::to_string(this->iteration),
+                           this->comm,
+                           this->iteration+1,
+                           std::move(cd),
+                           std::move(p2p_ghost_collisions),
+                           this->tracers0_cutoff);
+
+    // 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->getP2PComputer().reset_collision_counter(); // TODO should it be here ?
+    this->ps->completeLoop(this->dt);
+    const long int nbCollision = this->ps->getP2PComputer().get_collision_counter(); // TODO should it be here ?
+    this->iteration++;
+    return EXIT_SUCCESS;
+}
+
+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;
+}
+
+template <typename rnumber>
+int static_field<rnumber>::finalize(void)
+{
+    TIMEZONE("static_field::finalize");
+    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;
+    return EXIT_SUCCESS;
+}
+
+/** \brief Compute statistics.
+ *
+ */
+
+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;
+}
+
+template <typename rnumber>
+int static_field<rnumber>::read_parameters(void)
+{
+    TIMEZONE("static_field::read_parameters");
+    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;
+}
+
+template class static_field<float>;
+template class static_field<double>;
+
+

cpp/full_code/static_field_with_ghost_collisions.hpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 Max Planck Institute                                *
+*                 for Dynamics and Self-Organization                  *
+*                                                                     *
+*  This file is part of bfps.                                         *
+*                                                                     *
+*  bfps 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,            *
+*  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/>       *
+*                                                                     *
+* Contact: Cristian.Lalescu@ds.mpg.de                                 *
+*                                                                     *
+**********************************************************************/
+
+
+#ifndef STATIC_FIELD_WITH_GHOST_COLLISIONS_HPP
+#define STATIC_FIELD_WITH_GHOST_COLLISIONS_HPP
+
+#include <cstdlib>
+#include "base.hpp"
+#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"
+#include "particles/p2p_ghost_collisions.hpp"
+
+template <typename rnumber>
+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_with_p2p<long long int, double, p2p_ghost_collisions>> 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,
+                const std::string &simulation_name):
+            direct_numerical_simulation(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~static_field(){}
+
+        int initialize(void);
+        int step(void);
+        int finalize(void);
+
+        virtual int read_parameters(void);
+        int write_checkpoint(void);
+        int do_stats(void);
+};
+
+#endif//STATIC_FIELD_HPP
+

cpp/particles/abstract_particles_system_with_p2p.hpp

+#ifndef ABSTRACT_PARTICLES_SYSTEM_WITH_P2P_HPP
+#define ABSTRACT_PARTICLES_SYSTEM_WITH_P2P_HPP
+
+#include "abstract_particles_system.hpp"
+
+template <class partsize_t, class real_number, class p2p_computer_class>
+class abstract_particles_system_with_p2p : public abstract_particles_system<partsize_t,real_number> {
+public:
+    virtual p2p_computer_class& getP2PComputer() = 0;
+    virtual const p2p_computer_class& getP2PComputer() const = 0;
+};
+
+#endif // ABSTRACT_PARTICLES_SYSTEM_WITH_P2P_HPP

cpp/particles/p2p_computer.hpp

@@ -98,6 +98,9 @@ public:
         rhs_part2[3+IDXC_Z] += pos_part1[3+IDXC_Z] * 4 * ww * dot_product;
     }
 
+
+    void merge(const p2p_computer&){}
+
     bool isEnable() const {
         return isActive;
     }

cpp/particles/p2p_computer_empty.hpp

@@ -46,6 +46,8 @@ public:
                              const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
     }
 
+    void merge(const p2p_computer_empty&){}
+
     constexpr static bool isEnable() {
         return false;
     }

cpp/particles/p2p_cylinder_collisions.hpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
+*                                                                             *
+*  This file is part of bfps.                                                 *
+*                                                                             *
+*  bfps 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,                    *
+*  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/>               *
+*                                                                             *
+* Contact: Cristian.Lalescu@ds.mpg.de                                         *
+*                                                                             *
+******************************************************************************/
+#ifndef P2P_CYLINDER_COLLISIONS_HPP
+#define P2P_CYLINDER_COLLISIONS_HPP
+
+#include <cstring>
+
+template <class real_number, class partsize_t>
+class p2p_cylinder_collisions{
+    long int collision_counter;
+
+public:
+    p2p_cylinder_collisions() : collision_counter(0){}
+
+    // Copy constructor use a counter set to zero
+    p2p_cylinder_collisions(const p2p_cylinder_collisions&) : collision_counter(0){}
+
+    template <int size_particle_rhs>
+    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
+    }
+
+    template <int size_particle_rhs>
+    void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
+    }
+
+    template <int size_particle_positions, int size_particle_rhs>
+    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
+                             const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
+                             const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
+        // TODO check the orientation, etc....
+        if(check the test here){
+            collision_counter += 1;
+        }
+    }
+
+    void merge(const p2p_cylinder_collisions& other){
+        collision_counter += other.collision_counter;
+    }
+
+    constexpr static bool isEnable() {
+        return false;
+    }
+
+    long int get_collision_counter() const{
+        return collision_counter;
+    }
+
+    void reset_collision_counter() const{
+        collision_counter = 0;
+    }
+};
+
+
+#endif // P2P_CYLINDER_COLLISIONS_HPP

cpp/particles/p2p_ghost_collisions.hpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
+*                                                                             *
+*  This file is part of bfps.                                                 *
+*                                                                             *
+*  bfps 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,                    *
+*  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/>               *
+*                                                                             *
+* Contact: Cristian.Lalescu@ds.mpg.de                                         *
+*                                                                             *
+******************************************************************************/
+#ifndef P2P_GHOST_COLLISIONS_HPP
+#define P2P_GHOST_COLLISIONS_HPP
+
+#include <cstring>
+
+template <class real_number, class partsize_t>
+class p2p_ghost_collisions{
+    long int collision_counter;
+
+public:
+    p2p_ghost_collisions() : collision_counter(0){}
+
+    // Copy constructor use a counter set to zero
+    p2p_ghost_collisions(const p2p_ghost_collisions&) : collision_counter(0){}
+
+    template <int size_particle_rhs>
+    void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
+    }
+
+    template <int size_particle_rhs>
+    void reduce_particles_rhs(real_number /*rhs_dst*/[], const real_number /*rhs_src*/[], const partsize_t /*nbParticles*/) const{
+    }
+
+    template <int size_particle_positions, int size_particle_rhs>
+    void compute_interaction(const real_number /*pos_part1*/[], real_number /*rhs_part1*/[],
+                             const real_number /*pos_part2*/[], real_number /*rhs_part2*/[],
+                             const real_number /*dist_pow2*/,  const real_number /*cutoff*/,
+                             const real_number /*xshift_coef*/, const real_number /*yshift_coef*/, const real_number /*zshift_coef*/) const{
+        collision_counter += 1;
+    }
+
+    void merge(const p2p_ghost_collisions& other){
+        collision_counter += other.collision_counter;
+    }
+
+    constexpr static bool isEnable() {
+        return false;
+    }
+
+    long int get_collision_counter() const{
+        return collision_counter;
+    }
+
+    void reset_collision_counter() const{
+        collision_counter = 0;
+    }
+};
+
+
+#endif // P2P_GHOST_COLLISIONS_HPP

cpp/particles/particles_system.hpp

@@ -29,6 +29,7 @@
 #include <array>
 
 #include "abstract_particles_system.hpp"
+#include "abstract_particles_system_with_p2p.hpp"
 #include "particles_distr_mpi.hpp"
 #include "particles_output_hdf5.hpp"
 #include "particles_output_mpiio.hpp"
@@ -41,7 +42,7 @@
 
 template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours,
           int size_particle_positions, int size_particle_rhs, class p2p_computer_class, class particles_inner_computer_class>
-class particles_system : public abstract_particles_system<partsize_t, real_number> {
+class particles_system : public abstract_particles_system_with_p2p<partsize_t, real_number, p2p_computer_class> {
     static_assert(size_particle_positions >= 3, "There should be at least the positions X,Y,Z in the state");
 
     MPI_Comm mpi_com;
@@ -402,6 +403,14 @@ public:
             }
         }
     }
+
+    const p2p_computer_class& getP2PComputer() const final{
+        return computer_p2p;
+    }
+
+    p2p_computer_class& getP2PComputer() final{
+        return computer_p2p;
+    }
 };
 
 

cpp/particles/particles_system_builder.hpp

@@ -140,7 +140,7 @@ template <class partsize_t, class field_rnumber, field_backend be, field_compone
           class particles_inner_computer_class, int size_particle_positions, int size_particle_rhs>
 struct particles_system_build_container {
     template <const int interpolation_size, const int spline_mode>
-    static std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> instanciate(
+    static std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>> instanciate(
              const field<field_rnumber, be, fc>* fs_field, // (field object)
              const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
              const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
@@ -272,7 +272,7 @@ struct particles_system_build_container {
         part_sys->setParticleFileLayout(generator.getParticleFileLayout());
 
         // Return the created particles system
-        return std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>(part_sys);
+        return std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>>(part_sys);
     }
 };
 
@@ -305,7 +305,7 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
 template <class partsize_t, class field_rnumber, field_backend be, field_components fc,
           class p2p_computer_class, class particles_inner_computer_class,
           class particles_rnumber = double>
-inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>> particles_system_builder_with_p2p(
+inline std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>> particles_system_builder_with_p2p(
         const field<field_rnumber, be, fc>* fs_field, // (field object)
         const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
         const int nsteps, // to check coherency between parameters and hdf input file (nb rhs)
@@ -319,7 +319,7 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
         p2p_computer_class p2p_computer,
         particles_inner_computer_class inner_computer,
         const particles_rnumber cutoff){
-    return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
+    return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>>,
                        int, 1, 11, 1, // interpolation_size
                        int, 0, 3, 1, // spline_mode
                        particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,