From 3655957b6da5e720264c4ebd086f2b024021bc79 Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@ds.mpg.de>
Date: Tue, 27 Feb 2018 10:51:54 +0100
Subject: [PATCH] [broken] sample velocity gradient for complex particles
---
bfps/DNS.py | 24 +++++-----
bfps/cpp/full_code/NSVEcomplex_particles.cpp | 47 ++++++++++++--------
bfps/cpp/full_code/NSVEcomplex_particles.hpp | 21 +++++----
bfps/test/test_particles.py | 2 +-
setup.py | 2 +-
5 files changed, 56 insertions(+), 40 deletions(-)
diff --git a/bfps/DNS.py b/bfps/DNS.py
index 90988019..7eb15af6 100644
--- a/bfps/DNS.py
+++ b/bfps/DNS.py
@@ -127,7 +127,7 @@ class DNS(_code):
template_class = '{0}<{1}>::'.format(self.dns_type, rnumber),
template_prefix = 'template '.format(rnumber),
just_declaration = True) + '\n\n')
- if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output', 'NSVEparticlesP2P']:
+ if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output', 'NSVEcomplex_particles']:
outfile.write('template <typename rnumber> int NSVE<rnumber>::read_parameters(){return EXIT_SUCCESS;}\n')
outfile.write('template int NSVE<float>::read_parameters();\n')
outfile.write('template int NSVE<double>::read_parameters();\n\n')
@@ -167,9 +167,9 @@ class DNS(_code):
self.NSVEp_extra_parameters['tracers0_neighbours'] = int(1)
self.NSVEp_extra_parameters['tracers0_smoothness'] = int(1)
#self.extra_parameters = {}
- #for key in ['NSVE', 'NSVE_no_output', 'NSVEparticles', 'NSVEparticles_no_output', 'NSVEparticlesP2P']:
+ #for key in ['NSVE', 'NSVE_no_output', 'NSVEparticles', 'NSVEparticles_no_output', 'NSVEcomplex_particles']:
# self.extra_parameters[key] = {}
- #for key in ['NSVEparticles', 'NSVEparticles_no_output', 'NSVEparticlesP2P']:
+ #for key in ['NSVEparticles', 'NSVEparticles_no_output', 'NSVEcomplex_particles']:
# self.extra_parameters[key].update(self.NSVEp_extra_parameters)
return None
def get_kspace(self):
@@ -387,7 +387,7 @@ 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', 'NSVEparticlesP2P', 'NSVEparticles']:
+ if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles']:
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)
@@ -443,7 +443,7 @@ class DNS(_code):
for val in pbase_shape[1:]:
number_of_particles *= val
ncomponents = 3
- if self.dns_type in ['NSVEparticlesP2P']:
+ if self.dns_type in ['NSVEcomplex_particles']:
ncomponents = 6
with h5py.File(self.get_checkpoint_0_fname(), 'a') as ofile:
s = 0
@@ -622,8 +622,8 @@ class DNS(_code):
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
parser_NSVEp2p = subparsers.add_parser(
- 'NSVEparticlesP2P',
- help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
+ 'NSVEcomplex_particles',
+ help = 'plain Navier-Stokes vorticity formulation, with oriented active particles')
for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p']:
eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
@@ -668,7 +668,7 @@ class DNS(_code):
self.dns_type = opt.DNS_class
self.name = self.dns_type + '-' + self.fluid_precision + '-v' + bfps.__version__
# merge parameters if needed
- if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']:
+ if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']:
for k in self.NSVEp_extra_parameters.keys():
self.parameters[k] = self.NSVEp_extra_parameters[k]
if type(extra_parameters) != type(None):
@@ -728,7 +728,7 @@ 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 ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']:
+ if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']:
rhs_size = self.parameters['tracers0_integration_steps']
if type(opt.tracers0_integration_steps) != type(None):
rhs_size = opt.tracers0_integration_steps
@@ -932,9 +932,9 @@ class DNS(_code):
particle_file.create_group('tracers0/position')
particle_file.create_group('tracers0/velocity')
particle_file.create_group('tracers0/acceleration')
- if self.dns_type in ['NSVEparticlesP2P']:
+ if self.dns_type in ['NSVEcomplex_particles']:
particle_file.create_group('tracers0/orientation')
- particle_file.create_group('tracers0/vorticity')
+ particle_file.create_group('tracers0/velocity_gradient')
return None
def launch_jobs(
self,
@@ -994,7 +994,7 @@ class DNS(_code):
# particle_initial_condition[..., 2] += onedarray[None, :, None, None]
self.write_par(
particle_ic = None)
- if self.dns_type in ['NSVEparticles', 'NSVEparticlesP2P', 'NSVEparticles_no_output']:
+ if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']:
self.generate_particle_data(opt = opt)
self.run(
nb_processes = opt.nb_processes,
diff --git a/bfps/cpp/full_code/NSVEcomplex_particles.cpp b/bfps/cpp/full_code/NSVEcomplex_particles.cpp
index a6b7082e..93d0edc1 100644
--- a/bfps/cpp/full_code/NSVEcomplex_particles.cpp
+++ b/bfps/cpp/full_code/NSVEcomplex_particles.cpp
@@ -1,13 +1,13 @@
#include <string>
#include <cmath>
-#include "NSVEparticlesP2P.hpp"
+#include "NSVEcomplex_particles.hpp"
#include "scope_timer.hpp"
#include "particles/particles_sampling.hpp"
#include "particles/p2p_computer.hpp"
#include "particles/particles_inner_computer.hpp"
template <typename rnumber>
-int NSVEparticlesP2P<rnumber>::initialize(void)
+int NSVEcomplex_particles<rnumber>::initialize(void)
{
this->NSVE<rnumber>::initialize();
@@ -50,14 +50,18 @@ int NSVEparticlesP2P<rnumber>::initialize(void)
(this->simname + "_particles.h5"),
"tracers0",
"position/0");
- // TODO: remove the following testing initial condition, and use a proper
- // way to initialize with 0 (i.e. generate a 0 field as the initial condition).
- //*this->fs->cvorticity = 0.0;
+
+
+ /// allocate grad vel field
+ this->nabla_u = new field<rnumber, FFTW, THREExTHREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ this->fs->cvorticity->fftw_plan_rigor);
return EXIT_SUCCESS;
}
template <typename rnumber>
-int NSVEparticlesP2P<rnumber>::step(void)
+int NSVEcomplex_particles<rnumber>::step(void)
{
this->fs->compute_velocity(this->fs->cvorticity);
this->fs->cvelocity->ift();
@@ -74,7 +78,7 @@ int NSVEparticlesP2P<rnumber>::step(void)
}
template <typename rnumber>
-int NSVEparticlesP2P<rnumber>::write_checkpoint(void)
+int NSVEcomplex_particles<rnumber>::write_checkpoint(void)
{
this->NSVE<rnumber>::write_checkpoint();
this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
@@ -90,8 +94,9 @@ int NSVEparticlesP2P<rnumber>::write_checkpoint(void)
}
template <typename rnumber>
-int NSVEparticlesP2P<rnumber>::finalize(void)
+int NSVEcomplex_particles<rnumber>::finalize(void)
{
+ delete this->nabla_u;
delete this->particles_output_writer_mpi;
delete this->particles_sample_writer_mpi;
this->NSVE<rnumber>::finalize();
@@ -102,7 +107,7 @@ int NSVEparticlesP2P<rnumber>::finalize(void)
*/
template <typename rnumber>
-int NSVEparticlesP2P<rnumber>::do_stats()
+int NSVEcomplex_particles<rnumber>::do_stats()
{
/// perform fluid stats
this->NSVE<rnumber>::do_stats();
@@ -115,6 +120,7 @@ int NSVEparticlesP2P<rnumber>::do_stats()
/// allocate temporary data array
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>(
@@ -147,15 +153,20 @@ int NSVEparticlesP2P<rnumber>::do_stats()
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
- /// sample vorticity
- *this->tmp_vec_field = this->fs->cvorticity->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>(
+ /// sample velocity gradient
+ /// fs->cvelocity should contain the velocity in Fourier space
+ this->fs->compute_velocity(this->fs->cvorticity);
+ compute_gradient(
+ this->fs->kk,
+ this->fs->cvelocity,
+ this->nabla_u);
+ this->nabla_u->ift();
+ this->ps->sample_compute_field(*this->nabla_u, pdata2.get());
+ this->particles_sample_writer_mpi->template save_dataset<9>(
"tracers0",
- "vorticity",
+ "velocity_gradient",
pdata0.get(),
- &pdata1,
+ &pdata2,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
@@ -176,6 +187,6 @@ int NSVEparticlesP2P<rnumber>::do_stats()
return EXIT_SUCCESS;
}
-template class NSVEparticlesP2P<float>;
-template class NSVEparticlesP2P<double>;
+template class NSVEcomplex_particles<float>;
+template class NSVEcomplex_particles<double>;
diff --git a/bfps/cpp/full_code/NSVEcomplex_particles.hpp b/bfps/cpp/full_code/NSVEcomplex_particles.hpp
index b74169f5..2015ec5b 100644
--- a/bfps/cpp/full_code/NSVEcomplex_particles.hpp
+++ b/bfps/cpp/full_code/NSVEcomplex_particles.hpp
@@ -24,8 +24,8 @@
-#ifndef NSVEPARTICLESP2P_HPP
-#define NSVEPARTICLESP2P_HPP
+#ifndef NSVECOMPLEX_PARTICLES_HPP
+#define NSVECOMPLEX_PARTICLES_HPP
@@ -37,15 +37,18 @@
#include "particles/particles_output_hdf5.hpp"
#include "particles/particles_sampling.hpp"
-/** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
+/** \brief Navier-Stokes solver that includes complex particles.
*
* Child of Navier Stokes vorticity equation solver, this class calls all the
- * methods from `NSVE`, and in addition integrates simple Lagrangian tracers
+ * methods from `NSVE`, and in addition integrates `complex particles`
* in the resulting velocity field.
+ * By `complex particles` we mean neutrally buoyant, very small particles,
+ * which have an orientation and actively swim in that direction, and they may
+ * also interact with each other, trying to reorient to a common orientation.
*/
template <typename rnumber>
-class NSVEparticlesP2P: public NSVE<rnumber>
+class NSVEcomplex_particles: public NSVE<rnumber>
{
public:
@@ -67,16 +70,18 @@ class NSVEparticlesP2P: public NSVE<rnumber>
// TODO P2P use a reader with particle data
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;
+ // field for sampling velocity gradient
+ field<rnumber, FFTW, THREExTHREE> *nabla_u;
- NSVEparticlesP2P(
+ NSVEcomplex_particles(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
NSVE<rnumber>(
COMMUNICATOR,
simulation_name),
cutoff(10), inner_v0(1), enable_p2p(true), enable_inner(true), enable_vorticity_omega(true){}
- ~NSVEparticlesP2P(){}
+ ~NSVEcomplex_particles(){}
int initialize(void);
int step(void);
@@ -87,5 +92,5 @@ class NSVEparticlesP2P: public NSVE<rnumber>
int do_stats(void);
};
-#endif//NSVEPARTICLESP2P_HPP
+#endif//NSVECOMPLEX_PARTICLES_HPP
diff --git a/bfps/test/test_particles.py b/bfps/test/test_particles.py
index 40c428bd..6c12fac1 100644
--- a/bfps/test/test_particles.py
+++ b/bfps/test/test_particles.py
@@ -20,7 +20,7 @@ def main():
if sys.argv[2] == 'on':
c = DNS()
c.launch(
- ['NSVEparticlesP2P',
+ ['NSVEcomplex_particles',
'-n', '32',
'--src-simname', 'B32p1e4',
'--src-wd', bfps.lib_dir + '/test',
diff --git a/setup.py b/setup.py
index a160ce4b..ff73b945 100644
--- a/setup.py
+++ b/setup.py
@@ -88,7 +88,7 @@ print('This is bfps version ' + VERSION)
### lists of files and MANIFEST.in
-src_file_list = ['full_code/NSVEparticlesP2P',
+src_file_list = ['full_code/NSVEcomplex_particles',
'full_code/joint_acc_vel_stats',
'full_code/test',
'full_code/filter_test',
--
GitLab