diff --git a/CMakeLists.txt b/CMakeLists.txt
index 190fd619a9288ec805dfec5138d6d25e85dc1acc..78b106f95b3576d27ad8d461cb12b3ac7da02213 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -214,9 +214,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.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/ornstein_uhlenbeck_test.cpp
- ${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.cpp
- ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_ou_forcing.cpp)
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.cpp)
set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/code_base.hpp
@@ -293,9 +291,7 @@ set(hpp_for_lib
${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/ornstein_uhlenbeck_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.hpp
- ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_ou_forcing.hpp
)
#file(GLOB_RECURSE hpp_for_lib ${PROJECT_SOURCE_DIR}/*.hpp)
LIST(APPEND source_files ${hpp_for_lib} ${cpp_for_lib})
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 45e4ed0030e6094b0494bd546a3b517b363995f4..ca8e06b1c0d0033a545bba35b46efa49e8bb52bf 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -129,13 +129,6 @@ class TEST(_code):
pars['tracers0_integration_steps'] = int(4)
pars['tracers0_neighbours'] = int(1)
pars['tracers0_smoothness'] = int(1)
-
- # if dns_type == 'ornstein_uhlenbeck_test':
- # pars['ou_kmin'] = 0.0
- # pars['ou_kmax'] = 10.0
- # pars['ou_energy_amplitude'] = 1.0
- #TODO
-
return pars
def get_kspace(self):
kspace = {}
@@ -377,15 +370,6 @@ class TEST(_code):
ofile.require_group('tracers0')
for kk in ['position', 'velocity', 'vorticity', 'velocity_gradient']:
ofile['tracers0'].require_group(kk)
- if (self.dns_type == 'ornstein_uhlenbeck_test'):
- with h5py.File(os.path.join(self.work_dir, self.simname + '_ou_output.h5'), 'a') as ofile:
- ofile.require_group('spectra')
- ofile.require_group('moments')
- ofile.require_group('histograms')
- ofile['spectra'].create_dataset('ou_spectra', (1,300,3,3),dtype=np.float)
- ofile['moments'].create_dataset('ou', (1,10,3),dtype=np.float)
- ofile['histograms'].create_dataset('ou', (1,128,3),dtype=np.float)
-
self.run(
diff --git a/TurTLE/ou_quick_plot.py b/TurTLE/ou_quick_plot.py
deleted file mode 100644
index fcf5766e9d7774029caff99768d7026e3160f21a..0000000000000000000000000000000000000000
--- a/TurTLE/ou_quick_plot.py
+++ /dev/null
@@ -1,64 +0,0 @@
-import numpy as np
-import h5py
-import matplotlib.pyplot as plt
-
-
-def plot_spectrum():
- f = h5py.File('test_ou_output.h5','r')
- dset = np.array(f['spectra/ou_spectra'])
- s = dset.reshape((300,3,3))
- print('shape of s is {}'.format(s.shape))
- trace = s[:,0,0] + s[:,1,1] + s[:,2,2]
- print('shape of trace is {}'.format(trace.shape))
-
- x = np.linspace(1,300,10000)
- fivethird = lambda t: t**(-5./3.)
-
- plt.figure(0)
- plt.plot(trace)
- plt.plot(x,fivethird(x))
-
-def plot_histogram():
- f = h5py.File('test_ou_output.h5','r')
- dset = np.array(f['histograms/ou'])
- x_vel = dset[0,:,0]
- y_vel = dset[0,:,1]
- z_vel = dset[0,:,2]
- plt.figure(1)
- plt.plot(x_vel)
- plt.plot(y_vel)
- plt.plot(z_vel)
-
-def plot_field(file_name):
- f = h5py.File(file_name,'r')
- dset = np.array(f['ou_field/real/0'])
- x_vel = dset[:,:,:,0]
- y_vel = dset[:,:,:,1]
- z_vel = dset[:,:,:,2]
- plt.figure(2)
- plt.imshow(x_vel[:,:,15])
- plt.colorbar()
-
-class plot_field:
- def __init__(self):
- self.count = 2;
-
- def plot(self,file_name):
- f = h5py.File(file_name,'r')
- dset = np.array(f['ou_field/real/0'])
- x_vel = dset[:,:,:,0]
- y_vel = dset[:,:,:,1]
- z_vel = dset[:,:,:,2]
- plt.figure(self.count)
- plt.imshow(x_vel[:,:,15])
- plt.colorbar()
- self.count+=1
-
-
-
-plot_spectrum()
-plot_histogram()
-p = plot_field()
-#p.plot('ou_field.h5')
-#p.plot('add_ou_test_field.h5')
-plt.show()
diff --git a/cpp/full_code/NSVE_ou_forcing.cpp b/cpp/full_code/NSVE_ou_forcing.cpp
deleted file mode 100644
index d80fb7332ce76a6a16a3f0479238dc56d6e0d8e8..0000000000000000000000000000000000000000
--- a/cpp/full_code/NSVE_ou_forcing.cpp
+++ /dev/null
@@ -1,36 +0,0 @@
-#include <string>
-#include "NSVE_ou_forcing.hpp"
-#include "fftw_tools.hpp"
-
-
-
-template <typename rnumber>
-int NSVE_ou_forcing<rnumber>::read_parameters(void)
-{
- this->NSVE<rnumber>::read_parameters();
- return EXIT_SUCCESS;
-}
-
-template <typename rnumber>
-int NSVE_ou_forcing<rnumber>::initialize(void)
-{
- this->NSVE<rnumber>::initialize();
- this->fs = new ou_vorticity_equation<rnumber,FFTW>(
- this->simname.c_str(),
- this->nx, this->ny, this->nz,
- this->dkx, this->dky, this->dkz,
- this->ou_kmin, this->ou_kmax, this->ou_energy_amplitude,
- fftw_planner_string_to_flag[this->fftw_plan_rigor]);
-
- return EXIT_SUCCESS;
-}
-
-template <typename rnumber>
-int NSVE_ou_forcing<rnumber>::step(void){
- this->NSVE<rnumber>::step();
- //TODO STEP OF OU????
- return EXIT_SUCCESS;
-}
-
-template class NSVE_ou_forcing<float>;
-template class NSVE_ou_forcing<double>;
diff --git a/cpp/full_code/NSVE_ou_forcing.hpp b/cpp/full_code/NSVE_ou_forcing.hpp
deleted file mode 100644
index a93b636fc7bac20b6a3b813989269f40ed72544d..0000000000000000000000000000000000000000
--- a/cpp/full_code/NSVE_ou_forcing.hpp
+++ /dev/null
@@ -1,37 +0,0 @@
-#ifndef NSVE_OU_FORCING_HPP
-#define NSVE_OU_FORCING_HPP
-
-#include <cstdlib>
-#include "base.hpp"
-#include "ou_vorticity_equation.hpp"
-#include "NSVE.hpp"
-
-template <typename rnumber>
-class NSVE_ou_forcing: public NSVE<rnumber>
-{
- public:
-
- double ou_kmax;
- double ou_kmin;
- double ou_energy_amplitude;
-
- ou_vorticity_equation<rnumber,FFTW> *fs;
-
-
- NSVE_ou_forcing(
- const MPI_Comm COMMUNICATOR,
- const std::string &simulation_name):
- NSVE<rnumber>(
- COMMUNICATOR,
- simulation_name),
- ou_kmax(5.0), ou_kmin(2.0), ou_energy_amplitude(1.0){}
-
- ~NSVE_ou_forcing();
-
- int read_parameters(void);
- int step(void);
- int initialize(void);
-
-};
-
-#endif
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index a73e48e465ca221855e7459f72b34c661df29ba1..b6af42d324a27cddedd56aefbf7c57632ce186c6 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -28,7 +28,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
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(); //have it complex all zero TODO:Ask Chichi
+ this->ou_field->dft();
this->ou_field_vort = new field<rnumber,be,THREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
@@ -46,7 +46,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
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); //I can avoid one if statement later in the step method
+ assert(this->kk->kM2 >= this->ou_kmax_squ);
gen.resize(omp_get_max_threads());
@@ -110,9 +110,6 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
}
-// template <class rnumber,field_backend be>
-// void fill_number_of_modes<rnumber,be>::fill_number_of_modes(
-// double dk, double kmax
template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
@@ -137,9 +134,6 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
}
else
{
- //sigma =
- // sqrt(4*this->gamma(kabs)*this->energy(kabs)
- // /this->kk->nshell[(int)(kabs/this->kk->dk)]);
sigma =
sqrt(4*this->gamma(kabs)*this->energy(kabs)
/(4*M_PI*ksqu));
@@ -181,42 +175,6 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
});
}
-template <class rnumber, field_backend be>
-void ornstein_uhlenbeck_process<rnumber,be>::find_number_of_modes()
-{
-
- 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 kabs = sqrt(k2);
-
- ptrdiff_t pos = std::distance(std::find(aom.begin(), aom.end(), kabs), aom.begin());
-
- if (pos>=aom.size())
- {
- aom.push_back(kabs);
- nom.push_back(1);
- }
-
- else
- {
- nom[pos] += 1;
- }
-
-
- }
-
- }
-
- );
-}
template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
@@ -292,55 +250,6 @@ void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_replace(
);
}
-template <class rnumber, field_backend be>
-void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_sharp(
- 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)
- {
- src->cval(cindex,0,0) += this->ou_field->cval(cindex,0,0);
- src->cval(cindex,1,0) += this->ou_field->cval(cindex,1,0);
- src->cval(cindex,2,0) += this->ou_field->cval(cindex,2,0);
-
- src->cval(cindex,0,1) += this->ou_field->cval(cindex,0,1);
- src->cval(cindex,1,1) += this->ou_field->cval(cindex,1,1);
- src->cval(cindex,2,1) += this->ou_field->cval(cindex,2,1);
- }
-
- }
-
- );
-}
-
-template <class rnumber, field_backend be>
-void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_gaussian(
- field<rnumber, be, THREE> *src, double param)
-{
- 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)
- {
-
- }
-
- }
-
- );
-}
template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_vorticity(void)
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index 7710b46702ead3502ab93a7a6056be20d62ed542..a1aeaa98dd948ca39b2c96177341890b345c9953 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -6,7 +6,6 @@
#include <stdlib.h>
#include <iostream>
#include "field.hpp"
-// #include <unordered_map>
#include <random>
extern int myrank, nprocs;
@@ -28,11 +27,7 @@ class ornstein_uhlenbeck_process{
field<rnumber,be,THREExTHREE> *B;
kspace<be,SMOOTH> *kk;
- std::vector<double> nom; //number of modes
- std::vector<double> aom; //absolute k value of modes
-
- //std::random_device rd{};
- std::vector<std::mt19937_64> gen; //32-bit, maybe fix for double prec TODO
+ std::vector<std::mt19937_64> gen;
std::normal_distribution<double> d{0,1};
ornstein_uhlenbeck_process(
@@ -57,8 +52,8 @@ class ornstein_uhlenbeck_process{
inline double gamma(double kabs)
{
- //return kabs*sqrt(this->energy(kabs));
- return kolmogorov_constant*sqrt(kolmogorov_constant/this->ou_energy_amplitude);
+ return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
+ this->ou_energy_amplitude);
}
@@ -66,10 +61,6 @@ class ornstein_uhlenbeck_process{
void let_converge(void);
- void add_to_field_gaussian(
- field<rnumber,be,THREE> *src, double param);
- void add_to_field_sharp(
- field<rnumber,be,THREE> *src);
void add_to_field_replace(
field<rnumber,be,THREE> *src, std::string uv);
@@ -78,8 +69,6 @@ class ornstein_uhlenbeck_process{
void calc_ou_vorticity(void);
- void find_number_of_modes(void);
-
};
diff --git a/cpp/full_code/ornstein_uhlenbeck_test.cpp b/cpp/full_code/ornstein_uhlenbeck_test.cpp
deleted file mode 100644
index aaaa2fba3722c6030f8559ba45d363b57826cea4..0000000000000000000000000000000000000000
--- a/cpp/full_code/ornstein_uhlenbeck_test.cpp
+++ /dev/null
@@ -1,123 +0,0 @@
-#include <string>
-#include <cmath>
-#include "ornstein_uhlenbeck_test.hpp"
-#include "scope_timer.hpp"
-
-
-template <typename rnumber>
-int ornstein_uhlenbeck_test<rnumber>::initialize(void)
-{
- TIMEZONE("ornstein_uhlenbeck_test::initialize");
- this->read_parameters();
- this->ou = new ornstein_uhlenbeck_process<rnumber,FFTW>(
- "test",
- this->nx,this->ny,this->nz,
- this->ou_kmin, this->ou_kmax, this->ou_energy_amplitude,
- this->dkx,this->dky,this->dkz,
- FFTW_ESTIMATE);
- if (this->myrank == 0)
- {
- hid_t stat_file = H5Fopen(
- (this->simname + std::string(".h5")).c_str(),
- H5F_ACC_RDWR,
- H5P_DEFAULT);
- this->ou->kk->store(stat_file);
- H5Fclose(stat_file);
- }
- return EXIT_SUCCESS;
-}
-
-template <typename rnumber>
-int ornstein_uhlenbeck_test<rnumber>::finalize(void)
-{
- TIMEZONE("ornstein_uhlenbeck_test::finalize");
- delete this->ou;
- return EXIT_SUCCESS;
-}
-
-template <typename rnumber>
-int ornstein_uhlenbeck_test<rnumber>::read_parameters()
-{
- TIMEZONE("ornstein_uhlenbeck_test::read_parameters");
- this->test::read_parameters();
- // hid_t parameter_file = H5Fopen(
- // (this->simname + std::string(".h5")).c_str(),
- // H5F_ACC_RDONLY,
- // H5P_DEFAULT);
- // this->ou_kmin = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_kmin");
- // this->ou_kmax = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_kmax");
- // this->ou_energy_amplitude = hdf5_tools::read_value<double>(parameter_file, "/parameters/ou_energy_amplitude");
- // H5Fclose(parameter_file);
- //TODO
- this->ou_kmin = 3;
- this->ou_kmax = 14;
- this->ou_energy_amplitude = 1;
-
-
- return EXIT_SUCCESS;
-}
-
-template <typename rnumber>
-int ornstein_uhlenbeck_test<rnumber>::do_work(void)
-{
- TIMEZONE("ornstein_uhlenbeck_test::do_work");
-
- for (int step = 0; step < 10000; step++)
- {
- this->ou->step(1e-3);
- }
-
- this->calc_stats();
- this->test_summation();
-
- return EXIT_SUCCESS;
-}
-
-template <typename rnumber>
-void ornstein_uhlenbeck_test<rnumber>::calc_stats(void)
-{
- TIMEZONE("ornstein_uhlenbeck_test::calc_stats");
- field<rnumber,FFTW,THREE> *tmp_field = new field<rnumber,FFTW,THREE>(
- this->nx,this->ny,this->nz, MPI_COMM_WORLD, FFTW_ESTIMATE);
- *tmp_field=0.0;
- tmp_field->dft();
- *tmp_field = *this->ou->ou_field;
- std::cerr << tmp_field << " " << this->ou->ou_field << std::endl;
- std::string filename = this->simname + std::string("_fields.h5");
- hid_t out_file;
-
- if (this->myrank == 0) {
- out_file = H5Fopen("test_ou_output.h5", H5F_ACC_RDWR, H5P_DEFAULT);
- }
-
- this->ou->kk->template cospectrum<rnumber,THREE>(
- tmp_field->get_cdata(),
- tmp_field->get_cdata(),
- out_file,"ou_spectra",0);
-
- tmp_field->symmetrize(); //TODO necessary, not needed, introduce error?
- tmp_field->ift();
- std::vector<double> me;
- me.resize(3);
- std::fill(me.begin(),me.end(),3);
- tmp_field->compute_rspace_stats(out_file, "ou", 0, me);
- tmp_field->io("ou_field.h5", "ou_field", 0, false);
-}
-
-template <typename rnumber>
-void ornstein_uhlenbeck_test<rnumber>::test_summation()
-{
- field<rnumber,FFTW,THREE> *test_f = new field<rnumber,FFTW,THREE>(
- this->nx,this->ny,this->nz, MPI_COMM_WORLD, FFTW_ESTIMATE);
- *test_f = 0.0;
- test_f->dft();
- this->ou->add_to_field_replace(test_f, "velocity");
- test_f->symmetrize();
- test_f->ift();
- test_f->io("add_ou_test_field.h5", "ou_field", 0, false);
-
-}
-
-
-template class ornstein_uhlenbeck_test<float>;
-template class ornstein_uhlenbeck_test<double>;
diff --git a/cpp/full_code/ornstein_uhlenbeck_test.hpp b/cpp/full_code/ornstein_uhlenbeck_test.hpp
deleted file mode 100644
index 2df847ff959607b8f98de98cd506c33eeef4447e..0000000000000000000000000000000000000000
--- a/cpp/full_code/ornstein_uhlenbeck_test.hpp
+++ /dev/null
@@ -1,45 +0,0 @@
-#ifndef ORNSTEIN_UHLENBECK_TEST_HPP
-#define ORNSTEIN_UHLENBECK_TEST_HPP
-
-
-
-#include <cstdlib>
-#include "base.hpp"
-#include "kspace.hpp"
-#include "field.hpp"
-#include "full_code/test.hpp"
-#include "ornstein_uhlenbeck_process.hpp"
-
-
-
-template <typename rnumber>
-class ornstein_uhlenbeck_test: public test {
-
- public:
-
- double ou_kmin,ou_kmax;
- double ou_energy_amplitude;
-
- ornstein_uhlenbeck_process<rnumber,FFTW> *ou;
-
- ornstein_uhlenbeck_test(
- const MPI_Comm COMMUNICATOR,
- const std::string &simulation_name):
- test(
- COMMUNICATOR,
- simulation_name){}
- ~ornstein_uhlenbeck_test(){}
-
- int initialize(void);
- int do_work(void);
- int finalize(void);
- int read_parameters(void);
-
- void calc_stats(void);
- void test_summation(void);
-
-
-
-};
-
-#endif
diff --git a/cpp/full_code/ou_vorticity_equation.cpp b/cpp/full_code/ou_vorticity_equation.cpp
index f3bfe8c45db11c42eab7faf8ad6361803a1b2f39..37d3570d6f0c19849f5b584d4eaac5c7ed60c43f 100644
--- a/cpp/full_code/ou_vorticity_equation.cpp
+++ b/cpp/full_code/ou_vorticity_equation.cpp
@@ -38,7 +38,6 @@ void ou_vorticity_equation<rnumber, be>::step(double dt)
}
}
);
- //this->impose_forcing(this->v[1], this->v[0]);
this->omega_nonlin(1);
this->kk->CLOOP_K2(
@@ -60,7 +59,6 @@ void ou_vorticity_equation<rnumber, be>::step(double dt)
}
}
);
- //this->impose_forcing(this->v[2], this->v[0]);
this->omega_nonlin(2);
// store old vorticity
@@ -83,7 +81,6 @@ void ou_vorticity_equation<rnumber, be>::step(double dt)
}
}
);
- //this->impose_forcing(this->v[0], this->v[1]);
this->kk->template force_divfree<rnumber>(this->cvorticity->get_cdata());
this->cvorticity->symmetrize();
@@ -99,7 +96,7 @@ void ou_vorticity_equation<rnumber, be>::omega_nonlin(int src)
this->compute_velocity(this->v[src]);
this->add_ou_forcing_velocity();
- // TIME STEP?!?!?!!?! TODO
+ // TIME STEP TODO
/* get fields from Fourier space to real space */
this->u->ift();
@@ -143,7 +140,6 @@ void ou_vorticity_equation<rnumber, be>::omega_nonlin(int src)
this->u->cval(cindex, cc, i) = tmp[cc][i];
}
);
- //this->add_forcing(this->u, this->v[src]); TODO not necessary for ou?
this->kk->template force_divfree<rnumber>(this->u->get_cdata());
this->u->symmetrize();
}