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>;