From 6df4f5ac20ebef5c315a701e2088e32edd677508 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 3 Apr 2018 17:45:35 +0200
Subject: [PATCH] add test interpolation code that compiles

---
 bfps/cpp/full_code/test_interpolation.cpp | 162 ++++++++++++++++++++++
 bfps/cpp/full_code/test_interpolation.hpp |  53 +++++++
 setup.py                                  |   1 +
 3 files changed, 216 insertions(+)
 create mode 100644 bfps/cpp/full_code/test_interpolation.cpp
 create mode 100644 bfps/cpp/full_code/test_interpolation.hpp

diff --git a/bfps/cpp/full_code/test_interpolation.cpp b/bfps/cpp/full_code/test_interpolation.cpp
new file mode 100644
index 00000000..e93cd3eb
--- /dev/null
+++ b/bfps/cpp/full_code/test_interpolation.cpp
@@ -0,0 +1,162 @@
+#include "full_code/test_interpolation.hpp"
+
+
+template <typename rnumber>
+int test_interpolation<rnumber>::read_parameters(void)
+{
+    this->test::read_parameters();
+    hid_t parameter_file;
+    hid_t dset, memtype, space;
+    char fname[256];
+    hsize_t dims[1];
+    char *string_data;
+    sprintf(fname, "%s.h5", this->simname.c_str());
+    this->nparticles = hdf5_tools::read_value<int>(
+            parameter_file, "/parameters/nparticles");
+    this->tracers0_integration_steps = hdf5_tools::read_value<int>(
+            parameter_file, "/parameters/tracers0_integration_steps");
+    this->tracers0_neighbours = hdf5_tools::read_value<int>(
+            parameter_file, "/parameters/tracers0_neighbours");
+    this->tracers0_smoothness = hdf5_tools::read_value<int>(
+            parameter_file, "/parameters/tracers0_smoothness");
+    H5Fclose(parameter_file);
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int test_interpolation<rnumber>::initialize(void)
+{
+    this->test::initialize();
+
+    this->vorticity = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            DEFAULT_FFTW_FLAG);
+
+    this->velocity = new field<rnumber, FFTW, THREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            DEFAULT_FFTW_FLAG);
+
+    this->nabla_u = new field<rnumber, FFTW, THREExTHREE>(
+            this->nx, this->ny, this->nz,
+            this->comm,
+            DEFAULT_FFTW_FLAG);
+
+    this->kk = new kspace<FFTW, SMOOTH>(
+            this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+
+    if (this->myrank == 0)
+    {
+        hid_t stat_file = H5Fopen(
+                (this->simname + std::string(".h5")).c_str(),
+                H5F_ACC_RDWR,
+                H5P_DEFAULT);
+        this->kk->store(stat_file);
+        H5Fclose(stat_file);
+    }
+
+    this->ps = particles_system_builder(
+                this->velocity,                // (field object)
+                this->kk,                      // (kspace object, contains dkx, dky, dkz)
+                this->tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
+                (long long int)nparticles,      // to check coherency between parameters and hdf input file
+                this->simname + "_input.h5",    // particles input filename
+                std::string("/tracers0/state/0"), // dataset name for initial input
+                std::string("/tracers0/rhs/0")  , // dataset name for initial input
+                this->tracers0_neighbours,        // parameter (interpolation no neighbours)
+                this->tracers0_smoothness,        // parameter
+                this->comm,
+                1);
+    this->particles_output_writer_mpi = new particles_output_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                this->tracers0_integration_steps);
+    this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+        long long int, double, 3>(
+                MPI_COMM_WORLD,
+                this->ps->getGlobalNbParticles(),
+                (this->simname + "_output.h5"),
+                "tracers0",
+                "position/0");
+
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int test_interpolation<rnumber>::finalize(void)
+{
+    delete this->nabla_u;
+    delete this->velocity;
+    delete this->vorticity;
+    this->ps.release();
+    this->test::finalize();
+    delete this->kk;
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int test_interpolation<rnumber>::do_work()
+{
+    std::string fname = this->simname + std::string("_input.h5");
+    // read vorticity field
+    this->vorticity->io(
+            fname,
+            "vorticity",
+            0, true);
+
+    // compute velocity
+    invert_curl(this->kk, this->vorticity, this->velocity);
+
+    // compute velocity gradient
+    compute_gradient(this->kk, this->velocity, this->nabla_u);
+
+    // go to real space
+    this->vorticity->ift();
+    this->velocity->ift();
+    this->nabla_u->ift();
+
+    // allocate interpolation arrays
+    std::unique_ptr<double[]> p3data(new double[3*this->ps->getLocalNbParticles()]);
+    std::unique_ptr<double[]> p9data(new double[9*this->ps->getLocalNbParticles()]);
+
+    /// sample velocity at particles' position
+    this->ps->sample_compute_field(*this->velocity, p3data.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "velocity",
+            this->ps->getParticlesState(),
+            &p9data,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+    /// sample vorticity at particles' position
+    this->ps->sample_compute_field(*this->vorticity, p3data.get());
+    this->particles_sample_writer_mpi->template save_dataset<3>(
+            "tracers0",
+            "vorticity",
+            this->ps->getParticlesState(),
+            &p9data,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+    /// sample velocity gradient at particles' position
+    this->ps->sample_compute_field(*this->nabla_u, p9data.get());
+    this->particles_sample_writer_mpi->template save_dataset<9>(
+            "tracers0",
+            "velocity_gradient",
+            this->ps->getParticlesState(),
+            &p9data,
+            this->ps->getParticlesIndexes(),
+            this->ps->getLocalNbParticles(),
+            this->ps->get_step_idx()-1);
+
+    // no need to deallocate because we used "unique_ptr"
+    return EXIT_SUCCESS;
+}
+
+template class test_interpolation<float>;
+template class test_interpolation<double>;
+
diff --git a/bfps/cpp/full_code/test_interpolation.hpp b/bfps/cpp/full_code/test_interpolation.hpp
new file mode 100644
index 00000000..e8be247c
--- /dev/null
+++ b/bfps/cpp/full_code/test_interpolation.hpp
@@ -0,0 +1,53 @@
+#ifndef TEST_INTERPOLATION_HPP
+#define TEST_INTERPOLATION_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "full_code/test.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+
+/** \brief Interpolation tester.
+ *
+ */
+
+template <typename rnumber>
+class test_interpolation: public test
+{
+    public:
+        int nparticles;
+        int tracers0_integration_steps;
+        int tracers0_neighbours;
+        int tracers0_smoothness;
+
+        std::unique_ptr<abstract_particles_system<long long int, double>> ps;
+
+        particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
+        particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+
+        field<rnumber, FFTW, THREE> *velocity, *vorticity;
+        field<rnumber, FFTW, THREExTHREE> *nabla_u;
+
+        kspace<FFTW, SMOOTH> *kk;
+
+        test_interpolation(
+                const MPI_Comm COMMUNICATOR,
+                const std::string &simulation_name):
+            test(
+                    COMMUNICATOR,
+                    simulation_name){}
+        ~test_interpolation(){}
+
+        int initialize(void);
+        int do_work(void);
+        int finalize(void);
+
+        int read_parameters(void);
+};
+
+#endif//TEST_INTERPOLATION_HPP
+
diff --git a/setup.py b/setup.py
index 9511a92d..350bccf2 100644
--- a/setup.py
+++ b/setup.py
@@ -94,6 +94,7 @@ src_file_list = ['full_code/NSVEcomplex_particles',
                  'full_code/filter_test',
                  'full_code/field_test',
                  'full_code/field_output_test',
+                 'full_code/test_interpolation',
                  'hdf5_tools',
                  'full_code/get_rfields',
                  'full_code/resize',
-- 
GitLab