From 7e187b63992c81d7bb8b45e33baa073c2e86739c Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 12 May 2020 23:04:24 +0200
Subject: [PATCH] add basic test for tracer set solver

---
 CMakeLists.txt                                |   9 ++
 TurTLE/TEST.py                                |   6 +-
 cpp/full_code/test_tracer_set.cpp             | 102 ++++++++++++++++++
 cpp/full_code/test_tracer_set.hpp             |  60 +++++++++++
 cpp/particles/abstract_particle_rhs.hpp       |   4 +-
 .../interpolation/abstract_particle_set.hpp   |  14 +--
 .../interpolation/field_tinterpolator.hpp     |  14 +--
 cpp/particles/interpolation/particle_set.hpp  |  41 ++++++-
 cpp/particles/particle_solver.cpp             |   2 +-
 cpp/particles/particles_system.hpp            |   2 +-
 10 files changed, 230 insertions(+), 24 deletions(-)
 create mode 100644 cpp/full_code/test_tracer_set.cpp
 create mode 100644 cpp/full_code/test_tracer_set.hpp

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 7be5bcb1..0119e047 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -271,6 +271,7 @@ set(cpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/full_code/Gauss_field_test.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.cpp
+    ${PROJECT_SOURCE_DIR}/cpp/full_code/test_tracer_set.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/write_rpressure.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/bandpass_stats.cpp
@@ -299,6 +300,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/rhs/tracer_rhs.cpp
     ${PROJECT_SOURCE_DIR}/cpp/particles/particle_solver.cpp
     ${PROJECT_SOURCE_DIR}/cpp/scope_timer.cpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp
@@ -324,6 +326,7 @@ set(hpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/full_code/write_filtered_test.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.hpp
+    ${PROJECT_SOURCE_DIR}/cpp/full_code/test_tracer_set.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/write_rpressure.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/bandpass_stats.hpp
@@ -355,6 +358,7 @@ set(hpp_for_lib
     ${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/rhs/tracer_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
@@ -463,6 +467,11 @@ if (BUILD_TESTING)
         NAME test_Parseval
         COMMAND turtle.test_Parseval
         WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+    ### basic particle functionality
+    add_test(
+        NAME test_tracer_set
+        COMMAND turtle TEST test_tracer_set  -n 32 --np 2 --ntpp 2
+        WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
     ### compare DNS output to stored results
     add_test(
         NAME test_NSVEparticles
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 59757f64..e2c42bf9 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -366,6 +366,9 @@ class TEST(_code):
         parser_ornstein_uhlenbeck_test = subparsers.add_parser(
                 'ornstein_uhlenbeck_test',
                 help = 'test ornstein uhlenbeck process')
+        parser_test_tracer_set = subparsers.add_parser(
+                'test_tracer_set',
+                help = 'basic test of tracer set solver')
 
         for parser in ['Gauss_field_test',
                        'filter_test',
@@ -374,7 +377,8 @@ class TEST(_code):
                        'symmetrize_test',
                        'field_output_test',
                        'test_interpolation',
-                       'ornstein_uhlenbeck_test']:
+                       'ornstein_uhlenbeck_test',
+                       'test_tracer_set']:
             eval('self.simulation_parser_arguments(parser_' + parser + ')')
             eval('self.job_parser_arguments(parser_' + parser + ')')
             eval('self.parameters_to_parser_arguments(parser_' + parser + ')')
diff --git a/cpp/full_code/test_tracer_set.cpp b/cpp/full_code/test_tracer_set.cpp
new file mode 100644
index 00000000..92dae47b
--- /dev/null
+++ b/cpp/full_code/test_tracer_set.cpp
@@ -0,0 +1,102 @@
+/******************************************************************************
+*                                                                             *
+*  Copyright 2019 Max Planck Institute for Dynamics and Self-Organization     *
+*                                                                             *
+*  This file is part of bfps.                                                 *
+*                                                                             *
+*  bfps is free software: you can redistribute it and/or modify               *
+*  it under the terms of the GNU General Public License as published          *
+*  by the Free Software Foundation, either version 3 of the License,          *
+*  or (at your option) any later version.                                     *
+*                                                                             *
+*  bfps is distributed in the hope that it will be useful,                    *
+*  but WITHOUT ANY WARRANTY; without even the implied warranty of             *
+*  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the              *
+*  GNU General Public License for more details.                               *
+*                                                                             *
+*  You should have received a copy of the GNU General Public License          *
+*  along with bfps.  If not, see <http://www.gnu.org/licenses/>               *
+*                                                                             *
+* Contact: Cristian.Lalescu@ds.mpg.de                                         *
+*                                                                             *
+******************************************************************************/
+
+
+
+#include <string>
+#include <cmath>
+#include <random>
+#include "test_tracer_set.hpp"
+#include "particles/rhs/tracer_rhs.hpp"
+#include "particles/interpolation/particle_set.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int test_tracer_set<rnumber>::initialize(void)
+{
+    TIMEZONE("test_tracer_set::initialize");
+    this->read_parameters();
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int test_tracer_set<rnumber>::finalize(void)
+{
+    TIMEZONE("test_tracer_set::finalize");
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int test_tracer_set<rnumber>::read_parameters()
+{
+    TIMEZONE("test_tracer_set::read_parameters");
+    this->test::read_parameters();
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int test_tracer_set<rnumber>::do_work(void)
+{
+    TIMEZONE("test_tracer_set::do_work");
+    // allocate
+    field<rnumber, FFTW, THREE> *vec_field = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            FFTW_ESTIMATE);
+    field_tinterpolator<rnumber, FFTW, THREE, NONE> *fti = new field_tinterpolator<rnumber, FFTW, THREE, NONE>();
+    tracer_rhs<rnumber, FFTW, NONE> *trhs = new tracer_rhs<rnumber, FFTW, NONE>();
+
+    //std::array<double, 3> blad;
+    //blad[0] = 1.2;
+    //blad[1] = 1.2;
+    //blad[2] = 1.2;
+
+    //std::array<size_t, 3> blai;
+    //blai[0] = 32;
+    //blai[1] = 32;
+    //blai[2] = 32;
+
+    //particle_set<3, 2, 1> *pset = new particle_set<3, 2, 1>(
+    //        this->comm,
+    //        blad,
+    //        blad,
+    //        blad,
+    //        blai,
+    //        blai,
+    //        blai);
+
+    fti->set_field(vec_field);
+    trhs->setVelocity(fti);
+
+    // deallocate
+    //delete pset;
+    delete fti;
+    delete trhs;
+    delete vec_field;
+    return EXIT_SUCCESS;
+}
+
+template class test_tracer_set<float>;
+template class test_tracer_set<double>;
+
diff --git a/cpp/full_code/test_tracer_set.hpp b/cpp/full_code/test_tracer_set.hpp
new file mode 100644
index 00000000..b6d0ad4a
--- /dev/null
+++ b/cpp/full_code/test_tracer_set.hpp
@@ -0,0 +1,60 @@
+/**********************************************************************
+*                                                                     *
+*  Copyright 2017 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 TEST_TRACER_SET_HPP
+#define TEST_TRACER_SET_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "field.hpp"
+#include "full_code/test.hpp"
+
+/** \brief A class for testing basic `particle_set` functionality.
+ */
+
+template <typename rnumber>
+class test_tracer_set: public test
+{
+    public:
+        test_tracer_set(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            test(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~test_tracer_set(){}
+
+        int initialize(void);
+        int do_work(void);
+        int finalize(void);
+        int read_parameters(void);
+};
+
+#endif//TEST_TRACER_SET_HPP
+
diff --git a/cpp/particles/abstract_particle_rhs.hpp b/cpp/particles/abstract_particle_rhs.hpp
index a959a6a5..0a5fddc6 100644
--- a/cpp/particles/abstract_particle_rhs.hpp
+++ b/cpp/particles/abstract_particle_rhs.hpp
@@ -31,7 +31,7 @@
 // children of this class must provide a method to compute the rhs for a given particle set
 class abstract_particle_rhs
 {
-    private:
+    protected:
         using particle_rnumber = abstract_particle_set::particle_rnumber;
     public:
         // destructor
@@ -40,7 +40,7 @@ class abstract_particle_rhs
         // important bit
         virtual int operator()(
                 double t,
-                const abstract_particle_set &pset,
+                abstract_particle_set &pset,
                 particle_rnumber *result) const = 0;
 };
 
diff --git a/cpp/particles/interpolation/abstract_particle_set.hpp b/cpp/particles/interpolation/abstract_particle_set.hpp
index dfac8568..c452a50a 100644
--- a/cpp/particles/interpolation/abstract_particle_set.hpp
+++ b/cpp/particles/interpolation/abstract_particle_set.hpp
@@ -92,39 +92,39 @@ class abstract_particle_set
          * Optional: also redistribute a list of arrays of the same shape.
          */
         virtual int redistribute(
-                std::vector<std::unique_ptr<particle_rnumber>> &additional_data) = 0;
+                std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data) = 0;
 
         // sample field values at particle location
         virtual int sample(
                 const field<float,
                             FFTW,
                             ONE>         &field_to_sample,
-                particle_rnumber *result) const = 0;
+                particle_rnumber *result) = 0;
         virtual int sample(
                 const field<float,
                             FFTW,
                             THREE>       &field_to_sample,
-                particle_rnumber *result) const = 0;
+                particle_rnumber *result) = 0;
         virtual int sample(
                 const field<float,
                             FFTW,
                             THREExTHREE> &field_to_sample,
-                particle_rnumber *result) const = 0;
+                particle_rnumber *result) = 0;
         virtual int sample(
                 const field<double,
                             FFTW,
                             ONE>         &field_to_sample,
-                particle_rnumber *result) const = 0;
+                particle_rnumber *result) = 0;
         virtual int sample(
                 const field<double,
                             FFTW,
                             THREE>       &field_to_sample,
-                particle_rnumber *result) const = 0;
+                particle_rnumber *result) = 0;
         virtual int sample(
                 const field<double,
                             FFTW,
                             THREExTHREE> &field_to_sample,
-                particle_rnumber *result) const = 0;
+                particle_rnumber *result) = 0;
 };
 
 #endif//ABSTRACT_PARTICLE_SET_HPP
diff --git a/cpp/particles/interpolation/field_tinterpolator.hpp b/cpp/particles/interpolation/field_tinterpolator.hpp
index 5d77e336..dcf031a5 100644
--- a/cpp/particles/interpolation/field_tinterpolator.hpp
+++ b/cpp/particles/interpolation/field_tinterpolator.hpp
@@ -82,24 +82,24 @@ class field_tinterpolator
         // t is a fraction between 0 and 1.
         int operator()(
                 double t,
-                const abstract_particle_set *pset,
+                abstract_particle_set &pset,
                 particle_rnumber *result)
         {
             switch(tt)
             {
                 case NONE:
                     {
-                        pset->sample(*(this->field_list[0]), result);
+                        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++)
+                        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++)
                             {
diff --git a/cpp/particles/interpolation/particle_set.hpp b/cpp/particles/interpolation/particle_set.hpp
index 4a6bd4a6..5475c362 100644
--- a/cpp/particles/interpolation/particle_set.hpp
+++ b/cpp/particles/interpolation/particle_set.hpp
@@ -162,7 +162,7 @@ class particle_set: public abstract_particle_set
             return std::vector<hsize_t>(this->particle_file_layout);
         }
 
-        std::unique_ptr<particle_rnumber[]> extractParticlesState(
+        std::unique_ptr<particle_rnumber[]> extractFromParticleState(
                 const int firstState,
                 const int lastState) const
         {
@@ -182,7 +182,7 @@ class particle_set: public abstract_particle_set
         template <typename field_rnumber,
                   field_backend be,
                   field_components fc>
-        int sample(const field<field_rnumber, be, fc> &field_to_sample,
+        int template_sample(const field<field_rnumber, be, fc> &field_to_sample,
                    particle_rnumber *result)
         {
             this->pdistributor.template compute_distr<interpolator_class,
@@ -191,14 +191,45 @@ class particle_set: public abstract_particle_set
                                                       ncomp(fc)>(
                     this->pinterpolator,
                     field_to_sample,
-                    &this->local_number_of_particles,
-                    &this->local_state,
+                    this->number_particles_per_partition.get(),
+                    this->local_state.get(),
                     result,
                     neighbours);
             return EXIT_SUCCESS;
         }
 
-        int redistribute(std::vector<std::unique_ptr<particle_rnumber[]>> additional_data)
+        int sample(const field<float, FFTW, ONE> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            return template_sample<float, FFTW, ONE>(field_to_sample, result);
+        }
+        int sample(const field<float, FFTW, THREE> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            return template_sample<float, FFTW, THREE>(field_to_sample, result);
+        }
+        int sample(const field<float, FFTW, THREExTHREE> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            return template_sample<float, FFTW, THREExTHREE>(field_to_sample, result);
+        }
+        int sample(const field<double, FFTW, ONE> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            return template_sample<double, FFTW, ONE>(field_to_sample, result);
+        }
+        int sample(const field<double, FFTW, THREE> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            return template_sample<double, FFTW, THREE>(field_to_sample, result);
+        }
+        int sample(const field<double, FFTW, THREExTHREE> &field_to_sample,
+                   particle_rnumber *result)
+        {
+            return template_sample<double, FFTW, THREExTHREE>(field_to_sample, result);
+        }
+
+        int redistribute(std::vector<std::unique_ptr<particle_rnumber[]>> &additional_data)
         {
             this->pdistributor.template redistribute<interpolator_class,
                                                      state_size,
diff --git a/cpp/particles/particle_solver.cpp b/cpp/particles/particle_solver.cpp
index 7950bef4..e90e1673 100644
--- a/cpp/particles/particle_solver.cpp
+++ b/cpp/particles/particle_solver.cpp
@@ -44,7 +44,7 @@ int particle_solver::Euler(
                 });
 
         // update particle set
-        std::vector<std::unique_ptr<particle_rnumber>> temporary;
+        std::vector<std::unique_ptr<particle_rnumber[]>> temporary;
         this->pset->setParticleState(xx);
         this->pset->redistribute(temporary);
 
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index 2be42334..2f21ac42 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -361,7 +361,7 @@ public:
     //- Not generic to enable sampling begin
     void sample_compute_field(const field<float, FFTW, ONE>& sample_field,
                                 real_number sample_rhs[]) final {
-        // sample_compute<decltype(sample_field), 1>(sample_field, sample_rhs);
+        sample_compute<decltype(sample_field), 1>(sample_field, sample_rhs);
     }
     void sample_compute_field(const field<float, FFTW, THREE>& sample_field,
                                 real_number sample_rhs[]) final {
-- 
GitLab