Merge branch 'feature/random-particles' into develop

TurTLE/DNS.py

@@ -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)
@@ -1057,7 +1058,8 @@ class DNS(_code):
     def generate_particle_data(
             self,
             opt = None):
-        if self.parameters['nparticles'] > 0:
+        if (self.parameters['nparticles'] > 0):
+            if (self.parameters['cpp_random_particles'] == 0):
                 self.generate_tracer_state(
                         species = 0,
                         rseed = opt.particle_rand_seed)

cpp/full_code/NSVEparticles.cpp

@@ -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");

cpp/full_code/NSVEparticles.hpp

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

cpp/particles/particles_system_builder.hpp

@@ -40,7 +40,7 @@
 #include "kspace.hpp"
 
 const int MAX_INTERPOLATION_NEIGHBOURS=5;
-const int MAX_INTERPOLATION_SMOOTHNESS=2;
+const int MAX_INTERPOLATION_SMOOTHNESS=3;
 
 //////////////////////////////////////////////////////////////////////////////
 ///

tests/DNS/test_scaling.py

@@ -57,6 +57,7 @@ def get_DNS_parameters(
         DNS_parameters += [
                 '--precision', 'double']
     # check that source sim exists
+    print('looking for ', os.path.join(src_dirname, src_simname + '.h5'))
     assert(os.path.exists(os.path.join(src_dirname, src_simname + '.h5')))
     # check that source checkpoint exists
     dns_src = TurTLE.DNS(simname = src_simname, work_dir = src_dirname)
@@ -87,7 +88,7 @@ def get_DNS_parameters(
         DNS_parameters += [
                 '--tracers0_neighbours', '{0}'.format(nneighbours),
                 '--tracers0_smoothness', '{0}'.format(smoothness),
-                '--particle-rand-seed', '2']
+                '--cpp_random_particles', '2']
     if no_submit:
         DNS_parameters += ['--no-submit']
     DNS_parameters += ['--environment', environment,