diff --git a/CMakeLists.txt b/CMakeLists.txt index cf9556bdf30cf44a77cd0846b4bc006a60784ebf..ad41381e3495549cba520baa5a5f6badf647d47d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py index 21ef90ededf297c557b9395eff86d7b5e7c5e5e8..fb6369b2e4ce676936a5aeae4b10b5e0a151ea2a 100644 --- a/TurTLE/DNS.py +++ b/TurTLE/DNS.py @@ -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( diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp index b8b8158968df0239717f0836d416c49a62f63901..5005000b8374aef4f743322f9d3bf6647dee89d1 100644 --- a/cpp/full_code/Gauss_field_test.cpp +++ b/cpp/full_code/Gauss_field_test.cpp @@ -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; } diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp index 3854bb1f48f8b3a724a67d77fe370372901e299a..d14c3cfb5bde59f5d1ffa80405ce0d700d25aa49 100644 --- a/cpp/full_code/Gauss_field_test.hpp +++ b/cpp/full_code/Gauss_field_test.hpp @@ -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 { diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp new file mode 100644 index 0000000000000000000000000000000000000000..a0f81f52ad1fed439ed246d021eae681db2e20b3 --- /dev/null +++ b/cpp/full_code/kraichnan_field.cpp @@ -0,0 +1,263 @@ +/****************************************************************************** +* * +* 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>; + diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp new file mode 100644 index 0000000000000000000000000000000000000000..93b044c5249f2c53efb113eab41e7421161fe56b --- /dev/null +++ b/cpp/full_code/kraichnan_field.hpp @@ -0,0 +1,88 @@ +/********************************************************************** +* * +* 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 +