Commit 0e6891a2 authored by Chichi Lalescu's avatar Chichi Lalescu
Browse files

addinitial version of pressure Hessian sampler

parent 3c3a23a0
Pipeline #30565 passed with stage
in 14 minutes and 42 seconds
......@@ -128,11 +128,11 @@ 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', 'NSVEcomplex_particles']:
if self.dns_type in ['NSVEparticles', 'NSVE_no_output', 'NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEp_extra_sampling']:
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')
if self.dns_type in ['NSVEparticles_no_output']:
if self.dns_type in ['NSVEparticles_no_output', 'NSVEp_extra_sampling']:
outfile.write('template <typename rnumber> int NSVEparticles<rnumber>::read_parameters(){return EXIT_SUCCESS;}\n')
outfile.write('template int NSVEparticles<float>::read_parameters();\n')
outfile.write('template int NSVEparticles<double>::read_parameters();\n\n')
......@@ -682,7 +682,11 @@ class DNS(_code):
'NSVEcomplex_particles',
help = 'plain Navier-Stokes vorticity formulation, with oriented active particles')
for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p']:
parser_NSVEp_extra = subparsers.add_parser(
'NSVEp_extra_sampling',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers, that sample velocity gradient, as well as pressure and its derivatives.')
for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra']:
eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
eval('self.job_parser_arguments({0})'.format('parser_' + parser))
eval('self.particle_parser_arguments({0})'.format('parser_' + parser))
......@@ -725,7 +729,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', 'NSVEcomplex_particles', 'NSVEparticles_no_output']:
if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
for k in self.NSVEp_extra_parameters.keys():
self.parameters[k] = self.NSVEp_extra_parameters[k]
if type(extra_parameters) != type(None):
......@@ -790,7 +794,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', 'NSVEcomplex_particles', 'NSVEparticles_no_output']:
if self.dns_type in ['NSVEparticles', 'NSVEcomplex_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
......@@ -997,6 +1001,11 @@ 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 ['NSVEp_extra_sampling']:
particle_file.create_group('tracers0/velocity_gradient')
particle_file.create_group('tracers0/pressure')
particle_file.create_group('tracers0/pressure_gradient')
particle_file.create_group('tracers0/pressure_Hessian')
return None
def launch_jobs(
self,
......@@ -1056,7 +1065,7 @@ class DNS(_code):
# particle_initial_condition[..., 2] += onedarray[None, :, None, None]
self.write_par(
particle_ic = None)
if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output']:
if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
self.generate_particle_data(opt = opt)
self.run(
nb_processes = opt.nb_processes,
......
#include "full_code/NSVEp_extra_sampling.hpp"
template <typename rnumber>
int NSVEp_extra_sampling<rnumber>::initialize(void)
{
this->NSVEparticles<rnumber>::initialize();
/// 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);
this->pressure = new field<rnumber, FFTW, ONE>(
this->nx, this->ny, this->nz,
this->comm,
this->fs->cvorticity->fftw_plan_rigor);
this->nabla_p = new field<rnumber, FFTW, THREE>(
this->nx, this->ny, this->nz,
this->comm,
this->fs->cvorticity->fftw_plan_rigor);
this->Hessian_p = 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 NSVEp_extra_sampling<rnumber>::finalize(void)
{
delete this->nabla_u;
delete this->pressure;
delete this->nabla_p;
delete this->Hessian_p;
this->NSVEparticles<rnumber>::finalize();
return EXIT_SUCCESS;
}
template <typename rnumber>
int NSVEp_extra_sampling<rnumber>::do_stats()
{
this->NSVEparticles<rnumber>::do_stats();
if (!(this->iteration % this->niter_part == 0))
return EXIT_SUCCESS;
/// 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->fs->compute_pressure(this->pressure);
compute_gradient(
this->fs->kk,
this->pressure,
this->nabla_p);
compute_gradient(
this->fs->kk,
this->nabla_p,
this->Hessian_p);
this->pressure->ift();
this->nabla_p->ift();
this->Hessian_p->ift();
// sample velocity gradient
std::unique_ptr<double[]> pdata(new double[9*this->ps->getLocalNbParticles()]);
std::fill_n(pdata.get(), 9*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->nabla_u, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<9>(
"tracers0",
"velocity_gradient",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
// sample pressure
std::fill_n(pdata.get(), this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->pressure, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<1>(
"tracers0",
"pressure",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
// sample pressure gradient
std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->nabla_p, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"pressure_gradient",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
// sample pressure gradient
std::fill_n(pdata.get(), 9*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->Hessian_p, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<9>(
"tracers0",
"pressure_Hessian",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
return EXIT_SUCCESS;
}
template class NSVEp_extra_sampling<float>;
template class NSVEp_extra_sampling<double>;
#ifndef NSVEP_EXTRA_SAMPLING_HPP
#define NSVEP_EXTRA_SAMPLING_HPP
#include <cstdlib>
#include "base.hpp"
#include "vorticity_equation.hpp"
#include "full_code/NSVEparticles.hpp"
#include "particles/particles_system_builder.hpp"
#include "particles/particles_output_hdf5.hpp"
#include "particles/particles_sampling.hpp"
/** \brief Navier-Stokes solver with tracers that sample velocity gradient
* and pressure Hessian.
*
*/
template <typename rnumber>
class NSVEp_extra_sampling: public NSVEparticles<rnumber>
{
public:
/* other stuff */
field<rnumber, FFTW, ONE> *pressure;
field<rnumber, FFTW, THREE> *nabla_p;
field<rnumber, FFTW, THREExTHREE> *nabla_u;
field<rnumber, FFTW, THREExTHREE> *Hessian_p;
NSVEp_extra_sampling(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
NSVEparticles<rnumber>(
COMMUNICATOR,
simulation_name){}
~NSVEp_extra_sampling(){}
int initialize(void);
int finalize(void);
int read_parameters(void);
int do_stats(void);
};
#endif//NSVEP_EXTRA_SAMPLING_HPP
......@@ -133,7 +133,8 @@ src_file_list = ['full_code/NSVEcomplex_particles',
'spline_n10',
'Lagrange_polys',
'scope_timer',
'full_code/NSVEparticles']
'full_code/NSVEparticles',
'full_code/NSVEp_extra_sampling']
particle_headers = [
'cpp/particles/abstract_particles_input.hpp',
......
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