diff --git a/CMakeLists.txt b/CMakeLists.txt
index 89048a4ba41af69a7e35ceb76762576b7adc4911..aa5e8c13fe89ef2610e5bcf1a1581089737485ad 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -300,6 +300,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_Stokes_particles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.cpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -348,6 +349,7 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_Stokes_particles.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer_2nd_order.hpp
@@ -430,41 +432,56 @@ enable_testing()
if (BUILD_TESTING)
file(MAKE_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
### basic functionality
- add_test(NAME test_fftw
- COMMAND turtle.test_fftw
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_Parseval
- COMMAND turtle.test_Parseval
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_fftw
+ COMMAND turtle.test_fftw
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_Parseval
+ COMMAND turtle.test_Parseval
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
### compare DNS output to stored results
- add_test(NAME test_NSVEparticles
- COMMAND turtle.test_NSVEparticles
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_NSVEparticles
+ COMMAND turtle.test_NSVEparticles
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
### simple runs of post-processing tools
- add_test(NAME test_pp_single_to_double
- COMMAND turtle PP field_single_to_double --simname dns_nsveparticles --iter0 32 --iter1 32
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_pp_get_rfields
- COMMAND turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_pp_joint_acc_vel_stats
- COMMAND turtle PP joint_acc_vel_stats --simname dns_nsveparticles --iter0 0 --iter1 64
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_pp_resize
- COMMAND turtle PP resize --simname dns_nsveparticles --new_nx 96 --new_ny 96 --new_nz 96 --new_simname dns_nsveparticles_resized
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_single_to_double
+ COMMAND turtle PP field_single_to_double --simname dns_nsveparticles --iter0 32 --iter1 32
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_get_rfields
+ COMMAND turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_joint_acc_vel_stats
+ COMMAND turtle PP joint_acc_vel_stats --simname dns_nsveparticles --iter0 0 --iter1 64
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_resize
+ COMMAND turtle PP resize --simname dns_nsveparticles --new_nx 96 --new_ny 96 --new_nz 96 --new_simname dns_nsveparticles_resized
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
### simple runs of different DNS
- add_test(NAME test_NSVEp_extra_sampling
- COMMAND turtle DNS NSVEp_extra_sampling -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsvep_extra_sampling --nparticles 1000
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_NSVEcomplex_particles
- COMMAND turtle DNS NSVEcomplex_particles -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsvecomplex_particles --nparticles 1000
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_static_field
- COMMAND turtle DNS static_field --simname dns_static_field --nparticles 10000
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
- add_test(NAME test_kraichnan_field
- COMMAND turtle DNS kraichnan_field --simname dns_kraichnan_field --nparticles 10000
- WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_NSVEp_extra_sampling
+ COMMAND turtle DNS NSVEp_extra_sampling -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsvep_extra_sampling --nparticles 1000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_NSVEcomplex_particles
+ COMMAND turtle DNS NSVEcomplex_particles -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsvecomplex_particles --nparticles 1000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_static_field
+ COMMAND turtle DNS static_field --simname dns_static_field --nparticles 10000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_kraichnan_field
+ COMMAND turtle DNS kraichnan_field --simname dns_kraichnan_field --nparticles 10000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_NSVE_Stokes_particles
+ COMMAND turtle DNS NSVE_Stokes_particles -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsve_tokes_particles --nparticles 10000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
endif(BUILD_TESTING)
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index b8e4bd15ba1365db9bbca1a8f68ec6d8c50913e3..aa2a4bd29796252195f62870671289e08d94b8fc 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -439,7 +439,13 @@ class DNS(_code):
assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
assert (self.parameters['niter_todo'] % self.parameters['niter_out'] == 0)
assert (self.parameters['niter_out'] % self.parameters['niter_stat'] == 0)
- if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles', 'static_field', 'kraichnan_field']:
+ if self.dns_type in [
+ 'NSVEparticles_no_output',
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
+ 'NSVEparticles',
+ 'static_field',
+ 'kraichnan_field']:
assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
assert (self.parameters['niter_out'] % self.parameters['niter_part'] == 0)
_code.write_par(self, iter0 = iter0)
@@ -648,6 +654,10 @@ class DNS(_code):
'NSVEparticles',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
+ parser_NSVE_Stokes_particles = subparsers.add_parser(
+ 'NSVE_Stokes_particles',
+ help = 'plain Navier-Stokes vorticity formulation, with passive Stokes drag particles')
+
parser_NSVEp2p = subparsers.add_parser(
'NSVEcomplex_particles',
help = 'plain Navier-Stokes vorticity formulation, with oriented active particles')
@@ -660,6 +670,7 @@ class DNS(_code):
'NSVEparticles_no_output',
'NSVEp2',
'NSVEp2p',
+ 'NSVE_Stokes_particles',
'NSVEp_extra',
'static_field',
'kraichnan_field']:
@@ -673,6 +684,7 @@ class DNS(_code):
'NSVEparticles_no_output',
'NSVEp2',
'NSVEp2p',
+ 'NSVE_Stokes_particles',
'NSVEp_extra',
'static_field',
'kraichnan_field']:
@@ -726,7 +738,14 @@ class DNS(_code):
self.dns_type = opt.DNS_class
self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
# merge parameters if needed
- if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling', 'static_field', 'kraichnan_field']:
+ if self.dns_type in [
+ 'NSVEparticles',
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
+ 'NSVEparticles_no_output',
+ 'NSVEp_extra_sampling',
+ 'static_field',
+ 'kraichnan_field']:
for k in self.NSVEp_extra_parameters.keys():
self.parameters[k] = self.NSVEp_extra_parameters[k]
if type(extra_parameters) != type(None):
@@ -797,7 +816,14 @@ class DNS(_code):
# hardcoded FFTW complex representation size
field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
checkpoint_size = field_size
- if self.dns_type in ['kraichnan_field', 'static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
+ if self.dns_type in [
+ 'kraichnan_field',
+ 'static_field',
+ 'NSVEparticles',
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
+ 'NSVEparticles_no_output',
+ 'NSVEp_extra_sampling']:
rhs_size = self.parameters['tracers0_integration_steps']
if type(opt.tracers0_integration_steps) != type(None):
rhs_size = opt.tracers0_integration_steps
@@ -837,7 +863,9 @@ class DNS(_code):
integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps']
if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys():
integration_steps = self.parameters['tracers{0}_integration_steps'.format(species)]
- if self.dns_type == 'NSVEcomplex_particles' and species == 0:
+ if self.dns_type in [
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles'] and species == 0:
ncomponents = 6
with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
nn = self.parameters['nparticles']
@@ -1021,6 +1049,8 @@ class DNS(_code):
if self.dns_type in ['NSVEcomplex_particles']:
particle_file.create_group('tracers0/orientation')
particle_file.create_group('tracers0/velocity_gradient')
+ if self.dns_type in ['NSVE_Stokes_particles']:
+ particle_file.create_group('tracers0/momentum')
if self.dns_type in ['NSVEp_extra_sampling']:
particle_file.create_group('tracers0/velocity_gradient')
particle_file.create_group('tracers0/pressure')
@@ -1041,6 +1071,7 @@ class DNS(_code):
'static_field',
'NSVEparticles',
'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
'NSVEparticles_no_output',
'NSVEp_extra_sampling']:
if not os.path.exists(self.get_checkpoint_0_fname()):
@@ -1093,6 +1124,7 @@ class DNS(_code):
'static_field',
'NSVEparticles',
'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
'NSVEparticles_no_output',
'NSVEp_extra_sampling']:
self.generate_particle_data(opt = opt)
diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3606b162ce0be6b9904fa5a74ca63dcd1d011cf7
--- /dev/null
+++ b/cpp/full_code/NSVE_Stokes_particles.cpp
@@ -0,0 +1,215 @@
+/**********************************************************************
+* *
+* 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(0.1);
+
+ 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(particles_inner_computer_2nd_order_Stokes<double, long long int>()),
+ 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());
+ 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->completeLoop(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");
+ H5Fclose(parameter_file);
+ return EXIT_SUCCESS;
+}
+
+template class NSVE_Stokes_particles<float>;
+template class NSVE_Stokes_particles<double>;
+
diff --git a/cpp/full_code/NSVE_Stokes_particles.hpp b/cpp/full_code/NSVE_Stokes_particles.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..eec5247530a5b43647d847be262c99345495cd78
--- /dev/null
+++ b/cpp/full_code/NSVE_Stokes_particles.hpp
@@ -0,0 +1,89 @@
+/**********************************************************************
+* *
+* Copyright 2017 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 *
+* *
+**********************************************************************/
+
+
+
+#ifndef NSVE_STOKES_PARTICLES_HPP
+#define NSVE_STOKES_PARTICLES_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+
+/** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
+ *
+ * Child of Navier Stokes vorticity equation solver, this class calls all the
+ * methods from `NSVE`, and in addition integrates simple Lagrangian tracers
+ * in the resulting velocity field.
+ */
+
+template <typename rnumber>
+class NSVE_Stokes_particles: public NSVE<rnumber>
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ int niter_part;
+ int niter_part_fine_period;
+ int niter_part_fine_duration;
+ int nparticles;
+ int tracers0_integration_steps;
+ int tracers0_neighbours;
+ int tracers0_smoothness;
+ double tracers0_cutoff;
+
+
+ /* other stuff */
+ std::unique_ptr<abstract_particles_system<long long int, double>> ps;
+ field<rnumber, FFTW, ONE> *pressure;
+
+ particles_output_hdf5<long long int, double,6> *particles_output_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+
+
+ NSVE_Stokes_particles(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ ~NSVE_Stokes_particles(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//NSVE_STOKES_PARTICLES_HPP
+
diff --git a/cpp/particles/abstract_particles_system.hpp b/cpp/particles/abstract_particles_system.hpp
index 2f2f510f4bdad22b26b243607c7bbddcd2536771..e7f6cb250777ce7ca9498634fe0f8d2b329e9f17 100644
--- a/cpp/particles/abstract_particles_system.hpp
+++ b/cpp/particles/abstract_particles_system.hpp
@@ -61,6 +61,7 @@ public:
virtual void shift_rhs_vectors() = 0;
virtual void completeLoop(const real_number dt) = 0;
+ virtual void complete2ndOrderLoop(const real_number dt) = 0;
virtual void completeLoopWithVorticity(
const real_number dt,
diff --git a/cpp/particles/particles_inner_computer_2nd_order.hpp b/cpp/particles/particles_inner_computer_2nd_order.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..075914685ecb02d055947c8f1f0cd85fd8a97d26
--- /dev/null
+++ b/cpp/particles/particles_inner_computer_2nd_order.hpp
@@ -0,0 +1,84 @@
+/******************************************************************************
+* *
+* 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 *
+* *
+******************************************************************************/
+
+
+
+#ifndef PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
+#define PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
+
+#include <cstring>
+#include <cassert>
+
+template <class real_number, class partsize_t>
+class particles_inner_computer_2nd_order_Stokes{
+ double drag_coefficient;
+public:
+ template <int size_particle_state, int size_particle_rhs>
+ void compute_interaction(const partsize_t number_of_particles, real_number particle_state[], real_number particle_rhs[]) const{
+ }
+
+ template <int size_particle_state>
+ void enforce_unit_orientation(const partsize_t /*nb_particles*/, real_number /*pos_part*/[]) const{
+ }
+
+ template <int size_particle_state, int size_particle_rhs>
+ void add_Lagrange_multipliers(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[]) const{
+ }
+
+ template <int size_particle_state, int size_particle_rhs, int size_particle_rhs_extra>
+ void compute_interaction_with_extra(
+ const partsize_t number_of_particles,
+ real_number particle_state[],
+ real_number particle_rhs[],
+ const real_number sampled_velocity[]) const{
+ assert(size_particle_state == 6);
+ assert(size_particle_rhs == 6);
+ assert(size_particle_rhs_extra == 3);
+ #pragma omp parallel for
+ for(partsize_t idx_part = 0 ; idx_part < number_of_particles ; ++idx_part){
+ particle_rhs[idx_part*size_particle_rhs + IDXC_X] = particle_state[idx_part*size_particle_state + 3 + IDXC_X];
+ particle_rhs[idx_part*size_particle_rhs + IDXC_Y] = particle_state[idx_part*size_particle_state + 3 + IDXC_Y];
+ particle_rhs[idx_part*size_particle_rhs + IDXC_Z] = particle_state[idx_part*size_particle_state + 3 + IDXC_Z];
+ particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_X] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_X] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_X]);
+ particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Y] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Y] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Y]);
+ particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Z] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Z] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Z]);
+ }
+ }
+
+ constexpr static bool isEnable() {
+ return true;
+ }
+
+ void set_drag_coefficient(double mu)
+ {
+ this->drag_coefficient = mu;
+ }
+
+ double get_drag_coefficient()
+ {
+ return this->drag_coefficient;
+ }
+};
+
+#endif//PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
+
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index cdb523bb69de8c4b8d27ef319b33f29b8dc71db5..f276abfd6cbba028352027b43ebd98d0dfacc161 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -313,6 +313,21 @@ public:
shift_rhs_vectors();
}
+ void complete2ndOrderLoop(const real_number dt) final {
+ assert(size_particle_positions == 6);
+ assert(size_particle_rhs == 6);
+ std::unique_ptr<real_number[]> sampled_velocity(new real_number[getLocalNbParticles()*3]());
+ std::fill_n(sampled_velocity.get(), 3*getLocalNbParticles(), 0);
+ sample_compute_field(default_field, sampled_velocity.get());
+ this->computer_particules_inner.template compute_interaction_with_extra<size_particle_positions, size_particle_rhs, 3>(
+ my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get(),
+ sampled_velocity.get());
+ move(dt);
+ redistribute();
+ inc_step_idx();
+ shift_rhs_vectors();
+ }
+
void completeLoopWithVorticity(
const real_number dt,
const real_number particle_extra_field[]) final {
diff --git a/documentation/chapters/cpp_doxygen.rst b/documentation/chapters/cpp_doxygen.rst
index 1e4246a162083d4508f590d7370d2d29028cd598..51d1c9cca6c637fbce1335a3b8d6c358ddb6a2fd 100644
--- a/documentation/chapters/cpp_doxygen.rst
+++ b/documentation/chapters/cpp_doxygen.rst
@@ -21,3 +21,13 @@ CPP
.. doxygenclass:: NSVE
:project: TurTLE
:members:
+
+.. doxygenclass:: particles_distr_mpi
+ :project: TurTLE
+ :path: ...
+ :members: [...]
+ :protected-members:
+ :private-members:
+ :undoc-members:
+ :outline:
+ :no-link: