adding collision test to the pipeline

This adds the test of collisions and the corresponding particle to particle interaction to the pipeline.

TurTLE/test/collisions/NSVE_Stokes_growing_subset.cpp

+/**********************************************************************
+*                                                                     *
+*  Copyright 2022 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@mpcdf.mpg.de                              *
+*                                                                     *
+**********************************************************************/
+
+
+
+#include "scope_timer.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include <string>
+#include <cmath>
+
+template <typename rnumber>
+int NSVE_Stokes_growing_subset<rnumber>::initialize(void)
+{
+    TIMEZONE("NSVE_Stokes_growing_subset::intialize");
+    this->NSVE<rnumber>::initialize();
+    this->fti = new field_tinterpolator<rnumber, FFTW, THREE, NONE>();
+    this->srhs = new stokes_collisions_with_background_rhs<rnumber, FFTW, NONE>();
+    this->srhs->setDragCoefficient(0.1);
+    this->srhs->setStateSize(7);
+    this->srhs->getCollisionCounter().set_density_ratio(0.001);
+    this->srhs->getCollisionCounter().set_nu(this->fs->nu);
+    this->srhs->getCollisionCounter().set_dt(this->dt);
+    this->srhs->getCollisionCounter().set_nb_growing_particles(this->nb_growing_particles);
+    // We're not using Adams-Bashforth in this code.
+    // Euler integration for now, but this assert is there to ensure
+    // user knows what's happening.
+    assert(tracers0_integration_steps == 0);
+    // neighbours and smoothness are second and third template parameters of particle_set
+    assert(tracers0_neighbours == 3);
+    assert(tracers0_smoothness == 2);
+    this->pset = new particle_set<7, 3, 2>(
+            this->fs->cvelocity->rlayout,
+            this->dkx,
+            this->dky,
+            this->dkz,
+            this->tracers0_cutoff);
+
+    this->fti->set_field(this->fs->cvelocity);
+    this->srhs->setVelocity(this->fti);
+
+    particles_input_hdf5<long long int, double, 7, 7> particle_reader(
+            this->comm,
+            this->fs->get_current_fname(),
+            std::string("/tracers0/state/") + std::to_string(this->fs->iteration), // dataset name for initial input
+            std::string("/tracers0/rhs/")  + std::to_string(this->fs->iteration),  // dataset name for initial input
+            this->pset->getSpatialLowLimitZ(),
+            this->pset->getSpatialUpLimitZ());
+    this->pset->init(particle_reader);
+    DEBUG_MSG("local number of particles after input %d\n",this->pset->getLocalNumberOfParticles());
+    
+    this->psolver = new particle_solver(*(this->pset), 0);
+
+    this->particles_output_writer_mpi = new particles_output_hdf5<
+        long long int, double, 7>(
+                MPI_COMM_WORLD,
+                "tracers0",
+                nparticles,
+                tracers0_integration_steps);
+    this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+        long long int, double, double, 3>(
+                MPI_COMM_WORLD,
+                this->pset->getTotalNumberOfParticles(),
+                (this->simname + "_particles.h5"),
+                "tracers0",
+                "position/0");
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_growing_subset<rnumber>::step(void)
+{
+    TIMEZONE("NSVE_Stokes_growing_subset::step");
+    DEBUG_MSG("step() is called\n");
+    this->fs->compute_velocity(this->fs->cvorticity);
+    this->fs->cvelocity->ift();
+    DEBUG_MSG("real velocity field is computed\n");
+    this->psolver->Euler(this->dt, *(this->srhs));
+    DEBUG_MSG("particle timestep is performed\n");
+    DEBUG_MSG("local_number of collision pairs: %d\n", (this->srhs->getCollisionCounter()).get_local_collision_nb());
+    this->NSVE<rnumber>::step();
+    DEBUG_MSG("field step() is called\n");
+    this->psolver->setIteration(this->fs->iteration);
+    DEBUG_MSG("step is cpompleted\n");
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_growing_subset<rnumber>::write_checkpoint(void)
+{
+    TIMEZONE("NSVE_Stokes_growing_subset::write_checkpoint");
+    this->NSVE<rnumber>::write_checkpoint();
+    this->psolver->setIteration(this->fs->iteration);
+    this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+    this->psolver->template writeCheckpoint<7>(this->particles_output_writer_mpi);
+    this->particles_output_writer_mpi->close_file();
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_growing_subset<rnumber>::finalize(void)
+{
+    TIMEZONE("NSVE_Stokes_growing_subset::finalize");
+    delete this->particles_output_writer_mpi;
+    delete this->particles_sample_writer_mpi;
+    delete this->psolver;
+    delete this->pset;
+    delete this->srhs;
+    delete this->fti;
+    this->NSVE<rnumber>::finalize();
+    return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample fields at particle locations.
+ */
+
+template <typename rnumber>
+int NSVE_Stokes_growing_subset<rnumber>::do_stats()
+{
+    TIMEZONE("NSVE_Stokes_growing_subset::do_stats");
+    /// fluid stats go here
+    this->NSVE<rnumber>::do_stats();
+    this->srhs->getCollisionCounter().MPI_merge(
+            this->comm,
+            this->myrank,
+            this->nprocs);
+    if (this->myrank == 0)
+    {
+        // this->stat_file is only defined for rank 0,
+        // so only rank 0 can write the pairs.
+        // this is also a reason for calling MPI_merge beforehand.
+        hdf5_tools::write_particle_ID_pairs_with_single_rank<long long>(
+                this->srhs->getCollisionCounter().get_collision_pairs(
+                        this->comm,
+                        this->myrank,
+                        this->nprocs),
+                this->stat_file,
+                std::string("/statistics/collisions/pair_list_") + std::to_string(psolver->getIteration()).c_str());
+    }
+
+    /// either one of two conditions suffices to compute statistics:
+    /// 1) current iteration is a multiple of niter_part
+    /// 2) we are within niter_part_fine_duration/2 of a multiple of niter_part_fine_period
+    if (!(this->iteration % this->niter_part == 0))
+        return EXIT_SUCCESS;
+    // sample position
+    this->pset->writeStateTriplet(
+            0,
+            this->particles_sample_writer_mpi,
+            "tracers0",
+            "position",
+            psolver->getIteration());
+    if (!(this->iteration % this->niter_stat == 0))
+    {
+        // we need to compute velocity field manually, because it didn't happen in NSVE::do_stats()
+	this->fs->compute_velocity(this->fs->cvorticity);
+        *this->tmp_vec_field = this->fs->cvelocity->get_cdata();
+        this->tmp_vec_field->ift();
+    }
+    this->pset->writeSample(
+            this->tmp_vec_field,
+            this->particles_sample_writer_mpi,
+            "tracers0",
+            "velocity",
+            psolver->getIteration());
+    return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_growing_subset<rnumber>::read_parameters(void)
+{
+    TIMEZONE("NSVE_Stokes_growing_subset::read_parameters");
+    this->NSVE<rnumber>::read_parameters();
+    hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+    this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+    this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+    this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+    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_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/tracers0_cutoff");
+    this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
+    this->nb_growing_particles = hdf5_tools::read_value<int>(parameter_file, "parameters/nb_growing_particles");
+    H5Fclose(parameter_file);
+    MPI_Barrier(this->comm);
+    return EXIT_SUCCESS;
+}
+
+template class NSVE_Stokes_growing_subset<float>;
+template class NSVE_Stokes_growing_subset<double>;
+