Skip to content
Snippets Groups Projects
Commit c3776e3c authored by Armin Roediger's avatar Armin Roediger Committed by Cristian Lalescu
Browse files

Added generation of random particles by cpp in case cpp_random_particles is switched on.

parent bc33d02a
No related branches found
No related tags found
1 merge request!149Added generation of random particles by cpp in case cpp_random_particles is switched on.
......@@ -27,6 +27,8 @@
#include "scope_timer.hpp"
#include "particles/p2p/p2p_ghost_collisions.hpp"
#include "particles/inner/particles_inner_computer_2nd_order.hpp"
#include "particles/interpolation/particle_set.hpp"
#include "particles/particles_input_random.hpp"
#include <cmath>
......@@ -35,6 +37,7 @@ 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],
......@@ -42,6 +45,64 @@ int NSVE_Stokes_particles<rnumber>::initialize(void)
this->fs->cvelocity->rlayout->comm,
this->fs->cvelocity->fftw_plan_rigor);
this->particles_output_writer_mpi = new particles_output_hdf5<
long long int, double, 6>(
MPI_COMM_WORLD,
"tracers0",
nparticles,
tracers0_integration_steps);
// if this is the first iteration, and we are supposed to generate our own random particles,
// we do that here, and then we write them to file for future reference.
if ((this->cpp_random_particles > 0) &&
this->fs->iteration == 0)
{
// temporary particle set
// interpolation setup is irrelevant for output
particle_set<6, 1, 1> pset(
this->fs->cvelocity->rlayout,
this->fs->kk->dkx,
this->fs->kk->dky,
this->fs->kk->dkz);
particles_input_random<long long int, double, 6> pinput(
this->comm,
this->nparticles,
this->cpp_random_particles, // seed
pset.getSpatialLowLimitZ(),
pset.getSpatialUpLimitZ());
// initialize particle set
pset.init(pinput);
// allocate dummy rhs data
// we must fill these arrays with 0
// otherwise we get floating point exceptions because data is being written/read
std::vector<std::unique_ptr<double[]>> rhs_data;
rhs_data.resize(tracers0_integration_steps);
for (int counter = 0; counter < tracers0_integration_steps; counter++)
{
rhs_data[counter].reset(new double[pset.getLocalNumberOfParticles()*6]); // memory allocated here will be freed outside of the "if" block
std::fill_n(rhs_data[counter].get(), pset.getLocalNumberOfParticles()*6, 0);
}
// create dummy file_layout object
std::vector<hsize_t> file_layout;
file_layout.resize(1);
file_layout[0] = this->nparticles;
// write data
this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
this->particles_output_writer_mpi->update_particle_species_name("tracers0");
this->particles_output_writer_mpi->setParticleFileLayout(file_layout);
#ifdef USE_TIMING_OUTPUT
// barrier so that timings within save are more accurate.
MPI_Barrier(this->comm);
#endif
this->particles_output_writer_mpi->template save<6>(
pset.getParticleState(),
rhs_data.data(),
pset.getParticleIndices(),
pset.getLocalNumberOfParticles(),
0);
this->particles_output_writer_mpi->close_file();
}
particles_inner_computer_2nd_order_Stokes<double, long long int> current_particles_inner_computer;
current_particles_inner_computer.set_drag_coefficient(this->drag_coefficient);
this->ps = particles_system_builder_with_p2p(
......@@ -59,13 +120,6 @@ int NSVE_Stokes_particles<rnumber>::initialize(void)
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, double, 3>(
......@@ -212,6 +266,7 @@ int NSVE_Stokes_particles<rnumber>::read_parameters(void)
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<long long int>(parameter_file, "parameters/nparticles");
this->cpp_random_particles = hdf5_tools::read_value<int>(parameter_file, "parameters/cpp_random_particles");
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");
......
......@@ -52,6 +52,7 @@ class NSVE_Stokes_particles: public NSVE<rnumber>
int tracers0_integration_steps;
int tracers0_neighbours;
int tracers0_smoothness;
int cpp_random_particles;
double tracers0_cutoff;
double drag_coefficient;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment