Commit 3655957b authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

[broken] sample velocity gradient for complex particles

parent 0bb56546
Pipeline #25516 passed with stage
in 10 minutes and 44 seconds
......@@ -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,
......
#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>;
......@@ -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
......@@ -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',
......
......@@ -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',
......
Supports Markdown
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