Commit e9be1833 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

adds option for cpp random particle initialization

parent e4b8c9be
......@@ -140,6 +140,7 @@ class DNS(_code):
self.NSVEp_extra_parameters['niter_part_fine_period'] = int(10)
self.NSVEp_extra_parameters['niter_part_fine_duration'] = int(0)
self.NSVEp_extra_parameters['nparticles'] = int(10)
self.NSVEp_extra_parameters['cpp_random_particles'] = int(0)
self.NSVEp_extra_parameters['tracers0_integration_steps'] = int(4)
self.NSVEp_extra_parameters['tracers0_neighbours'] = int(1)
self.NSVEp_extra_parameters['tracers0_smoothness'] = int(1)
......@@ -914,7 +915,7 @@ class DNS(_code):
if self.dns_type == 'NSVE_Stokes_particles':
dset[cc*batch_size:(cc+1)*batch_size, 3:] = self.parameters['initial_particle_vel']*get_random_versors(batch_size)
else:
dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
nn -= batch_size
else:
dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
......@@ -1057,10 +1058,11 @@ class DNS(_code):
def generate_particle_data(
self,
opt = None):
if self.parameters['nparticles'] > 0:
self.generate_tracer_state(
species = 0,
rseed = opt.particle_rand_seed)
if (self.parameters['nparticles'] > 0):
if (self.parameters['cpp_random_particles'] == 0):
self.generate_tracer_state(
species = 0,
rseed = opt.particle_rand_seed)
if not os.path.exists(self.get_particle_file_name()):
with h5py.File(self.get_particle_file_name(), 'w') as particle_file:
particle_file.create_group('tracers0/position')
......@@ -1136,7 +1138,7 @@ class DNS(_code):
data = self.generate_vector_field(
write_to_file = False,
spectra_slope = 2.0,
amplitude = 0.05)
amplitude = 0.05)
f['vorticity/complex/{0}'.format(0)] = data
f.close()
if self.dns_type == 'kraichnan_field':
......
......@@ -28,6 +28,8 @@
#include <cmath>
#include "NSVEparticles.hpp"
#include "scope_timer.hpp"
#include "particles/interpolation/particle_set.hpp"
#include "particles/particles_input_random.hpp"
template <typename rnumber>
int NSVEparticles<rnumber>::initialize(void)
......@@ -41,6 +43,63 @@ int NSVEparticles<rnumber>::initialize(void)
this->fs->cvelocity->rlayout->comm,
this->fs->cvelocity->fftw_plan_rigor);
// initialize particle output object
// this may be needed before the particle system,
// if cpp_random_particles is true
this->particles_output_writer_mpi = new particles_output_hdf5<
long long int, double, 3>(
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<3, 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, 3> 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()*3]);
std::fill_n(rhs_data[counter].get(), pset.getLocalNumberOfParticles()*3, 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);
this->particles_output_writer_mpi->template save<3>(
pset.getParticleState(),
rhs_data.data(),
pset.getParticleIndex(),
pset.getLocalNumberOfParticles(),
0);
this->particles_output_writer_mpi->close_file();
}
DEBUG_MSG_WAIT(MPI_COMM_WORLD, "about to call particles_system_builder\n");
this->ps = particles_system_builder(
this->fs->cvelocity, // (field object)
......@@ -55,13 +114,8 @@ int NSVEparticles<rnumber>::initialize(void)
this->comm,
this->fs->iteration+1);
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, 3>(
MPI_COMM_WORLD,
"tracers0",
nparticles,
tracers0_integration_steps);
this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
// initialize sample object
this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
long long int, double, 3>(
MPI_COMM_WORLD,
......@@ -69,6 +123,9 @@ int NSVEparticles<rnumber>::initialize(void)
(this->simname + "_particles.h5"),
"tracers0",
"position/0");
// set particle file layout for output objects
this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
return EXIT_SUCCESS;
}
......@@ -198,6 +255,7 @@ int NSVEparticles<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<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");
......
......@@ -57,6 +57,7 @@ class NSVEparticles: public NSVE<rnumber>
int tracers0_integration_steps;
int tracers0_neighbours;
int tracers0_smoothness;
int cpp_random_particles;
/* other stuff */
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment