[wip] new class hierarchy for particles

preliminary. in particular solvers cannot be a simple function, they
should handle checkpointing --- especially if Adams-Bashforth particles
are to be reimplemented with the new classes.

CMakeLists.txt

@@ -297,6 +297,8 @@ set(cpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/spline_n9.cpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/spline_n10.cpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/Lagrange_polys.cpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/field_tinterpolator.cpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/particle_set.cpp
     ${PROJECT_SOURCE_DIR}/cpp/scope_timer.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp
@@ -348,6 +350,11 @@ set(hpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/spline_n9.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/spline_n10.hpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/Lagrange_polys.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/abstract_particle_set.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/particle_set.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/interpolation/field_tinterpolator.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particle_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
     ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.hpp

cpp/particles/abstract_particle_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 ABSTRACT_PARTICLE_RHS_HPP
+#define ABSTRACT_PARTICLE_RHS_HPP
+
+#include "particles/interpolation/abstract_particle_set.hpp"
+
+// children of this class must provide a method to compute the rhs for a given particle set
+class abstract_particle_rhs
+{
+    public:
+        // destructor
+        virtual ~abstract_particle_rhs(){}
+
+        // important bit
+        virtual int operator()(
+                double t,
+                const abstract_particle_set &pset,
+                particle_rnumber *result) const = 0;
+};
+
+#endif//ABSTRACT_PARTICLE_RHS_HPP
+

cpp/particles/interpolation/abstract_particle_set.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 ABSTRACT_PARTICLE_SET_HPP
+#define ABSTRACT_PARTICLE_SET_HPP
+
+#include <array>
+#include <cassert>
+#include <vector>
+#include "field.hpp"
+
+// the point here is to have a class that stores particle state data,
+// knows how to distribute this data among processes,
+// and knows how to compute interpolations.
+// the processes of MPI distribution and computing interpolations can't really be separated
+// so it makes sense to merge the functionalities.
+class abstract_particle_set
+{
+    protected:
+        using partsize_t = long long int;
+        using particle_rnumber = double;
+
+    public:
+        // destructor
+        virtual ~abstract_particle_set(){}
+
+        // extract particle information
+        virtual const particle_rnumber* getParticleState() const = 0;
+        virtual const partsize_t* getParticleIndex() const = 0;
+
+        virtual std::unique_ptr<particle_rnumber[]> extractFromParticleState(
+                const int firstComponent,
+                const int lastComponent) const = 0;
+
+        virtual partsize_t getLocalNumberOfParticles() const = 0;
+        virtual partsize_t getTotalNumberOfParticles() const = 0;
+        virtual int getStateSize() const = 0;
+
+        // reorganize particles within MPI domain
+        virtual int redistribute() = 0;
+
+        // sample field values at particle location
+        virtual int sample(
+                const field<float,
+                            FFTW,
+                            ONE>         &field_to_sample,
+                particle_rnumber *result) const = 0;
+        virtual int sample(
+                const field<float,
+                            FFTW,
+                            THREE>       &field_to_sample,
+                particle_rnumber *result) const = 0;
+        virtual int sample(
+                const field<float,
+                            FFTW,
+                            THREExTHREE> &field_to_sample,
+                particle_rnumber *result) const = 0;
+        virtual int sample(
+                const field<double,
+                            FFTW,
+                            ONE>         &field_to_sample,
+                particle_rnumber *result) const = 0;
+        virtual int sample(
+                const field<double,
+                            FFTW,
+                            THREE>       &field_to_sample,
+                particle_rnumber *result) const = 0;
+        virtual int sample(
+                const field<double,
+                            FFTW,
+                            THREExTHREE> &field_to_sample,
+                particle_rnumber *result) const = 0;
+};
+
+#endif//ABSTRACT_PARTICLE_SET_HPP
+

cpp/particles/interpolation/field_tinterpolator.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/interpolation/field_tinterpolator.hpp"
+
+template class field_tinterpolator<float,  FFTW, ONE,         NONE>;
+template class field_tinterpolator<float,  FFTW, THREE,       NONE>;
+template class field_tinterpolator<float,  FFTW, THREExTHREE, NONE>;
+
+template class field_tinterpolator<double, FFTW, ONE,         NONE>;
+template class field_tinterpolator<double, FFTW, THREE,       NONE>;
+template class field_tinterpolator<double, FFTW, THREExTHREE, NONE>;
+
+template class field_tinterpolator<float,  FFTW, ONE,         LINEAR>;
+template class field_tinterpolator<float,  FFTW, THREE,       LINEAR>;
+template class field_tinterpolator<float,  FFTW, THREExTHREE, LINEAR>;
+
+template class field_tinterpolator<double, FFTW, ONE,         LINEAR>;
+template class field_tinterpolator<double, FFTW, THREE,       LINEAR>;
+template class field_tinterpolator<double, FFTW, THREExTHREE, LINEAR>;
+

cpp/particles/interpolation/field_tinterpolator.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 FIELD_TINTERPOLATOR_HPP
+#define FIELD_TINTERPOLATOR_HPP
+
+#include <array>
+#include <cassert>
+#include <vector>
+#include "field.hpp"
+#include "particles/particles_distr_mpi.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/interpolation/abstract_particle_set.hpp"
+
+enum temporal_interpolation_type {NONE, LINEAR, CUBIC};
+
+template <typename rnumber,
+          field_backend be,
+          field_components fc,
+          temporal_interpolation_type tt>
+class field_tinterpolator
+{
+    private:
+        using particle_rnumber = double;
+        using partsize_t = long long int;
+
+        std::array<field<rnumber, be, fc>*, 4> field_list;  // list of fields for temporal interpolation
+
+    public:
+        field_tinterpolator(){
+            this->field_list[0] = NULL;
+            this->field_list[1] = NULL;
+            this->field_list[2] = NULL;
+            this->field_list[3] = NULL;
+        }
+
+        ~field_tinterpolator(){}
+
+        int set_field(
+                field<rnumber, be, fc> *field_src = NULL,
+                const int tindex = 0)
+        {
+            switch(tt)
+            {
+                case NONE:
+                    assert(tindex == 0);
+                    this->field_list[0] = field_src;
+                    break;
+                case LINEAR:
+                    assert(tindex == 0 || tindex == 1);
+                    this->field_list[tindex] = field_src;
+                    break;
+                case CUBIC:
+                    assert(false);
+                    break;
+            }
+            return EXIT_SUCCESS;
+        }
+
+        // t is a fraction between 0 and 1.
+        int operator()(
+                double t,
+                const abstract_particle_set *pset,
+                particle_rnumber *result)
+        {
+            switch(tt)
+            {
+                case NONE:
+                    {
+                        pset->sample(*(this->field_list[0]), result);
+                        break;
+                    }
+                case LINEAR:
+                    {
+                        particle_rnumber *result0, *result1;
+                        result0 = new particle_rnumber[pset->getLocalNumberOfParticles()*ncomp(fc)];
+                        result1 = new particle_rnumber[pset->getLocalNumberOfParticles()*ncomp(fc)];
+                        pset->sample(*(this->field_list[0]), result0);
+                        pset->sample(*(this->field_list[1]), result1);
+                        for (partsize_t idx_part = 0; idx_part < pset->getLocalNumberOfParticles(); idx_part++)
+                        {
+                            for (unsigned int cc = 0; cc < ncomp(fc); cc++)
+                            {
+                                result[idx_part*ncomp(fc) + cc] = (
+                                        (1-t)*result0[idx_part*ncomp(fc) + cc] +
+                                           t *result1[idx_part*ncomp(fc) + cc]);
+                            }
+                        }
+                        delete[] result0;
+                        delete[] result1;
+                        break;
+                    }
+                case CUBIC:
+                    assert(false);
+                    break;
+            }
+            return EXIT_SUCCESS;
+        }
+};
+
+#endif//FIELD_TINTERPOLATOR_HPP
+

cpp/particles/interpolation/particle_set.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/interpolation/particle_set.hpp"
+
+template class particle_set<3, 1, 0>;
+template class particle_set<3, 1, 1>;
+template class particle_set<3, 2, 0>;
+template class particle_set<3, 2, 1>;
+template class particle_set<3, 2, 2>;
+

cpp/particles/interpolation/particle_set.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 PARTICLE_SET_HPP
+#define PARTICLE_SET_HPP
+
+#include <array>
+#include <cassert>
+#include <vector>
+#include "field.hpp"
+#include "particles/particles_distr_mpi.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/interpolation/particles_field_computer.hpp"
+#include "particles/interpolation/particles_generic_interp.hpp"
+#include "particles/interpolation/abstract_particle_set.hpp"
+
+template <int state_size,
+          int neighbours,
+          int smoothness>
+class particle_set: public abstract_particle_set
+{
+    private:
+        // custom names for types
+        using interpolator_class = particles_field_computer<partsize_t,
+                                                            particle_rnumber,
+                                                            particles_generic_interp<particle_rnumber, neighbours, smoothness>,
+                                                            neighbours>;
+        using distributor_class = particles_distr_mpi<partsize_t, particle_rnumber>;
+
+        // MPI
+        MPI_Comm mpi_comm;
+
+        // data
+        const std::pair<int,int> current_partition_interval;
+        const int partition_interval_size;
+
+        std::unique_ptr<partsize_t[]> number_particles_per_partition;
+        std::unique_ptr<partsize_t[]> offset_particles_for_partition;
+
+        std::unique_ptr<particle_rnumber[]> local_state;
+        std::unique_ptr<partsize_t[]>       local_index;
+        partsize_t                          local_number_of_particles;
+        partsize_t                          total_number_of_particles;
+
+        // for I/O it helps to have custom array shape, we store it here
+        std::vector<hsize_t> particle_file_layout;
+
+        // raw objects for distributing data in MPI, and for interpolating fields
+        distributor_class   pdistributor;
+        particles_generic_interp<particle_rnumber, neighbours, smoothness> interpolation_formula;
+        interpolator_class  pinterpolator;
+
+    public:
+        particle_set(
+                const MPI_Comm in_mpi_comm,
+                const std::array<particle_rnumber, 3> &in_spatial_box_width,
+                const std::array<particle_rnumber, 3> &in_spatial_box_offset,
+                const std::array<particle_rnumber, 3> &in_spatial_partition_width,
+                const std::array<size_t, 3> &field_grid_dim,
+                const std::array<size_t, 3> &in_local_field_offset,
+                const std::array<size_t, 3> &in_local_field_dims):
+            mpi_comm(in_mpi_comm),
+            current_partition_interval(
+                    {in_local_field_offset[IDXC_Z],
+                     in_local_field_offset[IDXC_Z] + in_local_field_dims[IDXC_Z]}),
+            partition_interval_size(current_partition_interval.second - current_partition_interval.first),
+            pdistributor(
+                    in_mpi_comm,
+                    current_partition_interval,
+                    field_grid_dim),
+            interpolation_formula(),
+            pinterpolator(
+                    field_grid_dim,
+                    current_partition_interval,
+                    interpolation_formula,
+                    in_spatial_box_width,
+                    in_spatial_box_offset,
+                    in_spatial_partition_width)
+        {
+            this->number_particles_per_partition.reset(new partsize_t[partition_interval_size]);
+            this->offset_particles_for_partition.reset(new partsize_t[partition_interval_size+1]);
+        }
+
+        ~particle_set(){}
+
+        const particle_rnumber* getParticleState() const
+        {
+            return this->local_state.get();
+        }
+
+        const partsize_t* getParticleIndex() const
+        {
+            return this->local_index.get();
+        }
+
+        partsize_t getLocalNumberOfParticles() const
+        {
+            return this->local_number_of_particles;
+        }
+
+        partsize_t getTotalNumberOfParticles() const
+        {
+            return this->total_number_of_particles;
+        }
+
+        int setParticleFileLayout(std::vector<hsize_t> input_layout)
+        {
+            this->particle_file_layout.resize(input_layout.size());
+            for (unsigned int i=0; i<this->particle_file_layout.size(); i++)
+                this->particle_file_layout[i] = input_layout[i];
+            return EXIT_SUCCESS;
+        }
+
+        std::vector<hsize_t> getParticleFileLayout(void)
+        {
+            return std::vector<hsize_t>(this->particle_file_layout);
+        }
+
+        std::unique_ptr<particle_rnumber[]> extractParticlesState(
+                const int firstState,
+                const int lastState) const
+        {
+            const int numberOfStates = std::max(0,(std::min(lastState, state_size)-firstState));
+
+            std::unique_ptr<particle_rnumber[]> stateExtract(new particle_rnumber[local_number_of_particles*numberOfStates]);
+
+            for(partsize_t idx_part = 0 ; idx_part < this->local_number_of_particles ; ++idx_part){
+                for(int idxState = 0 ; idxState < numberOfStates ; ++idxState){
+                    stateExtract[idx_part*numberOfStates + idxState] = this->local_state[idx_part*state_size + idxState+firstState];
+                }
+            }
+
+            return stateExtract;
+        }
+
+        template <typename field_rnumber,
+                  field_backend be,
+                  field_components fc>
+        int sample(const field<field_rnumber, be, fc> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            this->pdistributor.template compute_distr<interpolator_class,
+                                                      field<field_rnumber, be, fc>,
+                                                      state_size,
+                                                      ncomp(fc)>(
+                    this->pinterpolator,
+                    field_to_sample,
+                    &this->local_number_of_particles,
+                    &this->local_state,
+                    result,
+                    neighbours);
+            return EXIT_SUCCESS;
+        }
+
+        int redistribute()
+        {
+            // todo: ask why the last template parameter is needed
+            this->pdistributor.template redistribute<interpolator_class,
+                                                     state_size,
+                                                     state_size,
+                                                     1>(
+                    this->pinterpolator,
+                    this->number_particles_per_partition.get(),
+                    &this->local_number_of_particles,
+                    &this->local_state,
+                    NULL,
+                    0,
+                    &this->local_index);
+            return EXIT_SUCCESS;
+        }
+};
+
+#endif//PARTICLE_SET_HPP
+

cpp/particles/interpolation/particles_field_computer.hpp

@@ -28,6 +28,7 @@
 
 #include <array>
 #include <utility>
+#include <cmath>
 
 #include "scope_timer.hpp"
 #include "particles/particles_utils.hpp"

cpp/particles/interpolation/spline.hpp

@@ -2,20 +2,20 @@
 *                                                                             *
 *  Copyright 2017 Max Planck Institute for Dynamics and Self-Organization     *
 *                                                                             *
-*  This file is part of bfps.                                                 *
+*  This file is part of TurTLE.                                               *
 *                                                                             *
-*  bfps is free software: you can redistribute it and/or modify               *
+*  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.                                     *
 *                                                                             *
-*  bfps is distributed in the hope that it will be useful,                    *
+*  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 bfps.  If not, see <http://www.gnu.org/licenses/>               *
+*  along with TurTLE.  If not, see <http://www.gnu.org/licenses/>             *
 *                                                                             *
 * Contact: Cristian.Lalescu@ds.mpg.de                                         *
 *                                                                             *
@@ -26,15 +26,15 @@
 #ifndef SPLINE_HPP
 #define SPLINE_HPP
 
-#include "spline_n1.hpp"
-#include "spline_n2.hpp"
-#include "spline_n3.hpp"
-#include "spline_n4.hpp"
-#include "spline_n5.hpp"
-#include "spline_n6.hpp"
-#include "spline_n7.hpp"
-#include "spline_n8.hpp"
-#include "spline_n9.hpp"
-#include "spline_n10.hpp"
+#include "particles/interpolation/spline_n1.hpp"
+#include "particles/interpolation/spline_n2.hpp"
+#include "particles/interpolation/spline_n3.hpp"
+#include "particles/interpolation/spline_n4.hpp"
+#include "particles/interpolation/spline_n5.hpp"
+#include "particles/interpolation/spline_n6.hpp"
+#include "particles/interpolation/spline_n7.hpp"
+#include "particles/interpolation/spline_n8.hpp"
+#include "particles/interpolation/spline_n9.hpp"
+#include "particles/interpolation/spline_n10.hpp"
 
 #endif

cpp/particles/particle_solver.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 PARTICLE_SOLVER_HPP
+#define PARTICLE_SOLVER_HPP
+
+#include <cassert>
+#include "particles/interpolation/abstract_particle_set.hpp"
+#include "particles/abstract_particle_rhs.hpp"
+
+namespace particle_solver
+{
+    template<int state_size>
+    int Euler(
+            double dt,
+            abstract_particle_set &pset,
+            const abstract_particle_rhs &rhs)
+    {
+        auto *rhs_val = new abstract_particle_set::particle_rnumber[pset.getLocalNumberOfParticles()*state_size];
+        auto *xx = new abstract_particle_set::particle_rnumber[pset.getLocalNumberOfParticles()*state_size];
+        rhs(0, pset, rhs_val);
+        for (abstract_particle_set::partsize_t idx_part = 0; idx_part < pset.getLocalNumberOfParticles(); idx_part++)
+        {
+            for (unsigned int cc = 0; cc < state_size; cc++)
+            {
+                xx[idx_part*state_size + cc] = (
+                        pset.getParticleState()[idx_part*state_size + cc] +
+                        dt * rhs_val[idx_part*state_size + cc]);
+            }
+        }
+        // TODO: replace pset particle state with xx values
+        delete[] rhs_val;
+        delete[] xx;
+        return EXIT_SUCCESS;
+    }
+};
+
+#endif//PARTICLE_SOLVER_HPP
+

cpp/particles/particles_system.hpp

@@ -29,6 +29,7 @@
 #include <array>
 #include <set>
 #include <algorithm>
+#include <cmath>
 
 #include "particles/abstract_particles_system.hpp"
 #include "particles/abstract_particles_system_with_p2p.hpp"
@@ -85,10 +86,12 @@ class particles_system : public abstract_particles_system_with_p2p<partsize_t, r
     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,
+    particles_system(const std::array<size_t,3>& field_grid_dim,
+                     const std::array<real_number,3>& in_spatial_box_width,
                      const std::array<real_number,3>& in_spatial_box_offset,
                      const std::array<real_number,3>& in_spatial_partition_width,
-                     const real_number in_my_spatial_low_limit, const real_number in_my_spatial_up_limit,
+                     const real_number in_my_spatial_low_limit,
+                     const real_number in_my_spatial_up_limit,
                      const std::array<size_t,3>& in_local_field_dims,
                      const std::array<size_t,3>& in_local_field_offset,
                      const field_class& in_field,
@@ -99,19 +102,36 @@ public:
                      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[IDXC_Z], in_local_field_offset[IDXC_Z] + in_local_field_dims[IDXC_Z]}),
+          current_partition_interval({in_local_field_offset[IDXC_Z],
+                                      in_local_field_offset[IDXC_Z] + in_local_field_dims[IDXC_Z]}),
           partition_interval_size(current_partition_interval.second - current_partition_interval.first),
           interpolator(),
-          particles_distr(in_mpi_com, current_partition_interval,field_grid_dim),
+          particles_distr(in_mpi_com,
+                          current_partition_interval,
+                          field_grid_dim),
           positions_updater(),
-          computer(field_grid_dim, current_partition_interval,
-                   interpolator, in_spatial_box_width, in_spatial_box_offset, in_spatial_partition_width),
+          computer(field_grid_dim,
+                   current_partition_interval,
+                   interpolator,
+                   in_spatial_box_width,
+                   in_spatial_box_offset,
+                   in_spatial_partition_width),
           default_field(in_field),
-          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),
-          computer_p2p(std::move(in_computer_p2p)), computer_particules_inner(std::move(in_computer_particules_inner)){
+          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),
+          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]);