diff --git a/CMakeLists.txt b/CMakeLists.txt
index 148b89dd01178eb88051a6feff89f9525434624a..a76751f64600ed3d5140ac69c80799745be011d2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -225,7 +225,10 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.cpp
- ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp)
+ ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.cpp)
+
set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/code_base.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.hpp
@@ -305,6 +308,8 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/omputils.hpp
${PROJECT_SOURCE_DIR}/cpp/shared_array.hpp
${PROJECT_SOURCE_DIR}/cpp/spline.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.hpp
)
#file(GLOB_RECURSE hpp_for_lib ${PROJECT_SOURCE_DIR}/*.hpp)
LIST(APPEND source_files ${hpp_for_lib} ${cpp_for_lib})
@@ -336,4 +341,3 @@ else()
install(CODE "execute_process(COMMAND ${CMAKE_COMMAND} -E copy ${PROJECT_SOURCE_DIR}/pc_host_info.py ${PROJECT_BINARY_DIR}/python/TurTLE/host_info.py)")
endif()
install(CODE "execute_process(COMMAND python3 ${PROJECT_SOURCE_DIR}/setup.py install --force --prefix=${CMAKE_INSTALL_PREFIX} WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/python/)")
-
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 9b490c778b5cdf97d0ef0b7cf2f0244e9ee68d45..1d5bea0ca4cb3e01b03b8562bd006d3af2c2ced3 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -2,20 +2,20 @@
# #
# Copyright 2015-2019 Max Planck Institute for Dynamics and Self-Organization #
# #
-# This file is part of TurTLE. #
+# This file is part of TurTLE. #
# #
-# TurTLE is free software: you can redistribute it and/or modify #
+# TurTLE 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. #
# #
-# TurTLE is distributed in the hope that it will be useful, #
+# TurTLE 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 TurTLE. If not, see <http://www.gnu.org/licenses/> #
+# along with TurTLE. If not, see <http://www.gnu.org/licenses/> #
# #
# Contact: Cristian.Lalescu@ds.mpg.de #
# #
@@ -657,7 +657,9 @@ class DNS(_code):
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.')
-
+ parser_NSVE_ou = subparsers.add_parser(
+ 'NSVE_ou_forcing',
+ help = 'plain Navier-Stokes vorticity formulation, with ornstein-uhlenbeck forcing')
for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra', 'static_field']:
eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
eval('self.job_parser_arguments({0})'.format('parser_' + parser))
@@ -1060,4 +1062,3 @@ class DNS(_code):
no_submit = opt.no_submit,
no_debug = opt.no_debug)
return None
-
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index dbf166d413dc2024e09e0e66c6e8b97d4c74df71..14f1e138ab3c2ce411d81990faa351557d820e94 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -275,12 +275,18 @@ class TEST(_code):
parser_test_interpolation = subparsers.add_parser(
'test_interpolation',
help = 'test velocity gradient interpolation')
+
+ parser_ornstein_uhlenbeck_test = subparsers.add_parser(
+ 'ornstein_uhlenbeck_test',
+ help = 'test ornstein uhlenbeck process')
+
for parser in ['parser_filter_test',
'parser_write_filtered_test',
'parser_field_test',
'parser_symmetrize_test',
'parser_field_output_test',
- 'parser_test_interpolation']:
+ 'parser_test_interpolation',
+ 'parser_ornstein_uhlenbeck_test']:
eval('self.simulation_parser_arguments(' + parser + ')')
eval('self.job_parser_arguments(' + parser + ')')
eval('self.parameters_to_parser_arguments(' + parser + ')')
@@ -368,6 +374,8 @@ class TEST(_code):
ofile.require_group('tracers0')
for kk in ['position', 'velocity', 'vorticity', 'velocity_gradient']:
ofile['tracers0'].require_group(kk)
+
+
self.run(
nb_processes = opt.nb_processes,
nb_threads_per_process = opt.nb_threads_per_process,
@@ -376,4 +384,3 @@ class TEST(_code):
minutes = opt.minutes % 60,
no_submit = opt.no_submit)
return None
-
diff --git a/cpp/full_code/NSVE.cpp b/cpp/full_code/NSVE.cpp
index 21788cc9f305570dfdcd066421a83bbc065120ae..909994305db5c9233e028f22aa343766f1912281 100644
--- a/cpp/full_code/NSVE.cpp
+++ b/cpp/full_code/NSVE.cpp
@@ -62,11 +62,11 @@ int NSVE<rnumber>::initialize(void)
}
this->fs = new vorticity_equation<rnumber, FFTW>(
simname.c_str(),
- nx, ny, nz,
- dkx, dky, dkz,
+ this->nx, this->ny, this->nz,
+ this->dkx, this->dky, this->dkz,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->tmp_vec_field = new field<rnumber, FFTW, THREE>(
- nx, ny, nz,
+ this->nx, this->ny, this->nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..b6af42d324a27cddedd56aefbf7c57632ce186c6
--- /dev/null
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -0,0 +1,280 @@
+#include "ornstein_uhlenbeck_process.hpp"
+#include <cmath>
+#include <cstring>
+#include <cassert>
+#include "scope_timer.hpp"
+#include <algorithm>
+#define NDEBUG
+
+
+template <class rnumber,field_backend be>
+ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
+ const char *NAME,
+ int nx,
+ int ny,
+ int nz,
+ double ou_kmin,
+ double ou_kmax,
+ double ou_energy_amplitude,
+ double DKX,
+ double DKY,
+ double DKZ,
+ unsigned FFTW_PLAN_RIGOR)
+{
+ TIMEZONE("ornstein_uhlenbeck_process::ornstein_uhlenbeck_process");
+ strncpy(this->name, NAME, 256);
+ this->name[255] = '\0';
+ this->iteration = 0;
+ this->ou_field = new field<rnumber,be,THREE>(
+ nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
+ *this->ou_field = 0.0;
+ this->ou_field->dft();
+
+ this->ou_field_vort = new field<rnumber,be,THREE>(
+ nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
+ *this->ou_field_vort = 0.0;
+ this->ou_field_vort->dft();
+
+ this->B = new field<rnumber,be,THREExTHREE>(
+ nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
+
+ this->kk = new kspace<be,SMOOTH>(
+ this->ou_field->clayout, DKX, DKY, DKZ);
+
+ this->ou_kmin_squ = pow(ou_kmin,2);
+ this->ou_kmax_squ = pow(ou_kmax,2);
+ this->ou_energy_amplitude = ou_energy_amplitude;
+ this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);
+
+ assert(this->kk->kM2 >= this->ou_kmax_squ);
+
+ gen.resize(omp_get_max_threads());
+
+ for(int thread=0;thread<omp_get_max_threads();thread++)
+ {
+ int current_seed =
+ this->ou_field->clayout->myrank*omp_get_max_threads() +
+ thread;
+ gen[thread].seed(current_seed);
+ }
+
+
+ this->initialize_B();
+
+}
+
+template <class rnumber,field_backend be>
+ornstein_uhlenbeck_process<rnumber,be>::~ornstein_uhlenbeck_process()
+{
+ TIMEZONE("ornstein_uhlenbeck_process::~ornstein_uhlenbeck_process");
+ delete this->kk;
+ delete this->ou_field;
+ delete this->ou_field_vort;
+ delete this->B;
+}
+
+template <class rnumber,field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::step(
+ double dt)
+{
+ // put in "CFL"-criterium TODO!!!
+ TIMEZONE("ornstein_uhlenbeck_process::step");
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ)
+ {
+ double g = this->gamma(sqrt(k2));
+ int tr = omp_get_thread_num();
+ double random[6] =
+ {
+ this->d(this->gen[tr]),this->d(this->gen[tr]),this->d(this->gen[tr]),
+ this->d(this->gen[tr]),this->d(this->gen[tr]),this->d(this->gen[tr])
+ };
+ for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+ this->ou_field->cval(cindex,cc,i) =
+ (1-dt*g) * this->ou_field->cval(cindex,cc,i) +
+ sqrt(dt) * (
+ this->B->cval(cindex,0,cc,0) * random[0+3*i] +
+ this->B->cval(cindex,1,cc,0) * random[1+3*i] +
+ this->B->cval(cindex,2,cc,0) * random[2+3*i] );
+
+ }
+ });
+
+ this->ou_field->symmetrize();
+ this->calc_ou_vorticity();
+
+}
+
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
+{
+ TIMEZONE("ornstein_uhlenbeck_process::initialize_B");
+ *this->B = 0.0;
+ this->kk->CLOOP(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex){
+ double ksqu = pow(this->kk->kx[xindex],2)+
+ pow(this->kk->ky[yindex],2)+
+ pow(this->kk->kz[zindex],2);
+ if (ksqu <= this->ou_kmax_squ && ksqu >= this->ou_kmin_squ)
+ {
+ double kabs = sqrt(ksqu);
+ double sigma;
+ if (ksqu == 0 || ksqu == this->ou_kmax_squ)
+ {
+ ksqu = 1; kabs = 1; sigma = 0;
+ }
+ else
+ {
+ sigma =
+ sqrt(4*this->gamma(kabs)*this->energy(kabs)
+ /(4*M_PI*ksqu));
+ }
+
+ for(int i=0;i<3;i++) {
+ for(int j=0;j<3;j++) {
+
+ if (i+j == 0)
+ this->B->cval(cindex,i,j,0) =
+ sigma/2. * (1-pow(this->kk->kx[xindex],2)/ksqu);
+
+ if (i+j == 4)
+ this->B->cval(cindex,i,j,0) =
+ sigma/2. * (1-pow(this->kk->kz[zindex],2)/ksqu);
+
+ if (i+j == 1)
+ this->B->cval(cindex,i,j,0) =
+ sigma/2. * (0-this->kk->kx[xindex]*this->kk->ky[yindex]/ksqu);
+
+ if (i+j == 3)
+ this->B->cval(cindex,i,j,0) =
+ sigma/2. * (0-this->kk->kz[zindex]*this->kk->ky[yindex]/ksqu);
+
+ if (i+j == 2) {
+
+ if(i==j)
+ this->B->cval(cindex,i,j,0) =
+ sigma/2. * (1-pow(this->kk->ky[yindex],2)/ksqu);
+
+ if(i!=j)
+ this->B->cval(cindex,i,j,0) =
+ sigma/2. * (0-this->kk->kx[xindex]*this->kk->kz[zindex]/ksqu);
+ }
+ }
+ }
+ }
+
+ });
+}
+
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
+{
+ //add some logic here TODO
+ for(int i=0; i<1000; i++)
+ {
+ this->step(2e-3);
+ }
+}
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber, be>::strip_from_field(
+ field<rnumber, be, THREE> *src)
+{
+ assert(src->real_space_representation==false);
+
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+
+ if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ){
+
+ for(int cc=0; cc < 3; cc++){
+ for(int imag=0; imag < 2; imag++){
+ src->cval(cindex,cc,imag) = 0;
+ }
+ }
+ }
+
+ }
+
+ );
+
+}
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_replace(
+ field<rnumber, be, THREE> *src, std::string uv)
+{
+ assert(src->real_space_representation==false);
+ assert((uv == "vorticity") || (uv == "velocity"));
+
+ field<rnumber, be, THREE> *field_to_replace;
+
+ if (uv == "vorticity") field_to_replace = this->ou_field_vort;
+ else field_to_replace = this->ou_field;
+
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+
+ if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ){
+
+ rnumber tmp;
+
+ for(int cc=0; cc < 3; cc++){
+ for(int imag=0; imag < 2; imag++){
+ tmp = field_to_replace->cval(cindex,cc,imag);
+ src->cval(cindex,cc,imag) = tmp;
+ }
+ }
+ }
+
+ }
+
+ );
+}
+
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_vorticity(void)
+{
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 <= this->kk->kM2)
+ {
+ this->ou_field_vort->cval(cindex,0,0) = -(this->kk->ky[yindex]*this->ou_field->cval(cindex,2,1) - this->kk->kz[zindex]*this->ou_field->cval(cindex,1,1));
+ this->ou_field_vort->cval(cindex,0,1) = (this->kk->ky[yindex]*this->ou_field->cval(cindex,2,0) - this->kk->kz[zindex]*this->ou_field->cval(cindex,1,0));
+ this->ou_field_vort->cval(cindex,1,0) = -(this->kk->kz[zindex]*this->ou_field->cval(cindex,0,1) - this->kk->kx[xindex]*this->ou_field->cval(cindex,2,1));
+ this->ou_field_vort->cval(cindex,1,1) = (this->kk->kz[zindex]*this->ou_field->cval(cindex,0,0) - this->kk->kx[xindex]*this->ou_field->cval(cindex,2,0));
+ this->ou_field_vort->cval(cindex,2,0) = -(this->kk->kx[xindex]*this->ou_field->cval(cindex,1,1) - this->kk->ky[yindex]*this->ou_field->cval(cindex,0,1));
+ this->ou_field_vort->cval(cindex,2,1) = (this->kk->kx[xindex]*this->ou_field->cval(cindex,1,0) - this->kk->ky[yindex]*this->ou_field->cval(cindex,0,0));
+ }
+ else
+ std::fill_n((rnumber*)(this->ou_field_vort->get_cdata()+3*cindex), 6, 0.0);
+ }
+ );
+ this->ou_field_vort->symmetrize();
+}
+
+template class ornstein_uhlenbeck_process<float,FFTW>;
+template class ornstein_uhlenbeck_process<double,FFTW>;
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a1aeaa98dd948ca39b2c96177341890b345c9953
--- /dev/null
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -0,0 +1,76 @@
+#ifndef ORNSTEIN_UHLENBECK_PROCESS_HPP
+#define ORNSTEIN_UHLENBECK_PROCESS_HPP
+
+#include <sys/stat.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iostream>
+#include "field.hpp"
+#include <random>
+
+extern int myrank, nprocs;
+
+template <typename rnumber, field_backend be>
+class ornstein_uhlenbeck_process{
+ public:
+ char name[256];
+
+ int iteration;
+ double ou_kmin_squ;
+ double ou_kmax_squ;
+ double ou_energy_amplitude;
+ double kolmogorov_constant = 2;
+ double epsilon;
+
+ field<rnumber,be,THREE> *ou_field;
+ field<rnumber,be,THREE> *ou_field_vort;
+ field<rnumber,be,THREExTHREE> *B;
+ kspace<be,SMOOTH> *kk;
+
+ std::vector<std::mt19937_64> gen;
+ std::normal_distribution<double> d{0,1};
+
+ ornstein_uhlenbeck_process(
+ const char *NAME,
+ int nx,
+ int ny,
+ int nz,
+ double ou_kmin,
+ double ou_kmax,
+ double ou_energy_amplitude,
+ double DKX = 1.0,
+ double DKY = 1.0,
+ double DKZ = 1.0,
+ unsigned FFTW_PLAN_RIGOR = FFTW_MEASURE);
+ ~ornstein_uhlenbeck_process(void);
+
+ void step(double dt);
+ inline double energy(double kabs)
+ {
+ return this->ou_energy_amplitude * pow(kabs,-5./3);
+ }
+
+ inline double gamma(double kabs)
+ {
+ return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
+ this->ou_energy_amplitude);
+ }
+
+
+ void initialize_B(void);
+
+ void let_converge(void);
+
+ void add_to_field_replace(
+ field<rnumber,be,THREE> *src, std::string uv);
+
+ void strip_from_field(
+ field<rnumber, be, THREE> *src);
+
+ void calc_ou_vorticity(void);
+
+
+};
+
+
+#endif
diff --git a/cpp/full_code/ou_vorticity_equation.cpp b/cpp/full_code/ou_vorticity_equation.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..37d3570d6f0c19849f5b584d4eaac5c7ed60c43f
--- /dev/null
+++ b/cpp/full_code/ou_vorticity_equation.cpp
@@ -0,0 +1,164 @@
+#include "ou_vorticity_equation.hpp"
+#include <cmath>
+#include <cstring>
+#include <cassert>
+#include "scope_timer.hpp"
+#define NDEBUG
+
+
+template <class rnumber, field_backend be>
+ou_vorticity_equation<rnumber,be>::~ou_vorticity_equation()
+{
+ delete this->ou;
+}
+
+template <class rnumber,
+ field_backend be>
+void ou_vorticity_equation<rnumber, be>::step(double dt)
+{
+ DEBUG_MSG("vorticity_equation::step\n");
+ TIMEZONE("vorticity_equation::step");
+ this->ou->add_to_field_replace(this->cvorticity, "vorticity");
+ *this->v[1] = 0.0;
+ this->omega_nonlin(0);
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 <= this->kk->kM2)
+ {
+ double factor0;
+ factor0 = exp(-this->nu * k2 * dt);
+ for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+ this->v[1]->cval(cindex,cc,i) = (
+ this->v[0]->cval(cindex,cc,i) +
+ dt*this->u->cval(cindex,cc,i))*factor0;
+ }
+ }
+ );
+
+ this->omega_nonlin(1);
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 <= this->kk->kM2)
+ {
+ double factor0, factor1;
+ factor0 = exp(-this->nu * k2 * dt/2);
+ factor1 = exp( this->nu * k2 * dt/2);
+ for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+ this->v[2]->cval(cindex, cc, i) = (
+ 3*this->v[0]->cval(cindex,cc,i)*factor0 +
+ ( this->v[1]->cval(cindex,cc,i) +
+ dt*this->u->cval(cindex,cc,i))*factor1)*0.25;
+ }
+ }
+ );
+
+ this->omega_nonlin(2);
+ // store old vorticity
+ *this->v[1] = *this->v[0];
+ this->kk->CLOOP_K2(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex,
+ double k2){
+ if (k2 <= this->kk->kM2)
+ {
+ double factor0;
+ factor0 = exp(-this->nu * k2 * dt * 0.5);
+ for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+ this->v[3]->cval(cindex,cc,i) = (
+ this->v[0]->cval(cindex,cc,i)*factor0 +
+ 2*(this->v[2]->cval(cindex,cc,i) +
+ dt*this->u->cval(cindex,cc,i)))*factor0/3;
+ }
+ }
+ );
+
+ this->kk->template force_divfree<rnumber>(this->cvorticity->get_cdata());
+ this->cvorticity->symmetrize();
+ this->iteration++;
+}
+
+
+template <class rnumber, field_backend be>
+void ou_vorticity_equation<rnumber, be>::omega_nonlin(int src)
+{
+ DEBUG_MSG("ou_vorticity_equation::omega_nonlin(%d)\n", src);
+ assert(src >= 0 && src < 3);
+ this->compute_velocity(this->v[src]);
+
+ this->add_ou_forcing_velocity();
+ // TIME STEP TODO
+
+ /* get fields from Fourier space to real space */
+ this->u->ift();
+ this->rvorticity->real_space_representation = false;
+ *this->rvorticity = this->v[src]->get_cdata();
+ this->rvorticity->ift();
+ /* compute cross product $u \times \omega$, and normalize */
+ this->u->RLOOP(
+ [&](ptrdiff_t rindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex){
+ rnumber tmp[3];
+ for (int cc=0; cc<3; cc++)
+ tmp[cc] = (this->u->rval(rindex,(cc+1)%3)*this->rvorticity->rval(rindex,(cc+2)%3) -
+ this->u->rval(rindex,(cc+2)%3)*this->rvorticity->rval(rindex,(cc+1)%3));
+ for (int cc=0; cc<3; cc++)
+ this->u->rval(rindex,cc) = tmp[cc] / this->u->npoints;
+ }
+ );
+ /* go back to Fourier space */
+ //this->clean_up_real_space(this->ru, 3);
+ this->u->dft();
+ this->kk->template dealias<rnumber, THREE>(this->u->get_cdata());
+ /* $\imath k \times Fourier(u \times \omega)$ */
+ this->kk->CLOOP(
+ [&](ptrdiff_t cindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex){
+ rnumber tmp[3][2];
+ {
+ tmp[0][0] = -(this->kk->ky[yindex]*this->u->cval(cindex,2,1) - this->kk->kz[zindex]*this->u->cval(cindex,1,1));
+ tmp[1][0] = -(this->kk->kz[zindex]*this->u->cval(cindex,0,1) - this->kk->kx[xindex]*this->u->cval(cindex,2,1));
+ tmp[2][0] = -(this->kk->kx[xindex]*this->u->cval(cindex,1,1) - this->kk->ky[yindex]*this->u->cval(cindex,0,1));
+ tmp[0][1] = (this->kk->ky[yindex]*this->u->cval(cindex,2,0) - this->kk->kz[zindex]*this->u->cval(cindex,1,0));
+ tmp[1][1] = (this->kk->kz[zindex]*this->u->cval(cindex,0,0) - this->kk->kx[xindex]*this->u->cval(cindex,2,0));
+ tmp[2][1] = (this->kk->kx[xindex]*this->u->cval(cindex,1,0) - this->kk->ky[yindex]*this->u->cval(cindex,0,0));
+ }
+ for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+ this->u->cval(cindex, cc, i) = tmp[cc][i];
+ }
+ );
+ this->kk->template force_divfree<rnumber>(this->u->get_cdata());
+ this->u->symmetrize();
+}
+
+template <class rnumber, field_backend be>
+void ou_vorticity_equation<rnumber, be>::add_ou_forcing_velocity()
+{
+ if (this->ou_forcing_type == "replace"){
+ this->ou->add_to_field_replace(this->u,"velocity");
+ }
+}
+
+template <class rnumber, field_backend be>
+void ou_vorticity_equation<rnumber, be>::add_ou_forcing_vorticity()
+{
+ if (this->ou_forcing_type == "replace"){
+ this->ou->add_to_field_replace(this->cvorticity,"vorticity");
+ }
+}
+
+template class ou_vorticity_equation<float, FFTW>;
+template class ou_vorticity_equation<double, FFTW>;
diff --git a/cpp/full_code/ou_vorticity_equation.hpp b/cpp/full_code/ou_vorticity_equation.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e36292c2271616794815d43d084810b9daeb3fc9
--- /dev/null
+++ b/cpp/full_code/ou_vorticity_equation.hpp
@@ -0,0 +1,53 @@
+#ifndef OU_VORTICITY_EQUATION_HPP
+#define OU_VORTICITY_EQUATION_HPP
+
+#include <cstdlib>
+#include "vorticity_equation.hpp"
+#include "ornstein_uhlenbeck_process.hpp"
+
+
+// vorticity_equation<double, FFTW> *test;
+
+template <typename rnumber, field_backend be>
+class ou_vorticity_equation : public vorticity_equation<rnumber, be>
+{
+ public:
+
+ std::string ou_forcing_type;
+
+ ornstein_uhlenbeck_process<rnumber,be> *ou;
+ ou_vorticity_equation(
+ const char *NAME,
+ int nx,
+ int ny,
+ int nz,
+ double ou_kmin,
+ double ou_kmax,
+ double ou_energy_amplitude,
+ double DKX = 1.0,
+ double DKY = 1.0,
+ double DKZ = 1.0,
+ unsigned FFTW_PLAN_RIGOR = FFTW_MEASURE):
+ vorticity_equation<rnumber,be>(
+ NAME, nx, ny, nz,
+ DKX,DKY,DKZ, FFTW_PLAN_RIGOR){
+ this->ou = new ornstein_uhlenbeck_process<rnumber,be>(
+ NAME,
+ nx, ny, nz,
+ ou_kmin, ou_kmax, ou_energy_amplitude,
+ DKX, DKY, DKZ, FFTW_PLAN_RIGOR);
+ this->ou_forcing_type = "replace";
+ }
+
+ ~ou_vorticity_equation();
+
+ void omega_nonlin(int src);
+ void step(double dt);
+
+ void add_ou_forcing_velocity(void);
+ void add_ou_forcing_vorticity(void);
+
+};
+
+
+#endif
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 5a3aef39caa2824f4d08e579d35734a1438ba5ec..7e6869d84f4cb6283721f1c321d2a7aaf258401b 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -1,9 +1,11 @@
#include "hdf5_tools.hpp"
+#include "scope_timer.hpp"
#include <cfloat>
#include <climits>
int hdf5_tools::require_size_single_dataset(hid_t dset, int tsize)
{
+ TIMEZONE("hdf5_tools::require_size_single_dataset");
int ndims;
hsize_t space;
space = H5Dget_space(dset);
@@ -24,6 +26,7 @@ int hdf5_tools::require_size_single_dataset(hid_t dset, int tsize)
int hdf5_tools::grow_single_dataset(hid_t dset, int tincrement)
{
+ TIMEZONE("hdf5_tools::grow_single_dataset");
int ndims;
hsize_t space;
space = H5Dget_space(dset);
@@ -48,6 +51,7 @@ herr_t hdf5_tools::require_size_dataset_visitor(
const H5O_info_t *info,
void *op_data)
{
+ TIMEZONE("hdf5_tools::require_size_dataset_visitor");
if (info->type == H5O_TYPE_DATASET)
{
hsize_t dset = H5Dopen(o_id, name, H5P_DEFAULT);
@@ -63,6 +67,7 @@ herr_t hdf5_tools::grow_dataset_visitor(
const H5O_info_t *info,
void *op_data)
{
+ TIMEZONE("hdf5_tools::grow_dataset_visitor");
if (info->type == H5O_TYPE_DATASET)
{
hsize_t dset = H5Dopen(o_id, name, H5P_DEFAULT);
@@ -78,6 +83,7 @@ int hdf5_tools::grow_file_datasets(
const std::string group_name,
int tincrement)
{
+ TIMEZONE("hdf5_tools::grow_file_datasets");
int file_problems = 0;
hid_t group;
@@ -98,6 +104,7 @@ int hdf5_tools::require_size_file_datasets(
const std::string group_name,
int tsize)
{
+ TIMEZONE("hdf5_tools::require_size_file_datasets");
int file_problems = 0;
hid_t group;
@@ -117,6 +124,7 @@ std::vector<number> hdf5_tools::read_vector(
const hid_t group,
const std::string dset_name)
{
+ TIMEZONE("hdf5_tools::read_vector");
std::vector<number> result;
hsize_t vector_length;
// first, read size of array
@@ -143,6 +151,7 @@ number hdf5_tools::read_value(
const hid_t group,
const std::string dset_name)
{
+ TIMEZONE("hdf5_tools::read_value");
number result;
hid_t dset;
hid_t mem_dtype;
@@ -177,6 +186,7 @@ std::vector<dtype> hdf5_tools::read_vector_with_single_rank(
const hid_t file_id,
const std::string dset_name)
{
+ TIMEZONE("hdf5_tools::read_vector_with_single_rank");
std::vector<dtype> data;
int vector_size;
if (myrank == rank_to_use)
@@ -208,6 +218,7 @@ std::string hdf5_tools::read_string(
const hid_t group,
const std::string dset_name)
{
+ TIMEZONE("hdf5_tools::read_string");
if (H5Lexists(group, dset_name.c_str(), H5P_DEFAULT))
{
hid_t dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
diff --git a/cpp/vorticity_equation.hpp b/cpp/vorticity_equation.hpp
index cbff223e2d8119b37ef30b4e6b739aa64eff7039..09119cc48db3c37c4faf3c54162dbe5ebaea99e9 100644
--- a/cpp/vorticity_equation.hpp
+++ b/cpp/vorticity_equation.hpp
@@ -87,8 +87,8 @@ class vorticity_equation
~vorticity_equation(void);
/* solver essential methods */
- void omega_nonlin(int src);
- void step(double dt);
+ virtual void omega_nonlin(int src);
+ virtual void step(double dt);
void impose_zero_modes(void);
/** \brief Method that computes force and adds it to the right hand side of the NS equations.
@@ -166,4 +166,3 @@ class vorticity_equation
};
#endif//VORTICITY_EQUATION
-