From 32e9dbfb6d55958ac8dae05d270141a7637d315c Mon Sep 17 00:00:00 2001
From: Cristian C Lalescu <Cristian.Lalescu@mpcdf.mpg.de>
Date: Thu, 17 Dec 2020 15:08:53 +0100
Subject: [PATCH] adds simple Kraichan scalar model for profiling
---
TurTLE/test/profile_scalar.py | 33 +-
TurTLE/test/profiler/Kraichnan_scalar_v1.cpp | 311 +++++++++++++++++++
TurTLE/test/profiler/Kraichnan_scalar_v1.hpp | 80 +++++
setup.py | 2 +
4 files changed, 411 insertions(+), 15 deletions(-)
create mode 100644 TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
create mode 100644 TurTLE/test/profiler/Kraichnan_scalar_v1.hpp
diff --git a/TurTLE/test/profile_scalar.py b/TurTLE/test/profile_scalar.py
index ce7c6eda..135ca33e 100644
--- a/TurTLE/test/profile_scalar.py
+++ b/TurTLE/test/profile_scalar.py
@@ -18,7 +18,9 @@ cpp_location = os.path.join(os.path.join(
class ADNS(TurTLE.DNS):
def __init__(
self,
+ dns_type = 'scalar_evolution',
**kwargs):
+ self.dns_type = dns_type
TurTLE.DNS.__init__(self, **kwargs)
return None
def write_src(
@@ -56,13 +58,12 @@ class ADNS(TurTLE.DNS):
['#include ' + hh
for hh in self.include_list])
self.name = 'scalar_evolution'
- self.dns_type = 'scalar_evolution'
self.definitions += open(
os.path.join(
- cpp_location, 'scalar_evolution.hpp'), 'r').read()
+ cpp_location, self.dns_type + '.hpp'), 'r').read()
self.definitions += open(
os.path.join(
- cpp_location, 'scalar_evolution.cpp'), 'r').read()
+ cpp_location, self.dns_type + '.cpp'), 'r').read()
with open(self.name + '.cpp', 'w') as outfile:
outfile.write(self.version_message + '\n\n')
outfile.write(self.includes + '\n\n')
@@ -84,12 +85,14 @@ class ADNS(TurTLE.DNS):
self.parameters['dt'] = float(0.0001)
self.parameters['histogram_bins'] = int(64)
self.parameters['max_value_estimate'] = float(1)
+ self.parameters['field_random_seed'] = int(1)
self.parameters['nu'] = float(0.1)
- self.parameters['mu'] = float(0.2)
- self.parameters['injection_rate'] = float(0.4)
- self.parameters['fk0'] = float(2.0)
- self.parameters['fk1'] = float(4.0)
+ if self.dns_type == 'scalar_evolution':
+ self.parameters['mu'] = float(0.2)
+ self.parameters['injection_rate'] = float(0.4)
+ self.parameters['fk0'] = float(2.0)
+ self.parameters['fk1'] = float(4.0)
self.parameters['spectrum_dissipation'] = float(0.027)
self.parameters['spectrum_Lint'] = float(1.3)
@@ -108,7 +111,6 @@ class ADNS(TurTLE.DNS):
opt = parser.parse_args(args)
self.simname=opt.simname
self.set_precision(opt.precision)
- self.dns_type = 'scalar_evolution'
self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
if ((self.parameters['niter_todo'] % self.parameters['niter_out']) != 0):
self.parameters['niter_out'] = self.parameters['niter_todo']
@@ -126,12 +128,13 @@ class ADNS(TurTLE.DNS):
opt.ny = opt.n
if type(opt.nz) == type(None):
opt.nz = opt.n
- if type(opt.fk0) == type(None):
- opt.fk0 = self.parameters['fk0']
- if type(opt.fk1) == type(None):
- opt.fk1 = self.parameters['fk1']
- if type(opt.injection_rate) == type(None):
- opt.injection_rate = self.parameters['injection_rate']
+ if self.dns_type == 'scalar_evolution':
+ if type(opt.fk0) == type(None):
+ opt.fk0 = self.parameters['fk0']
+ if type(opt.fk1) == type(None):
+ opt.fk1 = self.parameters['fk1']
+ if type(opt.injection_rate) == type(None):
+ opt.injection_rate = self.parameters['injection_rate']
if type(opt.dealias_type) == type(None):
opt.dealias_type = self.parameters['dealias_type']
if (opt.nx > opt.n or
@@ -213,7 +216,7 @@ class ADNS(TurTLE.DNS):
return None
def main():
- bla = ADNS()
+ bla = ADNS('Kraichnan_scalar_v1')
bla.launch(sys.argv[1:])
return None
diff --git a/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp b/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
new file mode 100644
index 00000000..e180ca95
--- /dev/null
+++ b/TurTLE/test/profiler/Kraichnan_scalar_v1.cpp
@@ -0,0 +1,311 @@
+/******************************************************************************
+* *
+* 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_scalar_v1.hpp"
+#include "scope_timer.hpp"
+#include "fftw_tools.hpp"
+
+
+template <typename rnumber>
+int Kraichnan_scalar_v1<rnumber>::initialize(void)
+{
+ TIMEZONE("Kraichnan_scalar_v1::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->scalar = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+
+ this->tscal0 = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+
+ this->ux = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+ this->uy = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+ this->uz = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+
+ this->kk = new kspace<FFTW, SMOOTH>(
+ this->scalar->clayout,
+ this->dkx, this->dky, this->dkz);
+
+ make_gaussian_random_field<rnumber, FFTW, ONE, SMOOTH>(
+ this->kk,
+ this->scalar,
+ this->random_seed,
+ this->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_large_scale_const,
+ this->spectrum_small_scale_const);
+
+ this->scalar->symmetrize();
+ this->kk->template dealias<rnumber, ONE>(this->scalar->get_cdata());
+
+ if (this->myrank == 0 && this->iteration == 0)
+ this->kk->store(stat_file);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Kraichnan_scalar_v1<rnumber>::step(void)
+{
+ TIMEZONE("Kraichnan_scalar_v1::step");
+ this->iteration++;
+
+ // store scalar
+ *(this->tscal0) = *(this->scalar);
+
+ // generate random velocity field
+ make_gaussian_random_field<rnumber, FFTW, ONE, SMOOTH>(
+ this->kk,
+ this->ux,
+ 3*(this->random_seed+this->iteration),
+ this->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_large_scale_const,
+ this->spectrum_small_scale_const,
+ double(3)/2); // normalization required because we will impose incompressibility
+ make_gaussian_random_field<rnumber, FFTW, ONE, SMOOTH>(
+ this->kk,
+ this->uy,
+ 3*(this->random_seed+this->iteration)+1,
+ this->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_large_scale_const,
+ this->spectrum_small_scale_const,
+ double(3)/2);
+ make_gaussian_random_field<rnumber, FFTW, ONE, SMOOTH>(
+ this->kk,
+ this->uz,
+ 3*(this->random_seed+this->iteration)+2,
+ this->spectrum_dissipation,
+ this->spectrum_Lint,
+ this->spectrum_etaK,
+ this->spectrum_large_scale_const,
+ this->spectrum_small_scale_const,
+ double(3)/2);
+
+ // ensure velocity field is solenoidal
+ {
+ TIMEZONE("Kraichnan_scalar_v1::divfree");
+ this->kk->CLOOP_K2(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ if (k2 > 0)
+ {
+ const typename fftw_interface<rnumber>::complex tval = {
+ rnumber((this->kk->kx[xindex]*this->ux->cval(cindex, 0) +
+ this->kk->ky[yindex]*this->uy->cval(cindex, 0) +
+ this->kk->kz[zindex]*this->uz->cval(cindex, 0) ) / k2),
+ rnumber((this->kk->kx[xindex]*this->ux->cval(cindex, 1) +
+ this->kk->ky[yindex]*this->uy->cval(cindex, 1) +
+ this->kk->kz[zindex]*this->uz->cval(cindex, 1) ) / k2)};
+ for (unsigned int imag_part=0; imag_part<2; imag_part++)
+ {
+ this->ux->cval(cindex, imag_part) -= tval[imag_part]*this->kk->kx[xindex];
+ this->uy->cval(cindex, imag_part) -= tval[imag_part]*this->kk->ky[yindex];
+ this->uz->cval(cindex, imag_part) -= tval[imag_part]*this->kk->kz[zindex];
+ }
+ }
+ });
+ }
+
+ // take fields to real space
+ this->ux->ift();
+ this->uy->ift();
+ this->uz->ift();
+ this->tscal0->ift();
+
+ // compute nonlinear term
+ this->tscal0->RLOOP (
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ this->ux->rval(rindex) *= this->tscal0->rval(rindex) / this->scalar->npoints;
+ this->uy->rval(rindex) *= this->tscal0->rval(rindex) / this->scalar->npoints;
+ this->uz->rval(rindex) *= this->tscal0->rval(rindex) / this->scalar->npoints;
+ }
+ );
+
+ // take fields back to Fourier space
+ this->ux->dft();
+ this->uy->dft();
+ this->uz->dft();
+
+ this->kk->template dealias<rnumber, ONE>(this->ux->get_cdata());
+ this->kk->template dealias<rnumber, ONE>(this->uy->get_cdata());
+ this->kk->template dealias<rnumber, ONE>(this->uz->get_cdata());
+
+ // Euler step
+ this->kk->CLOOP_K2(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ if (k2 <= this->kk->kM2)
+ {
+ this->scalar->cval(cindex,0) -= this->dt*(
+ this->nu*k2 * this->scalar->cval(cindex, 0) + // linear term
+ this->kk->kx[xindex] * this->ux->cval(cindex, 1) +
+ this->kk->ky[yindex] * this->uy->cval(cindex, 1) +
+ this->kk->kz[zindex] * this->uz->cval(cindex, 1)
+ );
+ this->scalar->cval(cindex,0) += this->dt*(
+ -this->nu*k2 * this->scalar->cval(cindex, 0) + // linear term
+ this->kk->kx[xindex] * this->ux->cval(cindex, 0) +
+ this->kk->ky[yindex] * this->uy->cval(cindex, 0) +
+ this->kk->kz[zindex] * this->uz->cval(cindex, 0)
+ );
+ }
+ else
+ std::fill_n((rnumber*)(this->scalar->get_cdata()+cindex), 2, 0.0);
+ });
+
+ // symmetrize
+ this->scalar->symmetrize();
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Kraichnan_scalar_v1<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("Kraichnan_scalar_v1::write_checkpoint");
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Kraichnan_scalar_v1<rnumber>::finalize(void)
+{
+ TIMEZONE("Kraichnan_scalar_v1::finalize");
+ if (this->myrank == 0)
+ H5Fclose(this->stat_file);
+ delete this->tscal0;
+ delete this->ux;
+ delete this->uy;
+ delete this->uz;
+ delete this->scalar;
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Kraichnan_scalar_v1<rnumber>::do_stats()
+{
+ TIMEZONE("Kraichnan_scalar_v1::do_stats");
+ if (!(this->iteration % this->niter_stat == 0))
+ return EXIT_SUCCESS;
+ hid_t stat_group;
+ if (this->myrank == 0)
+ stat_group = H5Gopen(
+ this->stat_file,
+ "statistics",
+ H5P_DEFAULT);
+ else
+ stat_group = 0;
+
+ *(this->tscal0) = *(this->scalar);
+
+ this->tscal0->compute_stats(
+ this->kk,
+ stat_group,
+ "scalar",
+ this->iteration / niter_stat,
+ max_value_estimate);
+
+ if (this->myrank == 0)
+ H5Gclose(stat_group);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Kraichnan_scalar_v1<rnumber>::read_parameters(void)
+{
+ TIMEZONE("Kraichnan_scalar_v1::read_parameters");
+ this->direct_numerical_simulation::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ this->nu = hdf5_tools::read_value<double>(parameter_file, "parameters/nu");
+ this->dt = hdf5_tools::read_value<double>(parameter_file, "parameters/dt");
+ this->histogram_bins = hdf5_tools::read_value<int>(parameter_file, "parameters/histogram_bins");
+ this->max_value_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_value_estimate");
+ this->fftw_plan_rigor = hdf5_tools::read_string(parameter_file, "parameters/fftw_plan_rigor");
+
+ this->spectrum_dissipation = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_dissipation");
+ this->spectrum_Lint = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_Lint");
+ this->spectrum_etaK = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_etaK");
+ this->spectrum_large_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_large_scale_const");
+ this->spectrum_small_scale_const = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_small_scale_const");
+ this->random_seed = hdf5_tools::read_value<int>(parameter_file, "/parameters/field_random_seed");
+ H5Fclose(parameter_file);
+ return EXIT_SUCCESS;
+}
+
+template class Kraichnan_scalar_v1<float>;
+template class Kraichnan_scalar_v1<double>;
+
diff --git a/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp b/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp
new file mode 100644
index 00000000..e57735cd
--- /dev/null
+++ b/TurTLE/test/profiler/Kraichnan_scalar_v1.hpp
@@ -0,0 +1,80 @@
+/**********************************************************************
+* *
+* 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 SCALAR_EVOLUTION_HPP
+#define SCALAR_EVOLUTION_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "full_code/direct_numerical_simulation.hpp"
+
+template <typename rnumber>
+class Kraichnan_scalar_v1: public direct_numerical_simulation
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ double dt;
+ int histogram_bins;
+ double max_value_estimate;
+ double nu;
+ std::string fftw_plan_rigor;
+
+ int random_seed;
+ double spectrum_dissipation;
+ double spectrum_Lint;
+ double spectrum_etaK;
+ double spectrum_large_scale_const;
+ double spectrum_small_scale_const;
+
+ /* other stuff */
+ field<rnumber, FFTW, ONE> *scalar;
+ field<rnumber, FFTW, ONE> *tscal0;
+ field<rnumber, FFTW, ONE> *ux, *uy, *uz;
+ kspace<FFTW, SMOOTH> *kk;
+
+
+ Kraichnan_scalar_v1(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ direct_numerical_simulation(
+ COMMUNICATOR,
+ simulation_name){}
+ ~Kraichnan_scalar_v1(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ virtual int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//SCALAR_EVOLUTION_HPP
+
diff --git a/setup.py b/setup.py
index 21dc139f..a4f43074 100644
--- a/setup.py
+++ b/setup.py
@@ -82,6 +82,8 @@ setup(
'test/particle_set/NSVEparticle_set.cpp',
'test/profiler/scalar_evolution.hpp',
'test/profiler/scalar_evolution.cpp',
+ 'test/profiler/Kraichnan_scalar_v1.hpp',
+ 'test/profiler/Kraichnan_scalar_v1.cpp',
]},
entry_points = {
'console_scripts': [
--
GitLab