ornstein_uhlenbeck_process.hpp 2.11 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#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;
23
24
        double kolmogorov_constant = 2;
        double epsilon;
25
26

        field<rnumber,be,THREE> *ou_field;
27
        field<rnumber,be,THREE> *ou_field_vort;
28
29
        field<rnumber,be,THREExTHREE> *B;
        kspace<be,SMOOTH> *kk;
Niklas Schnierstein's avatar
Niklas Schnierstein committed
30
31
32

        std::vector<double> nom; //number of modes
        std::vector<double> aom; //absolute k value of modes
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52

        std::random_device rd{};
        std::mt19937 gen{rd()};
        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);
53
        inline double energy(double kabs)
54
        {
55
            return this->ou_energy_amplitude * pow(kabs,-5./3);
56
        }
57
58
59

        inline double gamma(double kabs)
        {
Niklas Schnierstein's avatar
Niklas Schnierstein committed
60
61
            //return kabs*sqrt(this->energy(kabs));
            return 0.3333;
62
63
64
        }
        

65
        void initialize_B(void);
66

67
68
        void let_converge(void);

69
70
71
72
73
        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(
74
          field<rnumber,be,THREE> *src, std::string uv);
75

sniklas142's avatar
sniklas142 committed
76
77
78
        void strip_from_field(
          field<rnumber, be, THREE> *src);

79
80
        void calc_ou_vorticity(void);

Niklas Schnierstein's avatar
Niklas Schnierstein committed
81
82
        void find_number_of_modes(void);

83
84
85
86
87

};


#endif