-
Niklas Schnierstein authoredNiklas Schnierstein authored
ornstein_uhlenbeck_process.hpp 2.21 KiB
#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 <unordered_map>
#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<double> nom; //number of modes
std::vector<double> aom; //absolute k value of modes
std::random_device rd{};
std::mt19937 gen{rd()}; //32-bit, maybe fix for double prec TODO
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 kabs*sqrt(this->energy(kabs));
return kolmogorov_constant*sqrt(kolmogorov_constant/this->ou_energy_amplitude);
}
void initialize_B(void);
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);
void strip_from_field(
field<rnumber, be, THREE> *src);
void calc_ou_vorticity(void);
void find_number_of_modes(void);
};
#endif