-
Chichi Lalescu authoredChichi Lalescu authored
NSVEparticles.cpp 6.28 KiB
#include <string>
#include <cmath>
#include "NSVEparticles.hpp"
#include "scope_timer.hpp"
template <typename rnumber>
int NSVEparticles<rnumber>::initialize(void)
{
TIMEZONE("NSVEparticles::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);
this->ps = particles_system_builder(
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);
this->particles_output_writer_mpi = new particles_output_hdf5<
long long int, double, 3>(
MPI_COMM_WORLD,
"tracers0",
nparticles,
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 + "_particles.h5"),
"tracers0",
"position/0");
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEparticles<rnumber>::step(void)
{
TIMEZONE("NSVEparticles::step");
this->fs->compute_velocity(this->fs->cvorticity);
this->fs->cvelocity->ift();
this->ps->completeLoop(this->dt);
this->NSVE<rnumber>::step();
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEparticles<rnumber>::write_checkpoint(void)
{
TIMEZONE("NSVEparticles::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<3>(
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 NSVEparticles<rnumber>::finalize(void)
{
TIMEZONE("NSVEparticles::finalize");
delete this->pressure;
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 NSVEparticles<rnumber>::do_stats()
{
TIMEZONE("NSVEparticles::do_stats");
/// fluid stats go here
this->NSVE<rnumber>::do_stats();
if (!(this->iteration % this->niter_part == 0))
return EXIT_SUCCESS;
// allocate temporary data array
std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
// copy position data
/// sample position
std::copy(this->ps->getParticlesState(),
this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
pdata.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"position",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
/// sample velocity
std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"velocity",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
/// compute acceleration and sample it
this->fs->compute_Lagrangian_acceleration(this->tmp_vec_field, this->pressure);
this->tmp_vec_field->ift();
std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->tmp_vec_field, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"acceleration",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
// deallocate temporary data array
// TODO: is it required/safe to call the release method here?
//pdata.release();
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEparticles<rnumber>::read_parameters(void)
{
TIMEZONE("NSVEparticles::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->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 class NSVEparticles<float>;
template class NSVEparticles<double>;