Commit 90ab0dd2 authored by Lukas Bentkamp's avatar Lukas Bentkamp
Browse files

kraichnan_field

parent b22d29b2
......@@ -188,6 +188,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/static_field.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/kraichnan_field.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.cpp
......@@ -231,6 +232,7 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/static_field.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/kraichnan_field.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.hpp
......
......@@ -440,7 +440,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', 'NSVEcomplex_particles', 'NSVEparticles', 'static_field']:
if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_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)
......@@ -647,6 +647,10 @@ class DNS(_code):
'static_field',
help = 'static field with basic fluid tracers')
parser_kraichnan_field = subparsers.add_parser(
'kraichnan_field',
help = 'Kraichnan field with basic fluid tracers')
parser_NSVEp2 = subparsers.add_parser(
'NSVEparticles',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
......@@ -659,7 +663,7 @@ class DNS(_code):
'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', 'static_field']:
for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra', 'static_field', 'kraichnan_field']:
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))
......@@ -702,7 +706,7 @@ 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']:
if self.dns_type in ['NSVEparticles', 'NSVEcomplex_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):
......@@ -767,7 +771,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 ['static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
if self.dns_type in ['kraichnan_field', 'static_field', '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
......@@ -1043,7 +1047,7 @@ class DNS(_code):
f['vorticity/complex/{0}'.format(0)] = data
f.close()
# now take care of particles' initial condition
if self.dns_type in ['static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
if self.dns_type in ['kraichnan_field', 'static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
self.generate_particle_data(opt = opt)
return None
def launch_jobs(
......
......@@ -39,12 +39,15 @@ template <typename rnumber,
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed = 0)
const int rseed,
const double slope,
const double k_cutoff,
const double coefficient)
{
TIMEZONE("make_gaussian_random_field");
// initialize a separate random number generator for each thread
std::vector<std::mt19937_64> rgen;
std::normal_distribution<rnumber> rdist;
std::uniform_real_distribution<rnumber> rdist;
rgen.resize(omp_get_max_threads());
// seed random number generators such that no seed is ever repeated if we change the value of rseed.
// basically use a multi-dimensional array indexing technique to come up with actual seed.
......@@ -71,27 +74,33 @@ int make_gaussian_random_field(
switch(fc)
{
case ONE:
output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) / pow(k2, 17./12);
output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) / pow(k2, 17./12);
{
double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]);
output_field->cval(cindex,0) = coefficient * cos(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff);
output_field->cval(cindex,1) = coefficient * sin(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff);
break;
}
case THREE:
for (int cc = 0; cc<3; cc++)
{
output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) / pow(k2, 17./12);
output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) / pow(k2, 17./12);
double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]);
output_field->cval(cindex,cc,0) = coefficient * cos(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff);
output_field->cval(cindex,cc,1) = coefficient * sin(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff);
}
break;
case THREExTHREE:
for (int cc = 0; cc<3; cc++)
for (int ccc = 0; ccc<3; ccc++)
{
output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) / pow(k2, 17./12);
output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) / pow(k2, 17./12);
double temp_phase = 2*M_PI*rdist(rgen[omp_get_thread_num()]);
output_field->cval(cindex,cc,ccc,0) = coefficient * cos(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff);
output_field->cval(cindex,cc,ccc,1) = coefficient * sin(temp_phase) * pow(k2, slope/2.) * exp(- k2/k_cutoff);
}
break;
}
});
output_field->symmetrize();
return EXIT_SUCCESS;
}
......
......@@ -39,6 +39,20 @@
*
*/
template <typename rnumber,
field_backend be,
field_components fc,
kspace_dealias_type dt>
int make_gaussian_random_field(
kspace<be, dt> *kk,
field<rnumber, be, fc> *output_field,
const int rseed = 0,
const double slope = -5./3.,
const double k_cutoff = 10.,
const double coefficient = 1.);
template <typename rnumber>
class Gauss_field_test: public test
{
......
/******************************************************************************
* *
* 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 "kraichnan_field.hpp"
#include "scope_timer.hpp"
#include "fftw_tools.hpp"
#include "Gauss_field_test.hpp"
template <typename rnumber>
int kraichnan_field<rnumber>::initialize(void)
{
TIMEZONE("kraichnan_file::initialize");
this->read_iteration();
this->read_parameters();
if (this->myrank == 0)
{
// set caching parameters
hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
variable_used_only_in_assert(cache_err);
DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
this->stat_file = H5Fopen(
(this->simname + ".h5").c_str(),
H5F_ACC_RDWR,
fapl);
}
int data_file_problem;
if (this->myrank == 0)
data_file_problem = this->grow_file_datasets();
MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
if (data_file_problem > 0)
{
std::cerr <<
data_file_problem <<
" problems growing file datasets.\ntrying to exit now." <<
std::endl;
return EXIT_FAILURE;
}
/*
this->vorticity = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->vorticity->real_space_representation = false;
*/
this->velocity = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->velocity->real_space_representation = false;
//read vorticity field
/*this->vorticity->io(
this->get_current_fname(),
"vorticity",
this->iteration,
true);*/
this->kk = new kspace<FFTW, SMOOTH>(
this->velocity->clayout, this->dkx, this->dky, this->dkz);
// IS THE VELOCITY LAYOUT THE SAME AS THE VORTICITY ONE?
/*
// compute the velocity field and store
invert_curl(
this->kk,
this->vorticity,
velocity);
*/
// transform velocity and vorticity fields to real space
//this->velocity->ift();
//this->vorticity->ift();
// initialize particles
this->ps = particles_system_builder(
this->velocity, // (field object)
this->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->get_current_fname(), // particles input filename
std::string("/tracers0/state/") + std::to_string(this->iteration), // dataset name for initial input
std::string("/tracers0/rhs/") + std::to_string(this->iteration), // dataset name for initial input
tracers0_neighbours, // parameter (interpolation no neighbours)
tracers0_smoothness, // parameter
this->comm,
this->iteration+1);
// initialize output objects
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());
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 kraichnan_field<rnumber>::step(void)
{
TIMEZONE("kraichnan_file::step");
make_gaussian_random_field(this->kk, this->velocity, this->iteration, -5./3., 10., 0.01);
this->velocity->ift();
//this->kk->template project_divfree<rnumber>(this->velocity->get_cdata(), true);
// What does the ->template do?
// What is cdata and rdata?
// ADJUST TOTAL ENERGY, SLOPE AND CUTOFF
DEBUG_MSG("Iteration: %d\n",
this->iteration);
DEBUG_MSG("Real space: %d\n",
this->velocity->real_space_representation);
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(0,0));
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(0,1));
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(0,2));
DEBUG_MSG("Field entries: %g\n",
this->velocity->rval(15,2));
this->ps->completeLoop(sqrt(this->dt));
this->iteration++;
return EXIT_SUCCESS;
}
template <typename rnumber>
int kraichnan_field<rnumber>::write_checkpoint(void)
{
TIMEZONE("kraichnan_file::write_checkpoint");
this->particles_output_writer_mpi->open_file(this->get_current_fname());
this->particles_output_writer_mpi->template save<3>(
this->ps->getParticlesState(),
this->ps->getParticlesRhs(),
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->iteration);
this->particles_output_writer_mpi->close_file();
return EXIT_SUCCESS;
}
template <typename rnumber>
int kraichnan_field<rnumber>::finalize(void)
{
TIMEZONE("kraichnan_field::finalize");
if (this->myrank == 0)
H5Fclose(this->stat_file);
// good practice rule number n+1: always delete in inverse order of allocation
delete this->ps.release();
delete this->particles_output_writer_mpi;
delete this->particles_sample_writer_mpi;
delete this->kk;
delete this->velocity;
return EXIT_SUCCESS;
}
/** \brief Compute statistics.
*
*/
template <typename rnumber>
int kraichnan_field<rnumber>::do_stats()
{
TIMEZONE("kraichnan_field::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
std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
/// copy position data
/// sample position
std::copy(this->ps->getParticlesState(),
this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
pdata.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"position",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
/// sample velocity
std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
this->ps->sample_compute_field(*this->velocity, pdata.get());
this->particles_sample_writer_mpi->template save_dataset<3>(
"tracers0",
"velocity",
this->ps->getParticlesState(),
&pdata,
this->ps->getParticlesIndexes(),
this->ps->getLocalNbParticles(),
this->ps->get_step_idx()-1);
// deallocate temporary data array
delete[] pdata.release();
return EXIT_SUCCESS;
}
template <typename rnumber>
int kraichnan_field<rnumber>::read_parameters(void)
{
TIMEZONE("kraichnan_field::read_parameters");
this->direct_numerical_simulation::read_parameters();
hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
this->fftw_plan_rigor = hdf5_tools::read_string(parameter_file, "parameters/fftw_plan_rigor");
this->dt = hdf5_tools::read_value<double>(parameter_file, "parameters/dt");
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");
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
template class kraichnan_field<float>;
template class kraichnan_field<double>;
/**********************************************************************
* *
* 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 KRAICHNAN_FIELD_HPP
#define KRAICHNAN_FIELD_HPP
#include <cstdlib>
#include "base.hpp"
#include "kspace.hpp"
#include "field.hpp"
#include "full_code/direct_numerical_simulation.hpp"
#include "particles/particles_system_builder.hpp"
#include "particles/particles_output_hdf5.hpp"
#include "particles/particles_sampling.hpp"
template <typename rnumber>
class kraichnan_field: public direct_numerical_simulation
{
public:
/* parameters that are read in read_parameters */
double dt;
std::string fftw_plan_rigor;
/* particle 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;
/* other stuff */
kspace<FFTW, SMOOTH> *kk;
field<rnumber, FFTW, THREE> *velocity;
/* other stuff */
std::unique_ptr<abstract_particles_system<long long int, double>> ps;
particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
kraichnan_field(
const MPI_Comm COMMUNICATOR,
const std::string &simulation_name):
direct_numerical_simulation(
COMMUNICATOR,
simulation_name){}
~kraichnan_field(){}
int initialize(void);
int step(void);
int finalize(void);
virtual int read_parameters(void);
int write_checkpoint(void);
int do_stats(void);
};
#endif//KRAICHNAN_FIELD_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