Skip to content
Snippets Groups Projects
NSVE_Stokes_particles.cpp 10.34 KiB
/**********************************************************************
*                                                                     *
*  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 "NSVE_Stokes_particles.hpp"
#include "scope_timer.hpp"
#include "particles/particles_sampling.hpp"
#include "particles/p2p_ghost_collisions.hpp"
#include "particles/particles_inner_computer_2nd_order.hpp"

template <typename rnumber>
int NSVE_Stokes_particles<rnumber>::initialize(void)
{
    TIMEZONE("NSVE_Stokes_particles::intialize");
    this->NSVE<rnumber>::initialize();
    this->pressure = new field<rnumber, FFTW, ONE>(
            this->fs->cvelocity->rlayout->sizes[2],
            this->fs->cvelocity->rlayout->sizes[1],
            this->fs->cvelocity->rlayout->sizes[0],
            this->fs->cvelocity->rlayout->comm,
            this->fs->cvelocity->fftw_plan_rigor);

    particles_inner_computer_2nd_order_Stokes<double, long long int> current_particles_inner_computer;
    current_particles_inner_computer.set_drag_coefficient(this->drag_coefficient);
    //DEBUG_MSG("drag coefficient is set to %f \n", current_particles_inner_computer.get_drag_coefficient());
    //DEBUG_MSG_WAIT(MPI_COMM_WORLD, "before call to particles_system_builder\n");
    this->ps = particles_system_builder_with_p2p(
                this->fs->cvelocity,                                                    // (field object)
                this->fs->kk,                                                           // (kspace object, contains dkx, dky, dkz)
                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->fs->get_current_fname(),                                          // particles input filename
                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
                tracers0_neighbours,                                                    // parameter (interpolation no neighbours)
                tracers0_smoothness,                                                    // parameter
                this->comm,
                this->fs->iteration+1,
                std::move(p2p_ghost_collisions<double, long long int>()),
                std::move(current_particles_inner_computer),
                this->tracers0_cutoff);
    //DEBUG_MSG_WAIT(MPI_COMM_WORLD, "after call to particles_system_builder\n");
    this->particles_output_writer_mpi = new particles_output_hdf5<
        long long int, double, 6>(
                MPI_COMM_WORLD,
                "tracers0",
                nparticles,
                tracers0_integration_steps);
    this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
    this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
        long long int, double, 3>(
                MPI_COMM_WORLD,
                this->ps->getGlobalNbParticles(),
                (this->simname + "_particles.h5"),
                "tracers0",
                "position/0");
    this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
    //DEBUG_MSG("drag coefficient is after initialization %f \n", current_particles_inner_computer.get_drag_coefficient());
    return EXIT_SUCCESS;
}

template <typename rnumber>
int NSVE_Stokes_particles<rnumber>::step(void)
{
    TIMEZONE("NSVE_Stokes_particles::step");
    this->fs->compute_velocity(this->fs->cvorticity);
    this->fs->cvelocity->ift();
    this->ps->complete2ndOrderLoop(this->dt);
    this->NSVE<rnumber>::step();
    return EXIT_SUCCESS;
}

template <typename rnumber>
int NSVE_Stokes_particles<rnumber>::write_checkpoint(void)
{
    TIMEZONE("NSVE_Stokes_particles::write_checkpoint");
    this->NSVE<rnumber>::write_checkpoint();
    this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
    this->particles_output_writer_mpi->template save<6>(
            this->ps->getParticlesState(),
            this->ps->getParticlesRhs(),
            this->ps->getParticlesIndexes(),
            this->ps->getLocalNbParticles(),
            this->fs->iteration);
    this->particles_output_writer_mpi->close_file();
    return EXIT_SUCCESS;
}

template <typename rnumber>
int NSVE_Stokes_particles<rnumber>::finalize(void)
{
    TIMEZONE("NSVE_Stokes_particles::finalize");
    delete this->pressure;
    delete this->ps.release();
    delete this->particles_output_writer_mpi;
    delete this->particles_sample_writer_mpi;
    this->NSVE<rnumber>::finalize();
    return EXIT_SUCCESS;
}

/** \brief Compute fluid stats and sample fields at particle locations.
 */

template <typename rnumber>
int NSVE_Stokes_particles<rnumber>::do_stats()
{
    TIMEZONE("NSVE_Stokes_particles::do_stats");
    /// fluid stats go here
    this->NSVE<rnumber>::do_stats();


    /// 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 ||
          ((this->iteration + this->niter_part_fine_duration/2) % this->niter_part_fine_period <=
           this->niter_part_fine_duration)))
        return EXIT_SUCCESS;

    /// allocate temporary data array
    /// initialize pdata0 with the positions, and pdata1 with the orientations
    std::unique_ptr<double[]> pdata0 = this->ps->extractParticlesState(0, 3);
    std::unique_ptr<double[]> pdata1 = this->ps->extractParticlesState(3, 6);
    std::unique_ptr<double[]> pdata2(new double[9*this->ps->getLocalNbParticles()]);

    /// sample position
    this->particles_sample_writer_mpi->template save_dataset<3>(
            "tracers0",
            "position",
            pdata0.get(), // we need to use pdata0.get() here, because it's 3D, whereas getParticlesState may give something else
            &pdata0,
            this->ps->getParticlesIndexes(),
            this->ps->getLocalNbParticles(),
            this->ps->get_step_idx()-1);

    /// sample particle momentum
    this->particles_sample_writer_mpi->template save_dataset<3>(
            "tracers0",
            "momentum",
            pdata0.get(),
            &pdata1,
            this->ps->getParticlesIndexes(),
            this->ps->getLocalNbParticles(),
            this->ps->get_step_idx()-1);

    /// sample velocity
    std::fill_n(pdata1.get(), 3*this->ps->getLocalNbParticles(), 0);
    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->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get());
    this->particles_sample_writer_mpi->template save_dataset<3>(
            "tracers0",
            "velocity",
            pdata0.get(),
            &pdata1,
            this->ps->getParticlesIndexes(),
            this->ps->getLocalNbParticles(),
            this->ps->get_step_idx()-1);

    // deallocate temporary data array
    delete[] pdata0.release();
    delete[] pdata1.release();

    return EXIT_SUCCESS;
}

template <typename rnumber>
int NSVE_Stokes_particles<rnumber>::read_parameters(void)
{
    TIMEZONE("NSVE_Stokes_particles::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_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
    this->tracers0_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/tracers0_cutoff");
    this->drag_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/drag_coefficient");
    H5Fclose(parameter_file);
    return EXIT_SUCCESS;
}

template class NSVE_Stokes_particles<float>;
template class NSVE_Stokes_particles<double>;