From 7d3bb09554c1981c7851d64f2cc06c502743fd1b Mon Sep 17 00:00:00 2001
From: Berenger Bramas <berenger.bramas@mpcdf.mpg.de>
Date: Thu, 12 Oct 2017 17:07:27 +0200
Subject: [PATCH] Update the kernels and make it possible to update parameters
 at runtime

---
 bfps/cpp/field.hpp                            |  4 +
 bfps/cpp/full_code/NSVEparticlesP2P.cpp       | 25 ++++--
 bfps/cpp/full_code/NSVEparticlesP2P.hpp       | 10 ++-
 .../particles/abstract_particles_system.hpp   | 16 ++++
 bfps/cpp/particles/p2p_computer.hpp           | 84 ++++++++++---------
 bfps/cpp/particles/p2p_computer_empty.hpp     | 14 ++--
 bfps/cpp/particles/p2p_distr_mpi.hpp          | 49 ++---------
 .../particles/particles_adams_bashforth.hpp   | 29 +++++--
 bfps/cpp/particles/particles_distr_mpi.hpp    | 72 +---------------
 .../particles/particles_field_computer.hpp    |  6 +-
 .../particles/particles_inner_computer.hpp    | 65 ++++++++++++++
 .../particles_inner_computer_empty.hpp        | 24 ++++++
 bfps/cpp/particles/particles_system.hpp       | 79 +++++++++++------
 .../particles/particles_system_builder.hpp    | 43 +++++++---
 14 files changed, 313 insertions(+), 207 deletions(-)
 create mode 100644 bfps/cpp/particles/particles_inner_computer.hpp
 create mode 100644 bfps/cpp/particles/particles_inner_computer_empty.hpp

diff --git a/bfps/cpp/field.hpp b/bfps/cpp/field.hpp
index 9a5ab1be..2d93a0c9 100644
--- a/bfps/cpp/field.hpp
+++ b/bfps/cpp/field.hpp
@@ -58,6 +58,10 @@ class field
     private:
         rnumber *__restrict__ data; /**< data array */
     public:
+        constexpr int nb_components() {
+            return  ncomp(fc);
+        }
+
         hsize_t npoints; /**< total number of grid points. Useful for normalization. */
         bool real_space_representation; /**< `true` if field is in real space representation. */
 
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.cpp b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
index 08326119..9574b226 100644
--- a/bfps/cpp/full_code/NSVEparticlesP2P.cpp
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.cpp
@@ -4,13 +4,20 @@
 #include "scope_timer.hpp"
 #include "particles/particles_sampling.hpp"
 #include "particles/p2p_computer.hpp"
+#include "particles/particles_inner_computer.hpp"
 
 template <typename rnumber>
 int NSVEparticlesP2P<rnumber>::initialize(void)
 {
     this->NSVE<rnumber>::initialize();
 
-    this->ps = particles_system_builder(
+    p2p_computer<rnumber, long long int> current_p2p_computer;
+    current_p2p_computer.setEnable(enable_p2p);
+
+    particles_inner_computer<rnumber, long long int> current_particles_inner_computer(inner_v0);
+    current_particles_inner_computer.setEnable(enable_inner);
+
+    this->ps = particles_system_builder_with_p2p(
                 this->fs->cvelocity,              // (field object)
                 this->fs->kk,                     // (kspace object, contains dkx, dky, dkz)
                 tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
@@ -21,10 +28,13 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
                 tracers0_neighbours,        // parameter (interpolation no neighbours)
                 tracers0_smoothness,        // parameter
                 this->comm,
-                this->fs->iteration+1);
-    // TODO P2P write particle data too
+                this->fs->iteration+1,
+                std::move(current_p2p_computer),
+                std::move(current_particles_inner_computer),
+                cutoff);
+
     this->particles_output_writer_mpi = new particles_output_hdf5<
-        long long int, double, 3, 3>(
+        long long int, double, 6, 6>(
                 MPI_COMM_WORLD,
                 "tracers0",
                 nparticles,
@@ -37,7 +47,12 @@ int NSVEparticlesP2P<rnumber>::step(void)
 {
     this->fs->compute_velocity(this->fs->cvorticity);
     this->fs->cvelocity->ift();
-    this->ps->completeLoop(this->dt);
+    if(enable_vorticity_omega){
+        this->ps->completeLoopWithVorticity(this->dt, *this->fs->cvorticity);
+    }
+    else{
+        this->ps->completeLoop(this->dt);
+    }
     this->NSVE<rnumber>::step();
     return EXIT_SUCCESS;
 }
diff --git a/bfps/cpp/full_code/NSVEparticlesP2P.hpp b/bfps/cpp/full_code/NSVEparticlesP2P.hpp
index 9b015659..73427bca 100644
--- a/bfps/cpp/full_code/NSVEparticlesP2P.hpp
+++ b/bfps/cpp/full_code/NSVEparticlesP2P.hpp
@@ -55,10 +55,15 @@ class NSVEparticlesP2P: public NSVE<rnumber>
         int tracers0_neighbours;
         int tracers0_smoothness;
 
+        double cutoff;
+        bool enable_p2p;
+        bool enable_inner;
+        bool enable_vorticity_omega;
+
         /* other stuff */
         std::unique_ptr<abstract_particles_system<long long int, double>> ps;
         // TODO P2P use a reader with particle data
-        particles_output_hdf5<long long int, double,3,3> *particles_output_writer_mpi;
+        particles_output_hdf5<long long int, double,6,6> *particles_output_writer_mpi;
 
 
         NSVEparticlesP2P(
@@ -66,7 +71,8 @@ class NSVEparticlesP2P: public NSVE<rnumber>
                 const std::string &simulation_name):
             NSVE<rnumber>(
                     COMMUNICATOR,
-                    simulation_name){}
+                    simulation_name),
+            cutoff(std::numeric_limits<double>::max()){}
         ~NSVEparticlesP2P(){}
 
         int initialize(void);
diff --git a/bfps/cpp/particles/abstract_particles_system.hpp b/bfps/cpp/particles/abstract_particles_system.hpp
index 26ce7198..b2d09566 100644
--- a/bfps/cpp/particles/abstract_particles_system.hpp
+++ b/bfps/cpp/particles/abstract_particles_system.hpp
@@ -16,6 +16,10 @@ public:
 
     virtual void compute_p2p() = 0;
 
+    virtual void compute_particles_inner() = 0;
+
+    virtual void compute_particles_inner(const real_number particle_extra_rhs[]) = 0;
+
     virtual void move(const real_number dt) = 0;
 
     virtual void redistribute() = 0;
@@ -26,6 +30,9 @@ public:
 
     virtual void completeLoop(const real_number dt) = 0;
 
+    virtual void completeLoopWithVorticity(const real_number dt,
+                                           const real_number particle_extra_rhs[]) = 0;
+
     virtual const real_number* getParticlesPositions() const = 0;
 
     virtual const std::unique_ptr<real_number[]>* getParticlesRhs() const = 0;
@@ -54,6 +61,15 @@ public:
     virtual void sample_compute_field(const field<double, FFTW, THREExTHREE>& sample_field,
                                 real_number sample_rhs[]) = 0;
     //- Not generic to enable sampling end
+
+    template <typename rnumber, field_backend be, field_components fc>
+    void completeLoopWithVorticity(const real_number dt,
+                                   const field<rnumber, be, fc>& in_field) {
+        static_assert(fc == THREE, "only THREE is supported for now");
+        std::unique_ptr<real_number[]> extra_rhs(new real_number[getLocalNbParticles()*3]());
+        sample_compute_field(in_field, extra_rhs.get());
+        completeLoopWithVorticity(dt, extra_rhs.get());
+    }
 };
 
 #endif
diff --git a/bfps/cpp/particles/p2p_computer.hpp b/bfps/cpp/particles/p2p_computer.hpp
index bd6a4001..46328c3a 100644
--- a/bfps/cpp/particles/p2p_computer.hpp
+++ b/bfps/cpp/particles/p2p_computer.hpp
@@ -2,27 +2,31 @@
 #define P2P_COMPUTER_HPP
 
 #include <cstring>
-
-/** \brief A simple distance weighting function.
- *
- *  This function returns 1 if a distance is smaller than a cut-off length,
- *  i.e. particle 1 interacts with particle 2 if particle 2 is inside a
- *  sphere of radius `cutoff' centered on particle 1.
- */
-double dumb_distance_weight(
-        double distance,
-        double cutoff)
-{
-    // this function should only be called for interacting particles,
-    // and particles interact if they are closer than cutoff.
-    assert(distance < cutoff);
-    return 1.0;
-}
+#include <cassert>
 
 template <class real_number, class partsize_t>
 class p2p_computer{
+
+    bool isActive;
+
+    /** \brief A simple distance weighting function.
+     *
+     *  This function returns 1 if a distance is smaller than a cut-off length,
+     *  i.e. particle 1 interacts with particle 2 if particle 2 is inside a
+     *  sphere of radius `cutoff' centered on particle 1.
+     */
+    static double dumb_distance_weight(
+            const double dist_pow2,
+            const double cutoff){
+        // this function should only be called for interacting particles,
+        // and particles interact if they are closer than cutoff.
+        assert(dist_pow2 < cutoff*cutoff);
+        return 1.0;
+    }
+
+
 public:
-    constexpr static int size_data = 3;
+    p2p_computer() : isActive(true){}
 
     template <int size_particle_rhs>
     void init_result_array(real_number rhs[], const partsize_t nbParticles) const{
@@ -31,25 +35,26 @@ public:
 
     template <int size_particle_rhs>
     void reduce_particles_rhs(real_number rhs_dst[], const real_number rhs_src[], const partsize_t nbParticles) const{
+        static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
         for(int idx_part = 0 ; idx_part < nbParticles ; ++idx_part){
-            for(int idx_rhs = 0 ; idx_rhs < size_particle_rhs ; ++idx_rhs){
+            // We merge only the values modified by the current kernel (3-5)
+            for(int idx_rhs = 3 ; idx_rhs < size_particle_rhs ; ++idx_rhs){
                 rhs_dst[idx_part*size_particle_rhs+idx_rhs] += rhs_src[idx_part*size_particle_rhs+idx_rhs];
             }
         }
     }
 
-    template <int size_particle_positions, int size_particle_data, int size_particle_rhs>
-    void compute_interaction(const real_number pos_part1[], const real_number data_part1[], real_number rhs_part1[],
-                             const real_number pos_part2[], const real_number data_part2[], real_number rhs_part2[],
-                             const real_number dist_pow2,
-                             const real_number xshift_coef, const real_number yshift_coef, const real_number zshift_coef) const{
-        static_assert(size_particle_positions == 3, "This kernel works only with 3 values for one particle's position");
-        static_assert(size_particle_rhs >= 6, "This kernel works only with more than 6 values per particle's rhs");
-        static_assert(size_particle_data == 3, "This kernel works only with 3 values per particle's' data");
+    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{
+        static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position+orientation");
+        static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
 
         // TODO: a reasonable way of choosing between different distance_weight functions should be thought of.
         // We need to ask Michael about how flexible this distance_weight needs to be.
-        const double ww = dumb_distance_weight(distance, max_distance);
+        const double ww = dumb_distance_weight(dist_pow2, cutoff);
         ///
         /// term in equation is:
         ///
@@ -57,20 +62,23 @@ public:
         ///     (4 / \tau) \sum_j W_\ell ( | x^i - x^j | ) (p^i \cdot p^j)p^j
         /// \f]
         ///
-        const double dot_product = (data_part1[0]*data_part2[0] +
-                              data_part1[1]*data_part2[1] +
-                              data_part1[2]*data_part2[2]);
-        rhs_part1[3] += data_part2[0] * 4 * ww * dot_product;
-        rhs_part1[4] += data_part2[1] * 4 * ww * dot_product;
-        rhs_part1[5] += data_part2[2] * 4 * ww * dot_product;
-        rhs_part2[3] += data_part1[0] * 4 * ww * dot_product;
-        rhs_part2[4] += data_part1[1] * 4 * ww * dot_product;
-        rhs_part2[5] += data_part1[2] * 4 * ww * dot_product;
+        const double dot_product = (pos_part1[3+IDX_X]*pos_part2[3+IDX_X] +
+                              pos_part1[3+IDX_Y]*pos_part2[3+IDX_Y] +
+                              pos_part1[3+IDX_Z]*pos_part2[3+IDX_Z]);
+        rhs_part1[3+IDX_X] += pos_part2[3+IDX_X] * 4 * ww * dot_product;
+        rhs_part1[3+IDX_Y] += pos_part2[3+IDX_Y] * 4 * ww * dot_product;
+        rhs_part1[3+IDX_Z] += pos_part2[3+IDX_Z] * 4 * ww * dot_product;
+        rhs_part2[3+IDX_X] += pos_part1[3+IDX_X] * 4 * ww * dot_product;
+        rhs_part2[3+IDX_Y] += pos_part1[3+IDX_Y] * 4 * ww * dot_product;
+        rhs_part2[3+IDX_Z] += pos_part1[3+IDX_Z] * 4 * ww * dot_product;
     }
 
+    bool isEnable() const {
+        return isActive;
+    }
 
-    constexpr static bool isEmpty() {
-        return false;
+    void setEnable(const bool inIsActive) {
+        isActive = inIsActive;
     }
 };
 
diff --git a/bfps/cpp/particles/p2p_computer_empty.hpp b/bfps/cpp/particles/p2p_computer_empty.hpp
index 7076061e..5e442b86 100644
--- a/bfps/cpp/particles/p2p_computer_empty.hpp
+++ b/bfps/cpp/particles/p2p_computer_empty.hpp
@@ -6,8 +6,6 @@
 template <class real_number, class partsize_t>
 class p2p_computer_empty{
 public:
-    constexpr int static size_data = 0;
-
     template <int size_particle_rhs>
     void init_result_array(real_number /*rhs*/[], const partsize_t /*nbParticles*/) const{
     }
@@ -16,15 +14,15 @@ public:
     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_data, int size_particle_rhs>
-    void compute_interaction(const real_number /*pos_part1*/[], const real_number /*data_part1*/[], real_number /*rhs_part1*/[],
-                             const real_number /*pos_part2*/[], const real_number /*data_part2*/[], real_number /*rhs_part2*/[],
-                             const real_number /*dist_pow2*/,
+    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{
     }
 
-    constexpr static bool isEmpty() {
-        return true;
+    constexpr static bool isEnable() {
+        return false;
     }
 };
 
diff --git a/bfps/cpp/particles/p2p_distr_mpi.hpp b/bfps/cpp/particles/p2p_distr_mpi.hpp
index ccc236d6..b009a57e 100644
--- a/bfps/cpp/particles/p2p_distr_mpi.hpp
+++ b/bfps/cpp/particles/p2p_distr_mpi.hpp
@@ -36,7 +36,6 @@ protected:
 
         std::unique_ptr<real_number[]> toRecvAndMerge;
         std::unique_ptr<real_number[]> toCompute;
-        std::unique_ptr<real_number[]> toData;
         std::unique_ptr<real_number[]> results;
     };
 
@@ -248,11 +247,10 @@ public:
         return (diff_x*diff_x) + (diff_y*diff_y) + (diff_z*diff_z);
     }
 
-    template <class computer_class, int size_particle_positions, int size_particle_data, int size_particle_rhs>
+    template <class computer_class, int size_particle_positions, int size_particle_rhs>
     void compute_distr(computer_class& in_computer,
                        const partsize_t current_my_nb_particles_per_partition[],
                        real_number particles_positions[],
-                       real_number particles_data[],
                        real_number particles_current_rhs[],
                        partsize_t inout_index_particles[]){
         TIMEZONE("compute_distr");
@@ -313,9 +311,6 @@ public:
                 permute_copy<real_number, size_particle_positions>(current_offset_particles_for_partition[idxPartition],
                                                                    current_my_nb_particles_per_partition[idxPartition],
                                                                    part_to_sort.data(), particles_positions, &buffer);
-                permute_copy<real_number, size_particle_data>(current_offset_particles_for_partition[idxPartition],
-                                                             current_my_nb_particles_per_partition[idxPartition],
-                                                             part_to_sort.data(), particles_data, &buffer);
                 permute_copy<real_number, size_particle_rhs>(current_offset_particles_for_partition[idxPartition],
                                                              current_my_nb_particles_per_partition[idxPartition],
                                                              part_to_sort.data(), particles_current_rhs, &buffer);
@@ -457,14 +452,6 @@ public:
                               descriptor.destProc, TAG_POSITION_PARTICLES,
                               current_com, &mpiRequests.back()));
 
-                    whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
-                    mpiRequests.emplace_back();
-                    assert(descriptor.nbParticlesToExchange*size_particle_data < std::numeric_limits<int>::max());
-                    AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_data[particles_offset_layers[my_nb_cell_levels-descriptor.nbLevelsToExchange]*size_particle_data]),
-                              int(descriptor.nbParticlesToExchange*size_particle_data), particles_utils::GetMpiType(real_number()),
-                              descriptor.destProc, TAG_POSITION_PARTICLES,
-                              current_com, &mpiRequests.back()));
-
                     assert(descriptor.toRecvAndMerge == nullptr);
                     descriptor.toRecvAndMerge.reset(new real_number[descriptor.nbParticlesToExchange*size_particle_rhs]);
                     whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
@@ -526,15 +513,6 @@ public:
                             AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions),
                                                 particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
                                                 current_com, &mpiRequests.back()));
-
-
-                            descriptor.toData.reset(new real_number[NbParticlesToReceive*size_particle_data]);
-                            whatNext.emplace_back(std::pair<Action,int>{COMPUTE_PARTICLES, releasedAction.second});
-                            mpiRequests.emplace_back();
-                            assert(NbParticlesToReceive*size_particle_data < std::numeric_limits<int>::max());
-                            AssertMpi(MPI_Irecv(descriptor.toData.get(), int(NbParticlesToReceive*size_particle_data),
-                                                particles_utils::GetMpiType(real_number()), destProc, TAG_POSITION_PARTICLES,
-                                                current_com, &mpiRequests.back()));
                         }
                     }
 
@@ -548,7 +526,6 @@ public:
                         const partsize_t NbParticlesToReceive = descriptor.nbParticlesToExchange;
 
                         assert(descriptor.toCompute != nullptr);
-                        assert(descriptor.toData != nullptr);
                         assert(descriptor.positionsReceived == true);
                         descriptor.results.reset(new real_number[NbParticlesToReceive*size_particle_rhs]);
                         in_computer.template init_result_array<size_particle_rhs>(descriptor.results.get(), NbParticlesToReceive);
@@ -589,14 +566,12 @@ public:
                                                                                                 particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
                                                                                                 shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
                                                 if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
-                                                    in_computer.template compute_interaction<size_particle_positions,size_particle_data, size_particle_rhs>(
+                                                    in_computer.template compute_interaction<size_particle_positions, size_particle_rhs>(
                                                                         &descriptor.toCompute[(idxPart+idx_p1)*size_particle_positions],
-                                                                        &descriptor.toData[(idxPart+idx_p1)*size_particle_data],
                                                                         &descriptor.results[(idxPart+idx_p1)*size_particle_rhs],
                                                                         &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
-                                                                        &particles_data[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_data],
                                                                         &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
-                                                                        dist_r2, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
+                                                                        dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
                                                 }
                                             }
                                         }
@@ -681,14 +656,12 @@ public:
                                                                             particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions + IDX_Z],
                                                                             0, 0, 0);
                             if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
-                                in_computer.template compute_interaction<size_particle_positions,size_particle_data,size_particle_rhs>(
+                                in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                     &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
-                                                    &particles_data[(intervals[idx_1].first+idx_p1)*size_particle_data],
                                                     &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
                                                     &particles_positions[(intervals[idx_1].first+idx_p2)*size_particle_positions],
-                                                    &particles_data[(intervals[idx_1].first+idx_p2)*size_particle_data],
                                                     &particles_current_rhs[(intervals[idx_1].first+idx_p2)*size_particle_rhs],
-                                                    dist_r2, 0, 0, 0);
+                                                    dist_r2, cutoff_radius_compute, 0, 0, 0);
                             }
                         }
                     }
@@ -705,14 +678,12 @@ public:
                                                                                 particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
                                                                                 0, 0, 0);
                                 if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
-                                    in_computer.template compute_interaction<size_particle_positions,size_particle_data,size_particle_rhs>(
+                                    in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                         &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
-                                                        &particles_data[(intervals[idx_1].first+idx_p1)*size_particle_data],
                                                         &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
                                                         &particles_positions[(intervals[idx_2].first+idx_p2)*size_particle_positions],
-                                                        &particles_data[(intervals[idx_2].first+idx_p2)*size_particle_data],
                                                         &particles_current_rhs[(intervals[idx_2].first+idx_p2)*size_particle_rhs],
-                                                        dist_r2, 0, 0, 0);
+                                                        dist_r2, cutoff_radius_compute, 0, 0, 0);
                                 }
                             }
                         }
@@ -741,14 +712,12 @@ public:
                                                                                         particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions + IDX_Z],
                                                                                         shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
                                         if(dist_r2 < cutoff_radius_compute*cutoff_radius_compute){
-                                            in_computer.template compute_interaction<size_particle_positions,size_particle_data,size_particle_rhs>(
+                                            in_computer.template compute_interaction<size_particle_positions,size_particle_rhs>(
                                                                 &particles_positions[(intervals[idx_1].first+idx_p1)*size_particle_positions],
-                                                                &particles_data[(intervals[idx_1].first+idx_p1)*size_particle_data],
                                                                 &particles_current_rhs[(intervals[idx_1].first+idx_p1)*size_particle_rhs],
                                                                 &particles_positions[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_positions],
-                                                                &particles_data[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_data],
                                                                 &particles_current_rhs[((*neighbors[idx_neighbor])[idx_2].first+idx_p2)*size_particle_rhs],
-                                                                dist_r2, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
+                                                                dist_r2, cutoff_radius_compute, shift[idx_neighbor][IDX_X], shift[idx_neighbor][IDX_Y], shift[idx_neighbor][IDX_Z]);
                                         }
                                     }
                                 }
diff --git a/bfps/cpp/particles/particles_adams_bashforth.hpp b/bfps/cpp/particles/particles_adams_bashforth.hpp
index 2fb61462..20c1ea50 100644
--- a/bfps/cpp/particles/particles_adams_bashforth.hpp
+++ b/bfps/cpp/particles/particles_adams_bashforth.hpp
@@ -7,11 +7,30 @@
 #include "scope_timer.hpp"
 #include "particles_utils.hpp"
 
-template <class partsize_t, class real_number, int size_particle_positions = 3, int size_particle_rhs = 3>
-class particles_adams_bashforth {
-    static_assert(size_particle_positions == size_particle_rhs,
-                  "Not having the same dimension for positions and rhs looks like a bug,"
-                  "otherwise comment this assertion.");
+template <class partsize_t, class real_number, int size_particle_positions, int size_particle_rhs>
+class particles_adams_bashforth;
+
+
+template <class partsize_t, class real_number>
+class particles_adams_bashforth<partsize_t,real_number,6,6>{
+    static const int size_particle_positions = 6;
+    static const int size_particle_rhs = 6;
+public:
+    static const int Max_steps = 6;
+
+    void move_particles(real_number*__restrict__ particles_positions,
+                        const partsize_t nb_particles,
+                        const std::unique_ptr<real_number[]> particles_rhs[],
+                        const int nb_rhs, const real_number dt) const{
+        // TODO
+    }
+};
+
+
+template <class partsize_t, class real_number>
+class particles_adams_bashforth<partsize_t,real_number,3,3>{
+    static const int size_particle_positions = 3;
+    static const int size_particle_rhs = 3;
 public:
     static const int Max_steps = 6;
 
diff --git a/bfps/cpp/particles/particles_distr_mpi.hpp b/bfps/cpp/particles/particles_distr_mpi.hpp
index f0c09fd9..79ff2e8b 100644
--- a/bfps/cpp/particles/particles_distr_mpi.hpp
+++ b/bfps/cpp/particles/particles_distr_mpi.hpp
@@ -35,9 +35,6 @@ protected:
         TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
         TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
 
-        TAG_LOW_UP_MOVED_PARTICLES_DATA,
-        TAG_UP_LOW_MOVED_PARTICLES_DATA,
-
         TAG_LOW_UP_MOVED_PARTICLES_RHS,
         TAG_LOW_UP_MOVED_PARTICLES_RHS_MAX = TAG_LOW_UP_MOVED_PARTICLES_RHS+MaxNbRhs,
 
@@ -506,14 +503,13 @@ public:
 
     ////////////////////////////////////////////////////////////////////////////
 
-    template <class computer_class, int size_particle_positions, int size_particle_data, int size_particle_rhs, int size_particle_index>
+    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,
-                      std::unique_ptr<partsize_t[]>* inout_index_particles,
-                      std::unique_ptr<real_number[]>* inout_data_particles){
+                      std::unique_ptr<partsize_t[]>* inout_index_particles){
         TIMEZONE("redistribute");
 
         // Some latest processes might not be involved
@@ -545,11 +541,6 @@ public:
                           (*inout_index_particles)[size_particle_index*idx2+idx_val]);
             }
 
-            for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
-                std::swap((*inout_data_particles)[size_particle_data*idx1+idx_val],
-                          (*inout_data_particles)[size_particle_data*idx2+idx_val]);
-            }
-
             for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                 for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
                     std::swap(inout_rhs_particles[idx_rhs][idx1*size_particle_rhs + idx_val],
@@ -576,11 +567,6 @@ public:
                           (*inout_index_particles)[size_particle_index*idx2+idx_val]);
             }
 
-            for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
-                std::swap((*inout_data_particles)[size_particle_data*idx1+idx_val],
-                          (*inout_data_particles)[size_particle_data*idx2+idx_val]);
-            }
-
             for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                 for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
                     std::swap(inout_rhs_particles[idx_rhs][idx1*size_particle_rhs + idx_val],
@@ -597,8 +583,6 @@ public:
         std::unique_ptr<real_number[]> newParticlesUp;
         std::unique_ptr<partsize_t[]> newParticlesLowIndexes;
         std::unique_ptr<partsize_t[]> newParticlesUpIndexes;
-        std::unique_ptr<real_number[]> newParticlesLowData;
-        std::unique_ptr<real_number[]> newParticlesUpData;
         std::vector<std::unique_ptr<real_number[]>> newParticlesLowRhs(in_nb_rhs);
         std::vector<std::unique_ptr<real_number[]>> newParticlesUpRhs(in_nb_rhs);
 
@@ -633,16 +617,6 @@ public:
                           (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
                           MPI_COMM_WORLD, &mpiRequests.back()));
 
-                if(size_particle_data){
-                    whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
-                    mpiRequests.emplace_back();
-                    assert(nbOutLower*size_particle_data < std::numeric_limits<int>::max());
-                    AssertMpi(MPI_Isend(&(*inout_data_particles)[0], int(nbOutLower*size_particle_data),
-                              particles_utils::GetMpiType(partsize_t()),
-                              (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_DATA,
-                              MPI_COMM_WORLD, &mpiRequests.back()));
-                }
-
                 for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                     whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
                     mpiRequests.emplace_back();
@@ -680,16 +654,6 @@ public:
                           particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
                           MPI_COMM_WORLD, &mpiRequests.back()));
 
-                if(size_particle_data){
-                    whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
-                    mpiRequests.emplace_back();
-                    assert(nbOutUpper*size_particle_data < std::numeric_limits<int>::max());
-                    AssertMpi(MPI_Isend(&(*inout_data_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_data],
-                              int(nbOutUpper*size_particle_data),
-                              particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_DATA,
-                              MPI_COMM_WORLD, &mpiRequests.back()));
-                }
-
                 for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                     whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
                     mpiRequests.emplace_back();
@@ -732,17 +696,6 @@ public:
                                   (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
                                   MPI_COMM_WORLD, &mpiRequests.back()));
 
-                        if(size_particle_data){
-                            newParticlesLowData.reset(new real_number[nbNewFromLow*size_particle_data]);
-                            whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
-                            mpiRequests.emplace_back();
-                            assert(nbNewFromLow*size_particle_data < std::numeric_limits<int>::max());
-                            AssertMpi(MPI_Irecv(&newParticlesLowData[0], int(nbNewFromLow*size_particle_data),
-                                      particles_utils::GetMpiType(real_number()),
-                                      (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_DATA,
-                                      MPI_COMM_WORLD, &mpiRequests.back()));
-                        }
-
                         for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                             newParticlesLowRhs[idx_rhs].reset(new real_number[nbNewFromLow*size_particle_rhs]);
                             whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -773,17 +726,6 @@ public:
                                   (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
                                   MPI_COMM_WORLD, &mpiRequests.back()));
 
-                        if(size_particle_data){
-                            newParticlesUpData.reset(new real_number[nbNewFromUp*size_particle_data]);
-                            whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
-                            mpiRequests.emplace_back();
-                            assert(nbNewFromUp*size_particle_data < std::numeric_limits<int>::max());
-                            AssertMpi(MPI_Irecv(&newParticlesUpData[0], int(nbNewFromUp*size_particle_data),
-                                      particles_utils::GetMpiType(real_number()),
-                                      (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_DATA,
-                                      MPI_COMM_WORLD, &mpiRequests.back()));
-                        }
-
                         for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                             newParticlesUpRhs[idx_rhs].reset(new real_number[nbNewFromUp*size_particle_rhs]);
                             whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
@@ -814,7 +756,6 @@ public:
 
             std::unique_ptr<real_number[]> newArray(new real_number[myTotalNewNbParticles*size_particle_positions]);
             std::unique_ptr<partsize_t[]> newArrayIndexes(new partsize_t[myTotalNewNbParticles*size_particle_index]);
-            std::unique_ptr<real_number[]> newArrayData(new real_number[myTotalNewNbParticles*size_particle_data]);
             std::vector<std::unique_ptr<real_number[]>> newArrayRhs(in_nb_rhs);
             for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                 newArrayRhs[idx_rhs].reset(new real_number[myTotalNewNbParticles*size_particle_rhs]);
@@ -825,7 +766,6 @@ public:
                 const particles_utils::fixed_copy fcp(0, 0, nbNewFromLow);
                 fcp.copy(newArray, newParticlesLow, size_particle_positions);
                 fcp.copy(newArrayIndexes, newParticlesLowIndexes, size_particle_index);
-                fcp.copy(newArrayData, newParticlesLowData, size_particle_data);
                 for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                     fcp.copy(newArrayRhs[idx_rhs], newParticlesLowRhs[idx_rhs], size_particle_rhs);
                 }
@@ -836,7 +776,6 @@ public:
                 const particles_utils::fixed_copy fcp(nbNewFromLow, nbOutLower, nbOldParticlesInside);
                 fcp.copy(newArray, (*inout_positions_particles), size_particle_positions);
                 fcp.copy(newArrayIndexes, (*inout_index_particles), size_particle_index);
-                fcp.copy(newArrayData, (*inout_data_particles), size_particle_data);
                 for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                     fcp.copy(newArrayRhs[idx_rhs], inout_rhs_particles[idx_rhs], size_particle_rhs);
                 }
@@ -847,7 +786,6 @@ public:
                 const particles_utils::fixed_copy fcp(nbNewFromLow+nbOldParticlesInside, 0, nbNewFromUp);
                 fcp.copy(newArray, newParticlesUp, size_particle_positions);
                 fcp.copy(newArrayIndexes, newParticlesUpIndexes, size_particle_index);
-                fcp.copy(newArrayData, newParticlesUpData, size_particle_data);
                 for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                     fcp.copy(newArrayRhs[idx_rhs], newParticlesUpRhs[idx_rhs], size_particle_rhs);
                 }
@@ -855,7 +793,6 @@ public:
 
             (*inout_positions_particles) = std::move(newArray);
             (*inout_index_particles) = std::move(newArrayIndexes);
-            (*inout_data_particles) = std::move(newArrayData);
             for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                 inout_rhs_particles[idx_rhs] = std::move(newArrayRhs[idx_rhs]);
             }
@@ -880,11 +817,6 @@ public:
                               (*inout_index_particles)[size_particle_index*idx2 + idx_val]);
                 }
 
-                for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
-                    std::swap((*inout_data_particles)[size_particle_data*idx1 + idx_val],
-                              (*inout_data_particles)[size_particle_data*idx2 + idx_val]);
-                }
-
                 for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
                     for(int idx_val = 0 ; idx_val < size_particle_rhs ; ++idx_val){
                         std::swap(inout_rhs_particles[idx_rhs][idx1*size_particle_rhs + idx_val],
diff --git a/bfps/cpp/particles/particles_field_computer.hpp b/bfps/cpp/particles/particles_field_computer.hpp
index f68f2fc0..70da4757 100644
--- a/bfps/cpp/particles/particles_field_computer.hpp
+++ b/bfps/cpp/particles/particles_field_computer.hpp
@@ -67,8 +67,9 @@ public:
                                    const real_number particles_positions[],
                                    real_number particles_current_rhs[],
                                    const partsize_t nb_particles) const {
+        static_assert(field_class::nb_components() <= size_particle_rhs, "Cannot store all the component in the given array");
         TIMEZONE("particles_field_computer::apply_computation");
-        //DEBUG_MSG("just entered particles_field_computer::apply_computation\n");
+
         for(partsize_t idxPart = 0 ; idxPart < nb_particles ; ++idxPart){
             const real_number reltv_x = get_norm_pos_in_cell(particles_positions[idxPart*3+IDX_X], IDX_X);
             const real_number reltv_y = get_norm_pos_in_cell(particles_positions[idxPart*3+IDX_Y], IDX_Y);
@@ -146,7 +147,8 @@ public:
                             const ptrdiff_t tindex = field.get_rindex_from_global(idx_x_pbc, idx_y_pbc, idx_z_pbc);
 
                             // getValue does not necessary return real_number
-                            for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){
+                            // size_particle_rhs is just for the leading dimension of the array
+                            for(int idx_rhs_val = 0 ; idx_rhs_val < field_class::nb_components() ; ++idx_rhs_val){
                                 particles_current_rhs[idxPart*size_particle_rhs+idx_rhs_val] += real_number(field.rval(tindex,idx_rhs_val))*coef;
                             }
                         }
diff --git a/bfps/cpp/particles/particles_inner_computer.hpp b/bfps/cpp/particles/particles_inner_computer.hpp
new file mode 100644
index 00000000..4d0a6678
--- /dev/null
+++ b/bfps/cpp/particles/particles_inner_computer.hpp
@@ -0,0 +1,65 @@
+#ifndef PARTICLES_INNER_COMPUTER_HPP
+#define PARTICLES_INNER_COMPUTER_HPP
+
+#include <cstring>
+#include <cassert>
+
+template <class real_number, class partsize_t>
+class particles_inner_computer{
+    bool isActive;
+    const real_number v0;
+
+public:
+    explicit particles_inner_computer(const real_number inV0) :  isActive(true), v0(inV0){
+    }
+
+    template <int size_particle_positions, int size_particle_rhs>
+    void compute_interaction(const partsize_t nb_particles, real_number pos_part[], real_number rhs_part[]) const{
+        static_assert(size_particle_positions == 6, "This kernel works only with 6 values for one particle's position");
+        static_assert(size_particle_rhs == 6, "This kernel works only with 6 values per particle's rhs");
+
+        #pragma omp parallel for
+        for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+            // Add attr × V0 to the field interpolation
+            rhs_part[idx_part*size_particle_rhs + IDX_X] += pos_part[idx_part*size_particle_positions + 3+IDX_X]*v0;
+            rhs_part[idx_part*size_particle_rhs + IDX_Y] += pos_part[idx_part*size_particle_positions + 3+IDX_Y]*v0;
+            rhs_part[idx_part*size_particle_rhs + IDX_Z] += pos_part[idx_part*size_particle_positions + 3+IDX_Z]*v0;
+
+            real_number alpha[3]= {0}; // TODO compute aplha
+
+            rhs_part[idx_part*size_particle_rhs + 3+IDX_X] += alpha[IDX_X];
+            rhs_part[idx_part*size_particle_rhs + 3+IDX_Y] += alpha[IDX_Y];
+            rhs_part[idx_part*size_particle_rhs + 3+IDX_Z] += alpha[IDX_Z];
+        }
+    }
+
+    template <int size_particle_positions, int size_particle_rhs, int size_particle_rhs_extra>
+    void compute_interaction_with_extra(const partsize_t nb_particles, real_number pos_part[], real_number rhs_part[],
+                                        const real_number rhs_part_extra[]) const{
+        static_assert(size_particle_rhs_extra == 3, "This kernel works only with 3 values for one particle's rhs extra");
+
+        compute_interaction<size_particle_positions, size_particle_rhs>(nb_particles, pos_part, rhs_part);
+
+        #pragma omp parallel for
+        for(partsize_t idx_part = 0 ; idx_part < nb_particles ; ++idx_part){
+            // Cross product vorticity/orientation
+            rhs_part[idx_part*size_particle_rhs + 3+IDX_X] += rhs_part_extra[idx_part*size_particle_rhs_extra + IDX_Y]*pos_part[idx_part*size_particle_positions + 3+IDX_Z]
+                                                             - rhs_part_extra[idx_part*size_particle_rhs_extra + IDX_Z]*pos_part[idx_part*size_particle_positions + 3+IDX_Y];
+            rhs_part[idx_part*size_particle_rhs + 3+IDX_Y] += rhs_part_extra[idx_part*size_particle_rhs_extra + IDX_Z]*pos_part[idx_part*size_particle_positions + 3+IDX_X]
+                                                             - rhs_part_extra[idx_part*size_particle_rhs_extra + IDX_X]*pos_part[idx_part*size_particle_positions + 3+IDX_Z];
+            rhs_part[idx_part*size_particle_rhs + 3+IDX_Z] += rhs_part_extra[idx_part*size_particle_rhs_extra + IDX_X]*pos_part[idx_part*size_particle_positions + 3+IDX_Y]
+                                                             - rhs_part_extra[idx_part*size_particle_rhs_extra + IDX_Y]*pos_part[idx_part*size_particle_positions + 3+IDX_X];
+        }
+    }
+
+
+    bool isEnable() const {
+        return isActive;
+    }
+
+    void setEnable(const bool inIsActive) {
+        isActive = inIsActive;
+    }
+};
+
+#endif
diff --git a/bfps/cpp/particles/particles_inner_computer_empty.hpp b/bfps/cpp/particles/particles_inner_computer_empty.hpp
new file mode 100644
index 00000000..1bd3b1ec
--- /dev/null
+++ b/bfps/cpp/particles/particles_inner_computer_empty.hpp
@@ -0,0 +1,24 @@
+#ifndef PARTICLES_INNER_COMPUTER_EMPTY_HPP
+#define PARTICLES_INNER_COMPUTER_EMPTY_HPP
+
+#include <cstring>
+#include <cassert>
+
+template <class real_number, class partsize_t>
+class particles_inner_computer_empty{
+public:
+    template <int size_particle_positions, int size_particle_rhs>
+    void compute_interaction(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[]) const{
+    }
+
+    template <int size_particle_positions, int size_particle_rhs, int size_particle_rhs_extra>
+    void compute_interaction_with_extra(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[],
+                             const real_number /*rhs_part_extra*/[]) const{
+    }
+
+    constexpr static bool isEnable() {
+        return false;
+    }
+};
+
+#endif
diff --git a/bfps/cpp/particles/particles_system.hpp b/bfps/cpp/particles/particles_system.hpp
index 12bf6c29..2b56cae8 100644
--- a/bfps/cpp/particles/particles_system.hpp
+++ b/bfps/cpp/particles/particles_system.hpp
@@ -15,7 +15,7 @@
 #include "p2p_distr_mpi.hpp"
 
 template <class partsize_t, class real_number, class field_rnumber, class field_class, class interpolator_class, int interp_neighbours,
-          int size_particle_rhs, class p2p_computer_class, int size_particle_data>
+          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> {
     MPI_Comm mpi_com;
 
@@ -26,7 +26,7 @@ class particles_system : public abstract_particles_system<partsize_t, real_numbe
 
     particles_distr_mpi<partsize_t, real_number> particles_distr;
 
-    particles_adams_bashforth<partsize_t, real_number, 3, size_particle_rhs> positions_updater;
+    particles_adams_bashforth<partsize_t, real_number, size_particle_positions, size_particle_rhs> positions_updater;
 
     using computer_class = particles_field_computer<partsize_t, real_number, interpolator_class, interp_neighbours>;
     computer_class computer;
@@ -51,7 +51,7 @@ class particles_system : public abstract_particles_system<partsize_t, real_numbe
 
     p2p_distr_mpi<partsize_t, real_number> distr_p2p;
     p2p_computer_class computer_p2p;
-    std::unique_ptr<real_number[]> my_particles_data;
+    particles_inner_computer_class computer_particules_inner;
 
 public:
     particles_system(const std::array<size_t,3>& field_grid_dim, const std::array<real_number,3>& in_spatial_box_width,
@@ -63,8 +63,10 @@ public:
                      const field_class& in_field,
                      MPI_Comm in_mpi_com,
                      const partsize_t in_total_nb_particles,
-                     const int in_current_iteration = 1,
-                     const real_number in_cutoff = 1.)
+                     const real_number in_cutoff,
+                     p2p_computer_class in_computer_p2p,
+                     particles_inner_computer_class in_computer_particules_inner,
+                     const int in_current_iteration = 1)
         : mpi_com(in_mpi_com),
           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),
@@ -77,7 +79,8 @@ public:
           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), total_nb_particles(in_total_nb_particles), step_idx(in_current_iteration),
-          distr_p2p(in_mpi_com, current_partition_interval,field_grid_dim, spatial_box_width, in_spatial_box_offset, in_cutoff){
+          distr_p2p(in_mpi_com, current_partition_interval,field_grid_dim, spatial_box_width, in_spatial_box_offset, in_cutoff),
+          computer_p2p(std::move(in_computer_p2p)), computer_particules_inner(std::move(in_computer_particules_inner)){
 
         current_my_nb_particles_per_partition.reset(new partsize_t[partition_interval_size]);
         current_offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
@@ -93,16 +96,14 @@ public:
         my_particles_positions_indexes = particles_input.getMyParticlesIndexes();
         my_particles_rhs = particles_input.getMyRhs();
         my_nb_particles = particles_input.getLocalNbParticles();
-        // TODO P2P get it from loader
-        my_particles_data.reset(new real_number[my_nb_particles*size_particle_data]());
 
         for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me
-            const int partition_level = computer.pbc_field_layer(my_particles_positions[idx_part*3+IDX_Z], IDX_Z);
+            const int partition_level = computer.pbc_field_layer(my_particles_positions[idx_part*size_particle_positions+IDX_Z], IDX_Z);
             assert(partition_level >= current_partition_interval.first);
             assert(partition_level < current_partition_interval.second);
         }
 
-        particles_utils::partition_extra_z<partsize_t, 3>(&my_particles_positions[0], my_nb_particles, partition_interval_size,
+        particles_utils::partition_extra_z<partsize_t, size_particle_positions>(&my_particles_positions[0], my_nb_particles, partition_interval_size,
                                               current_my_nb_particles_per_partition.get(), current_offset_particles_for_partition.get(),
         [&](const real_number& z_pos){
             const int partition_level = computer.pbc_field_layer(z_pos, IDX_Z);
@@ -117,10 +118,6 @@ public:
                               my_particles_rhs[idx_rhs][idx2*size_particle_rhs + idx_val]);
                 }
             }
-            for(int idx_val = 0 ; idx_val < size_particle_data ; ++idx_val){
-                std::swap(my_particles_data[idx1*size_particle_data + idx_val],
-                          my_particles_data[idx2*size_particle_data + idx_val]);
-            }
         });
 
         {// TODO remove
@@ -128,16 +125,15 @@ 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(computer.pbc_field_layer(my_particles_positions[idx*3+IDX_Z], IDX_Z)-current_partition_interval.first == idxPartition);
+                    assert(computer.pbc_field_layer(my_particles_positions[idx*size_particle_positions+IDX_Z], IDX_Z)-current_partition_interval.first == idxPartition);
                 }
             }
         }
     }
 
-
     void compute() final {
         TIMEZONE("particles_system::compute");
-        particles_distr.template compute_distr<computer_class, field_class, 3, size_particle_rhs>(
+        particles_distr.template compute_distr<computer_class, field_class, size_particle_positions, size_particle_rhs>(
                                computer, default_field,
                                current_my_nb_particles_per_partition.get(),
                                my_particles_positions.get(),
@@ -147,20 +143,37 @@ public:
 
     void compute_p2p() final {
         // TODO P2P
-        if(p2p_computer_class::isEmpty() == false){
+        if(computer_p2p.isEnable() == true){
             TIMEZONE("particles_system::compute_p2p");
-            distr_p2p.template compute_distr<p2p_computer_class, 3, size_particle_data, size_particle_rhs>(
+            distr_p2p.template compute_distr<p2p_computer_class, size_particle_positions, size_particle_rhs>(
                             computer_p2p, current_my_nb_particles_per_partition.get(),
-                            my_particles_positions.get(), my_particles_data.get(), my_particles_rhs.front().get(),
+                            my_particles_positions.get(), my_particles_rhs.front().get(),
                             my_particles_positions_indexes.get());
         }
     }
 
+    void compute_particles_inner() final {
+        if(computer_particules_inner.isEnable() == true){
+            TIMEZONE("particles_system::compute_particles_inner");
+            computer_particules_inner.template compute_interaction<size_particle_positions, size_particle_rhs>(
+                            my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get());
+        }
+    }
+
+    void compute_particles_inner(const real_number particle_extra_rhs[]) final {
+        if(computer_particules_inner.isEnable() == true){
+            TIMEZONE("particles_system::compute_particles_inner");
+            computer_particules_inner.template compute_interaction_with_extra<size_particle_positions, size_particle_rhs, 3>(
+                            my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get(),
+                            particle_extra_rhs);
+        }
+    }
+
     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>(
+        particles_distr.template compute_distr<computer_class, sample_field_class, size_particle_positions, sample_size_particle_rhs>(
                                computer, sample_field,
                                current_my_nb_particles_per_partition.get(),
                                my_particles_positions.get(),
@@ -204,14 +217,13 @@ public:
 
     void redistribute() final {
         TIMEZONE("particles_system::redistribute");
-        particles_distr.template redistribute<computer_class, 3, size_particle_data, size_particle_rhs, 1>(
+        particles_distr.template redistribute<computer_class, size_particle_positions, 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()),
-                              &my_particles_positions_indexes,
-                              &my_particles_data);
+                              &my_particles_positions_indexes);
     }
 
     void inc_step_idx() final {
@@ -237,6 +249,19 @@ public:
         TIMEZONE("particles_system::completeLoop");
         compute();
         compute_p2p();
+        compute_particles_inner();
+        move(dt);
+        redistribute();
+        inc_step_idx();
+        shift_rhs_vectors();
+    }
+
+    void completeLoopWithVorticity(const real_number dt,
+                                   const real_number particle_extra_rhs[]) final {
+        TIMEZONE("particles_system::completeLoop");
+        compute();
+        compute_p2p();
+        compute_particles_inner(particle_extra_rhs);
         move(dt);
         redistribute();
         inc_step_idx();
@@ -269,9 +294,9 @@ public:
 
     void checkNan() const { // TODO remove
         for(partsize_t idx_part = 0 ; idx_part < my_nb_particles ; ++idx_part){ // TODO remove me
-            assert(std::isnan(my_particles_positions[idx_part*3+IDX_X]) == false);
-            assert(std::isnan(my_particles_positions[idx_part*3+IDX_Y]) == false);
-            assert(std::isnan(my_particles_positions[idx_part*3+IDX_Z]) == false);
+            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDX_X]) == false);
+            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDX_Y]) == false);
+            assert(std::isnan(my_particles_positions[idx_part*size_particle_positions+IDX_Z]) == false);
 
             for(int idx_rhs = 0 ; idx_rhs < my_particles_rhs.size() ; ++idx_rhs){
                 for(int idx_rhs_val = 0 ; idx_rhs_val < size_particle_rhs ; ++idx_rhs_val){
diff --git a/bfps/cpp/particles/particles_system_builder.hpp b/bfps/cpp/particles/particles_system_builder.hpp
index f9ce512d..ba078000 100644
--- a/bfps/cpp/particles/particles_system_builder.hpp
+++ b/bfps/cpp/particles/particles_system_builder.hpp
@@ -9,6 +9,7 @@
 #include "particles_input_hdf5.hpp"
 #include "particles_generic_interp.hpp"
 #include "p2p_computer_empty.hpp"
+#include "particles_inner_computer_empty.hpp"
 
 #include "field.hpp"
 #include "kspace.hpp"
@@ -110,7 +111,8 @@ inline RetType evaluate(IterType1 value1, IterType2 value2, Args... args){
 ///
 //////////////////////////////////////////////////////////////////////////////
 
-template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber, class p2p_computer_class>
+template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class particles_rnumber, class p2p_computer_class,
+          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(
@@ -121,7 +123,10 @@ struct particles_system_build_container {
              const std::string& fname_input, // particles input filename
             const std::string& inDatanameState, const std::string& inDatanameRhs, // input dataset names
              MPI_Comm mpi_comm,
-            const int in_current_iteration){
+            const int in_current_iteration,
+            p2p_computer_class p2p_computer,
+            particles_inner_computer_class inner_computer,
+            const particles_rnumber cutoff = std::numeric_limits<particles_rnumber>::max()){
 
         // The size of the field grid (global size) all_size seems
         std::array<size_t,3> field_grid_dim;
@@ -197,8 +202,10 @@ struct particles_system_build_container {
         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),
-                                                       p2p_computer_class, p2p_computer_class::size_data>;
+                                                       interpolation_size,
+                                                       size_particle_positions, size_particle_rhs,
+                                                       p2p_computer_class,
+                                                       particles_inner_computer_class>;
         particles_system_type* part_sys = new particles_system_type(field_grid_dim,
                                                spatial_box_width,
                                                spatial_box_offset,
@@ -210,6 +217,9 @@ struct particles_system_build_container {
                                                (*fs_field),
                                                mpi_comm,
                                                nparticles,
+                                               cutoff,
+                                               p2p_computer,
+                                               inner_computer,
                                                in_current_iteration);
 
         // TODO P2P load particle data too
@@ -254,13 +264,19 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
     return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
                        int, 1, 11, 1, // interpolation_size
                        int, 0, 3, 1, // spline_mode
-                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber, p2p_computer_empty<particles_rnumber,partsize_t>>>(
+                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
+                                                        p2p_computer_empty<particles_rnumber,partsize_t>,
+                                                        particles_inner_computer_empty<particles_rnumber,partsize_t>,
+                                                        3,3>>(
                            interpolation_size, // template iterator 1
                            spline_mode, // template iterator 2
-                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration);
+                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration,
+                           p2p_computer_empty<particles_rnumber,partsize_t>(), particles_inner_computer_empty<particles_rnumber,partsize_t>());
 }
 
-template <class partsize_t, class field_rnumber, field_backend be, field_components fc, class p2p_computer_class, class particles_rnumber = double>
+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(
         const field<field_rnumber, be, fc>* fs_field, // (field object)
         const kspace<be, SMOOTH>* fs_kk, // (kspace object, contains dkx, dky, dkz)
@@ -271,14 +287,21 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
         const int interpolation_size,
         const int spline_mode,
         MPI_Comm mpi_comm,
-        const int in_current_iteration){
+        const int in_current_iteration,
+        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>>,
                        int, 1, 11, 1, // interpolation_size
                        int, 0, 3, 1, // spline_mode
-                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,p2p_computer_class>>(
+                       particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
+                                                        p2p_computer_class,
+                                                        particles_inner_computer_class,
+                                                        6,6>>(
                            interpolation_size, // template iterator 1
                            spline_mode, // template iterator 2
-                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration);
+                           fs_field,fs_kk, nsteps, nparticles, fname_input, inDatanameState, inDatanameRhs, mpi_comm, in_current_iteration,
+                           std::move(p2p_computer), std::move(inner_computer), cutoff);
 }
 
 
-- 
GitLab