diff --git a/CMakeLists.txt b/CMakeLists.txt
index eb67ee611dd7a61f8c9bea64ad247602510f77ef..60ac7739add4250ee20e90f6e032631fb88f2e8d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -299,6 +299,7 @@ set(cpp_for_lib
${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/particles/particle_solver.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
diff --git a/cpp/particles/abstract_particle_rhs.hpp b/cpp/particles/abstract_particle_rhs.hpp
index 92f94477eeebea57e569f253677db42a6d673188..a959a6a5414cec60f92eea487e7b25322e328f85 100644
--- a/cpp/particles/abstract_particle_rhs.hpp
+++ b/cpp/particles/abstract_particle_rhs.hpp
@@ -31,6 +31,8 @@
// children of this class must provide a method to compute the rhs for a given particle set
class abstract_particle_rhs
{
+ private:
+ using particle_rnumber = abstract_particle_set::particle_rnumber;
public:
// destructor
virtual ~abstract_particle_rhs(){}
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index da1ce30b62f4b32b98c4601564109c3fa2e62970..dfac856804ebfbdf8d08968be38aa85ab51b8e37 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -29,26 +29,36 @@
#include <array>
#include <cassert>
#include <vector>
+#include <memory>
#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.
+
+
+/** Brings together particle information with interpolation functionality.
+ *
+ * This is an abstract class that defines the functionality required by the
+ * particle solver code.
+ * We define methods for:
+ * - accessing particle set data and miscelaneous information.
+ * - updating particle information and redistributing data among MPI processes
+ * accordingly.
+ * - approximating fields at particle locations ("sampling").
+ *
+ */
+
class abstract_particle_set
{
- protected:
+ public:
using partsize_t = long long int;
using particle_rnumber = double;
- public:
- // destructor
+ // virtual destructor
virtual ~abstract_particle_set(){}
// extract particle information
virtual const particle_rnumber* getParticleState() const = 0;
virtual const partsize_t* getParticleIndex() const = 0;
+ virtual int setParticleState(particle_rnumber *) = 0;
virtual std::unique_ptr<particle_rnumber[]> extractFromParticleState(
const int firstComponent,
@@ -56,10 +66,33 @@ class abstract_particle_set
virtual partsize_t getLocalNumberOfParticles() const = 0;
virtual partsize_t getTotalNumberOfParticles() const = 0;
- virtual int getStateSize() const = 0;
+ virtual const int getStateSize() const = 0;
+
+ // generic loop over data
+ template <class func_type>
+ int LOOP(func_type expression)
+ {
+ for (partsize_t idx = 0; idx < this->getLocalNumberOfParticles()*this->getStateSize(); idx++)
+ expression(idx);
+ return EXIT_SUCCESS;
+ }
+ template <class func_type>
+ int LOOP_state(func_type expression)
+ {
+ for (partsize_t idx_part = 0; idx_part < this->getLocalNumberOfParticles(); idx_part++)
+ for (unsigned int cc = 0; cc < this->getStateSize(); cc++)
+ expression(idx_part, cc);
+ return EXIT_SUCCESS;
+ }
- // reorganize particles within MPI domain
- virtual int redistribute() = 0;
+ /** Reorganize particles within MPI domain.
+ *
+ * Based on current particle locations, redistribute the full state
+ * data among the MPI processes.
+ * Optional: also redistribute a list of arrays of the same shape.
+ */
+ virtual int redistribute(
+ std::vector<std::unique_ptr<particle_rnumber>> &additional_data) = 0;
// sample field values at particle location
virtual int sample(
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 55045e091750387a7bff1fe13e8188120adf97bc..84b958d0e4f55192c94f3da0a69a6e9540d53879 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -29,6 +29,7 @@
#include <array>
#include <cassert>
#include <vector>
+#include <memory>
#include "field.hpp"
#include "particles/particles_distr_mpi.hpp"
#include "particles/particles_output_hdf5.hpp"
@@ -36,6 +37,15 @@
#include "particles/interpolation/particles_generic_interp.hpp"
#include "particles/interpolation/abstract_particle_set.hpp"
+
+
+/** Practical implementation of particle sets.
+ *
+ * Child of `abstract_particle_set`. Brings together the functionality of
+ * `particles_field_computer` and `particles_distr_mpi`.
+ *
+ */
+
template <int state_size,
int neighbours,
int smoothness>
@@ -110,6 +120,15 @@ class particle_set: public abstract_particle_set
return this->local_state.get();
}
+ int setParticleState(particle_rnumber *src_local_state)
+ {
+ this->LOOP(
+ [&](const partsize_t idx){
+ this->local_state[idx] = src_local_state[idx];
+ });
+ return EXIT_SUCCESS;
+ }
+
const partsize_t* getParticleIndex() const
{
return this->local_index.get();
@@ -125,6 +144,11 @@ class particle_set: public abstract_particle_set
return this->total_number_of_particles;
}
+ const int getStateSize() const
+ {
+ return state_size;
+ }
+
int setParticleFileLayout(std::vector<hsize_t> input_layout)
{
this->particle_file_layout.resize(input_layout.size());
@@ -174,9 +198,8 @@ class particle_set: public abstract_particle_set
return EXIT_SUCCESS;
}
- int redistribute()
+ int redistribute(std::vector<std::unique_ptr<particle_rnumber>> additional_data)
{
- // todo: ask why the last template parameter is needed
this->pdistributor.template redistribute<interpolator_class,
state_size,
state_size,
@@ -185,8 +208,8 @@ class particle_set: public abstract_particle_set
this->number_particles_per_partition.get(),
&this->local_number_of_particles,
&this->local_state,
- NULL,
- 0,
+ additional_data.data(),
+ int(additional_data.size()),
&this->local_index);
return EXIT_SUCCESS;
}
diff --git a/cpp/particles/interpolation/particles_field_computer.hpp b/cpp/particles/interpolation/particles_field_computer.hpp
index 620b4685fdc140c999a9bce13784f4dc1d87a773..a80a37ceabedfd9e60f2bbd11ea5b11ab284560b 100644
--- a/cpp/particles/interpolation/particles_field_computer.hpp
+++ b/cpp/particles/interpolation/particles_field_computer.hpp
@@ -33,6 +33,12 @@
#include "scope_timer.hpp"
#include "particles/particles_utils.hpp"
+
+
+/** \brief Computes interpolation kernel sum operation.
+ *
+ * */
+
template <class partsize_t,
class real_number,
class interpolator_class,
diff --git a/cpp/particles/interpolation/particles_generic_interp.hpp b/cpp/particles/interpolation/particles_generic_interp.hpp
index e75ae9b3d92b34b170f13bde065eb24afa9dafe8..4966d5fdf66a434a3a51936cc3f00f72cc4d6adc 100644
--- a/cpp/particles/interpolation/particles_generic_interp.hpp
+++ b/cpp/particles/interpolation/particles_generic_interp.hpp
@@ -26,6 +26,10 @@
#ifndef PARTICLES_GENERIC_INTERP_HPP
#define PARTICLES_GENERIC_INTERP_HPP
+/** \brief Computes weights for interpolation formulas.
+ *
+ * */
+
template <class real_number, int interp_neighbours, int smoothness>
class particles_generic_interp;
diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7950bef404d3a64e57ee1f0315af6d903abdc866
--- /dev/null
+++ b/cpp/particles/particle_solver.cpp
@@ -0,0 +1,56 @@
+/******************************************************************************
+* *
+* 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/particle_solver.hpp"
+
+int particle_solver::Euler(
+ double dt,
+ const abstract_particle_rhs &rhs)
+ {
+ // temporary array for storing right-hand-side
+ auto *rhs_val = new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())];
+ // temporary array for new particle state
+ auto *xx = new particle_rnumber[this->pset->getLocalNumberOfParticles()*(this->pset->getStateSize())];
+
+ // compute right hand side
+ rhs(0, *(this->pset), rhs_val);
+
+ // compute new particle state
+ this->pset->LOOP(
+ [&](const partsize_t idx){
+ xx[idx] = this->pset->getParticleState()[idx] + dt * rhs_val[idx];
+ });
+
+ // update particle set
+ std::vector<std::unique_ptr<particle_rnumber>> temporary;
+ this->pset->setParticleState(xx);
+ this->pset->redistribute(temporary);
+
+ // deallocate temporary arrays
+ delete[] rhs_val;
+ delete[] xx;
+ return EXIT_SUCCESS;
+ }
+
diff --git a/cpp/particles/particle_solver.hpp b/cpp/particles/particle_solver.hpp
index a8ecff83fff82faec97f670b73e426c72880760a..e46677ba357e7436e342364d28a0fd3ed581f94c 100644
--- a/cpp/particles/particle_solver.hpp
+++ b/cpp/particles/particle_solver.hpp
@@ -30,31 +30,36 @@
#include "particles/interpolation/abstract_particle_set.hpp"
#include "particles/abstract_particle_rhs.hpp"
-namespace particle_solver
+
+
+/** Time-stepping and checkpointing functionality for particle systems.
+ *
+ * This class implements the universal functionality of a particle solver:
+ * time-stepping methods and checkpointing, together with the associated
+ * bookkeeping.
+ *
+ * We rely on predefined functionality of particle set objects, as well as
+ * "particle right-hand-side computers".
+ */
+
+class 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;
- }
+ private:
+ using particle_rnumber = abstract_particle_set::particle_rnumber;
+ using partsize_t = abstract_particle_set::partsize_t;
+
+ abstract_particle_set *pset;
+
+ public:
+ particle_solver(
+ abstract_particle_set &in_pset):
+ pset(&in_pset)
+ {}
+ ~particle_solver(){}
+
+ int Euler(
+ double dt,
+ const abstract_particle_rhs &rhs);
};
#endif//PARTICLE_SOLVER_HPP