[WIP] adds tracer_with_collision_counter

broken, because I'm doing something wrong with the arrays pased to
p2p_distr_mpi::compute_distr.

CMakeLists.txt

@@ -301,6 +301,7 @@ set(cpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/field_tinterpolator.cpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/particle_set.cpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_rhs.cpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_with_collision_counter_rhs.cpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/particle_solver.cpp
     ${PROJECT_SOURCE_DIR}/cpp/scope_timer.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp
@@ -359,6 +360,7 @@ set(hpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/field_tinterpolator.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particle_rhs.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_rhs.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/particle_solver.hpp
     ${PROJECT_SOURCE_DIR}/cpp/scope_timer.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.hpp

cpp/full_code/test_tracer_set.cpp

@@ -28,6 +28,7 @@
 #include <random>
 #include "test_tracer_set.hpp"
 #include "particles/rhs/tracer_rhs.hpp"
+#include "particles/rhs/tracer_with_collision_counter_rhs.hpp"
 #include "particles/interpolation/particle_set.hpp"
 #include "particles/particle_solver.hpp"
 #include "particles/particles_input_random.hpp"
@@ -71,7 +72,7 @@ int test_tracer_set<rnumber>::do_work(void)
     vec_field->real_space_representation = true;
 
     field_tinterpolator<rnumber, FFTW, THREE, NONE> *fti = new field_tinterpolator<rnumber, FFTW, THREE, NONE>();
-    tracer_rhs<rnumber, FFTW, NONE> *trhs = new tracer_rhs<rnumber, FFTW, NONE>();
+    tracer_with_collision_counter_rhs<rnumber, FFTW, NONE> *trhs = new tracer_with_collision_counter_rhs<rnumber, FFTW, NONE>();
 
     particle_set<3, 2, 1> pset(
             vec_field->rlayout,

cpp/particles/abstract_particle_rhs.hpp

@@ -42,6 +42,11 @@ class abstract_particle_rhs
                 double t,
                 abstract_particle_set &pset,
                 particle_rnumber *result) const = 0;
+
+        // for post-timestep actions such as
+        // normalization of orientation vector
+        virtual int imposeModelConstraints(
+                abstract_particle_set &pset) const = 0;
 };
 
 #endif//ABSTRACT_PARTICLE_RHS_HPP

cpp/particles/interpolation/abstract_particle_set.hpp

@@ -31,6 +31,7 @@
 #include <vector>
 #include <memory>
 #include "field.hpp"
+#include "particles/p2p/p2p_distr_mpi.hpp"
 
 
 
@@ -56,8 +57,8 @@ class abstract_particle_set
         virtual ~abstract_particle_set(){}
 
         // extract particle information
-        virtual const particle_rnumber* getParticleState() const = 0;
-        virtual const partsize_t* getParticleIndex() const = 0;
+        virtual particle_rnumber* getParticleState() const = 0;
+        virtual partsize_t* getParticleIndex() const = 0;
         virtual int setParticleState(particle_rnumber *) = 0;
 
         virtual std::unique_ptr<particle_rnumber[]> extractFromParticleState(
@@ -66,8 +67,12 @@ class abstract_particle_set
 
         virtual partsize_t getLocalNumberOfParticles() const = 0;
         virtual partsize_t getTotalNumberOfParticles() const = 0;
+        virtual partsize_t* getParticlesPerPartition() const = 0;
         virtual const int getStateSize() const = 0;
 
+        // get p2p computer
+        virtual p2p_distr_mpi<partsize_t, particle_rnumber>* getP2PComputer() = 0;
+
         // generic loop over data
         template <class func_type>
         int LOOP(func_type expression)
@@ -125,6 +130,28 @@ class abstract_particle_set
                             FFTW,
                             THREExTHREE> &field_to_sample,
                 particle_rnumber *result) = 0;
+
+        // apply p2p kernel
+        template <int state_size,
+                  class p2p_kernel_class>
+        int applyP2PKernel(
+                p2p_kernel_class p2p_kernel,
+                std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data)
+        {
+            // there must be at least one array with additional data,
+            // since the p2p kernel expects an array where to store the result of the computation
+            assert(int(additional_data.size()) > int(0));
+            this->getP2PComputer()->template compute_distr<p2p_kernel_class,
+                                                           state_size, // I'd prefer to use this->getStateSize() here. not possible because virtual function is not constexpr...
+                                                           state_size>(
+                    p2p_kernel,
+                    this->getParticlesPerPartition(),
+                    this->getParticleState(),
+                    additional_data.data(),
+                    int(additional_data.size()),
+                    this->getParticleIndex());
+            return EXIT_SUCCESS;
+        }
 };
 
 #endif//ABSTRACT_PARTICLE_SET_HPP

cpp/particles/interpolation/particle_set.hpp

@@ -38,6 +38,7 @@
 #include "particles/interpolation/particles_field_computer.hpp"
 #include "particles/interpolation/particles_generic_interp.hpp"
 #include "particles/interpolation/abstract_particle_set.hpp"
+#include "particles/p2p/p2p_distr_mpi.hpp"
 
 
 
@@ -80,22 +81,24 @@ class particle_set: public abstract_particle_set
         std::vector<hsize_t> particle_file_layout;
 
         // raw objects for distributing data in MPI, and for interpolating fields
-        distributor_class   pdistributor;
+        distributor_class   pDistributor;
         particles_generic_interp<particle_rnumber, neighbours, smoothness> interpolation_formula;
-        interpolator_class  pinterpolator;
+        interpolator_class  pInterpolator;
+        p2p_distr_mpi<partsize_t, particle_rnumber> p2pComputer;
 
     public:
         particle_set(
                 const abstract_field_layout *fl,
                 const double dkx,
                 const double dky,
-                const double dkz):
+                const double dkz,
+                const double p2p_cutoff = 1.0):
             mpi_comm(fl->getMPIComm()),
             current_partition_interval(
                     {fl->getStart(IDXV_Z),
                      fl->getStart(IDXV_Z) + fl->getSubSize(IDXV_Z)}),
             partition_interval_size(current_partition_interval.second - current_partition_interval.first),
-            pdistributor(
+            pDistributor(
                     mpi_comm,
                     current_partition_interval,
                     std::array<size_t, 3>({  //field_grid_dim from particles_system_builder
@@ -104,7 +107,7 @@ class particle_set: public abstract_particle_set
                         fl->getSize(IDXV_Z),
                         })),
             interpolation_formula(),
-            pinterpolator(
+            pInterpolator(
                     std::array<size_t, 3>({ //field_grid_dim from particles_system_builder
                         fl->getSize(IDXV_X),
                         fl->getSize(IDXV_Y),
@@ -126,7 +129,27 @@ class particle_set: public abstract_particle_set
                         (4*acos(0) / dkx) / fl->getSize(IDXV_X),
                         (4*acos(0) / dky) / fl->getSize(IDXV_Y),
                         (4*acos(0) / dkz) / fl->getSize(IDXV_Z)
-                        }))
+                        })),
+            p2pComputer(
+                    mpi_comm,
+                    current_partition_interval,
+                    std::array<size_t, 3>({  //field_grid_dim from particles_system_builder
+                        fl->getSize(IDXV_X),
+                        fl->getSize(IDXV_Y),
+                        fl->getSize(IDXV_Z),
+                        }),
+                    std::array<particle_rnumber, 3>({ // spatial_box_width from particles_system_builder
+                        4*acos(0) / dkx,
+                        4*acos(0) / dky,
+                        4*acos(0) / dkz,
+                        }),
+                    std::array<particle_rnumber, 3>({ // spatial_box_offset from particles_system_builder
+                        0,
+                        0,
+                        0,
+                        }),
+                    p2p_cutoff)
+
         {
             // if these assertions fail,
             // it means the constructor prototype is broken
@@ -145,7 +168,7 @@ class particle_set: public abstract_particle_set
 
         ~particle_set(){}
 
-        const particle_rnumber* getParticleState() const
+        particle_rnumber* getParticleState() const
         {
             return this->local_state.get();
         }
@@ -159,7 +182,7 @@ class particle_set: public abstract_particle_set
             return EXIT_SUCCESS;
         }
 
-        const partsize_t* getParticleIndex() const
+        partsize_t* getParticleIndex() const
         {
             return this->local_index.get();
         }
@@ -174,6 +197,11 @@ class particle_set: public abstract_particle_set
             return this->total_number_of_particles;
         }
 
+        partsize_t* getParticlesPerPartition() const
+        {
+            return this->number_particles_per_partition.get();
+        }
+
         const int getStateSize() const
         {
             return state_size;
@@ -217,11 +245,11 @@ class particle_set: public abstract_particle_set
         {
             // attention: compute_distr adds result on top of existing values.
             // please clean up result as appropriate before call.
-            this->pdistributor.template compute_distr<interpolator_class,
+            this->pDistributor.template compute_distr<interpolator_class,
                                                       field<field_rnumber, be, fc>,
                                                       state_size,
                                                       ncomp(fc)>(
-                    this->pinterpolator,
+                    this->pInterpolator,
                     field_to_sample,
                     this->number_particles_per_partition.get(),
                     this->local_state.get(),
@@ -263,11 +291,11 @@ class particle_set: public abstract_particle_set
 
         int redistribute(std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data)
         {
-            this->pdistributor.template redistribute<interpolator_class,
+            this->pDistributor.template redistribute<interpolator_class,
                                                      state_size,
                                                      state_size,
                                                      1>(
-                    this->pinterpolator,
+                    this->pInterpolator,
                     this->number_particles_per_partition.get(),
                     &this->local_number_of_particles,
                     &this->local_state,
@@ -292,7 +320,7 @@ class particle_set: public abstract_particle_set
                     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);
+                        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;
                     },
@@ -305,12 +333,17 @@ class particle_set: public abstract_particle_set
 
         particle_rnumber getSpatialLowLimitZ()
         {
-            return this->pinterpolator.getSpatialLowLimitZ();
+            return this->pInterpolator.getSpatialLowLimitZ();
         }
 
         particle_rnumber getSpatialUpLimitZ()
         {
-            return this->pinterpolator.getSpatialUpLimitZ();
+            return this->pInterpolator.getSpatialUpLimitZ();
+        }
+
+        p2p_distr_mpi<partsize_t, particle_rnumber> *getP2PComputer()
+        {
+            return &(this->p2pComputer);
         }
 };
 

cpp/particles/particle_solver.cpp

@@ -51,6 +51,9 @@ int particle_solver::Euler(
         // deallocate temporary arrays
         delete[] rhs_val;
         delete[] xx;
+
+        // ensure new state is consistent with physical model
+        rhs.imposeModelConstraints(*(this->pset));
         return EXIT_SUCCESS;
     }
 

cpp/particles/rhs/tracer_rhs.hpp

@@ -58,6 +58,12 @@ class tracer_rhs: public abstract_particle_rhs
             std::fill_n(result, pset.getLocalNumberOfParticles()*3, 0);
             return (*(this->velocity))(t, pset, result);
         }
+
+        int imposeModelConstraints(
+                abstract_particle_set &pset) const
+        {
+            return EXIT_SUCCESS;
+        }
 };
 
 #endif//TRACER_RHS_HPP

cpp/particles/rhs/tracer_with_collision_counter_rhs.cpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2020 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 "particles/rhs/tracer_with_collision_counter_rhs.hpp"
+
+
+template class tracer_with_collision_counter_rhs<float, FFTW, NONE>;
+template class tracer_with_collision_counter_rhs<double, FFTW, NONE>;
+

cpp/particles/rhs/tracer_with_collision_counter_rhs.hpp

+/******************************************************************************
+*                                                                             *
+*  Copyright 2020 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 TRACER_WITH_COLLISION_COUNTER_RHS_HPP
+#define TRACER_WITH_COLLISION_COUNTER_RHS_HPP
+
+#include "particles/interpolation/field_tinterpolator.hpp"
+#include "particles/abstract_particle_rhs.hpp"
+#include "particles/rhs/tracer_rhs.hpp"
+#include "particles/p2p/p2p_ghost_collisions.hpp"
+
+
+template <typename rnumber,
+          field_backend be,
+          temporal_interpolation_type tt>
+class tracer_with_collision_counter_rhs: public tracer_rhs<rnumber, be, tt>
+{
+    protected:
+        p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t> p2p_gc;
+
+    public:
+        tracer_with_collision_counter_rhs(){}
+        ~tracer_with_collision_counter_rhs(){}
+
+        int operator()(
+                double t,
+                abstract_particle_set &pset,
+                abstract_particle_set::particle_rnumber *result) const
+        {
+            // interpolation adds on top of existing values, so result must be cleared.
+            std::fill_n(result, pset.getLocalNumberOfParticles()*3, 0);
+            int interpolation_result = (*(this->velocity))(t, pset, result);
+            std::vector<std::unique_ptr<abstract_particle_set::particle_rnumber[]>> bla;
+            bla.resize(1);
+            bla[0].reset(new abstract_particle_set::particle_rnumber[pset.getLocalNumberOfParticles()]);
+            pset.template applyP2PKernel<3, p2p_ghost_collisions<abstract_particle_set::particle_rnumber, abstract_particle_set::partsize_t>>(
+                    this->p2p_gc,
+                    bla);
+            return interpolation_result;
+        }
+};
+
+#endif//TRACER_WITH_COLLISION_COUNTER_RHS_HPP
+