Commit 2f650c95 authored by sniklas142's avatar sniklas142
Browse files

clean up - ready to merge

parent ba71d3ac
......@@ -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})
......
......@@ -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(
......
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()
#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>;
#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
......@@ -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)
......
......@@ -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);
};
......
#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>;
#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
......@@ -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();
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment