Skip to content
Snippets Groups Projects
Commit 32e9dbfb authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

adds simple Kraichan scalar model for profiling

parent 2e575267
Branches
Tags
No related merge requests found
......@@ -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
......
/******************************************************************************
* *
* 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>;
/**********************************************************************
* *
* 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
......@@ -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': [
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment