ornstein_uhlenbeck_process.cpp 9.32 KB
Newer Older
1
2
3
4
5
#include "ornstein_uhlenbeck_process.hpp"
#include <cmath>
#include <cstring>
#include <cassert>
#include "scope_timer.hpp"
Niklas Schnierstein's avatar
Niklas Schnierstein committed
6
#include <algorithm>
7
#include <chrono>
8
9
10
11
12
13
14
15
16
17
18


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,
19
    double ou_gamma_factor,
20
21
22
23
24
25
26
27
28
29
30
31
    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;
32
    this->ou_field->dft();
33
34
35
36
37
38

    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();

39
40
41
42
43
44
45
46
47
    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;
48
    this->ou_gamma_factor = ou_gamma_factor;
49
50
    this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);

51
    assert(this->kk->kM2 >= this->ou_kmax_squ);
Niklas Schnierstein's avatar
Niklas Schnierstein committed
52

53
    gen.resize(omp_get_max_threads());
Niklas Schnierstein's avatar
Niklas Schnierstein committed
54

55
56
57
58
59
60
61
    long now;

    if (this->ou_field->clayout->myrank == 0){
	now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
    }
    MPI_Bcast(&now,1,MPI_LONG,0,this->ou_field->comm);

62
63
    for(int thread=0;thread<omp_get_max_threads();thread++)
    {
64
            long current_seed =
65
                    this->ou_field->clayout->myrank*omp_get_max_threads() +
66
                    thread+now;
67
68
69
            gen[thread].seed(current_seed);
    }

70

71
72
73
74
75
76
77
78
79
80
    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;
81
    delete this->ou_field_vort;
82
83
84
85
86
87
88
    delete this->B;
}

template <class rnumber,field_backend be>
void ornstein_uhlenbeck_process<rnumber,be>::step(
    double dt)
{
89
    // put in "CFL"-criterium TODO!!!
90
91
92
93
94
95
96
97
98
    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)
        {
99
            double g = this->gamma(sqrt(k2));
100
            int tr = omp_get_thread_num();
101
102
            double random[6] =
            {
103
104
                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])
105
106
107
108
109
110
111
112
113
114
115
            };
            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] );

        }
    });
116

Niklas Schnierstein's avatar
Niklas Schnierstein committed
117
    this->ou_field->symmetrize();
118
119
    this->calc_ou_vorticity();

120
121
122
123
124
125
126
}


template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
{
    TIMEZONE("ornstein_uhlenbeck_process::initialize_B");
127
    *this->B = 0.0;
128
129
130
131
132
133
134
135
    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);
136
137
138
            if (ksqu <= this->ou_kmax_squ && ksqu >= this->ou_kmin_squ)
            {
                double kabs = sqrt(ksqu);
Niklas Schnierstein's avatar
Niklas Schnierstein committed
139
                double sigma;
sniklas142's avatar
sniklas142 committed
140
                if (ksqu == 0 || ksqu == this->ou_kmax_squ)
Niklas Schnierstein's avatar
Niklas Schnierstein committed
141
142
143
144
145
146
                {
                    ksqu = 1; kabs = 1; sigma = 0;
                }
                else
                {
                    sigma =
147
                        sqrt(4*this->gamma(kabs)*this->energy(kabs)
Niklas Schnierstein's avatar
Niklas Schnierstein committed
148
                        /(4*M_PI*ksqu));
Niklas Schnierstein's avatar
Niklas Schnierstein committed
149
                }
150
151
152
153
154
155
156

                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);
157

158
159
160
                        if (i+j == 4)
                            this->B->cval(cindex,i,j,0) =
                                sigma/2. * (1-pow(this->kk->kz[zindex],2)/ksqu);
161

162
163
164
                        if (i+j == 1)
                            this->B->cval(cindex,i,j,0) =
                                sigma/2. * (0-this->kk->kx[xindex]*this->kk->ky[yindex]/ksqu);
165

166
167
168
                        if (i+j == 3)
                            this->B->cval(cindex,i,j,0) =
                                sigma/2. * (0-this->kk->kz[zindex]*this->kk->ky[yindex]/ksqu);
169

170
                        if (i+j == 2) {
171

172
173
174
                            if(i==j)
                                this->B->cval(cindex,i,j,0) =
                                    sigma/2. * (1-pow(this->kk->ky[yindex],2)/ksqu);
175

176
177
178
179
                            if(i!=j)
                                this->B->cval(cindex,i,j,0) =
                                    sigma/2. * (0-this->kk->kx[xindex]*this->kk->kz[zindex]/ksqu);
                        }
180
181
182
                    }
                }
            }
183

184
185
186
        });
}

187
188
189
190

template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
{
191
192
193
194
195
  double ou_kmin = sqrt(this->ou_kmin_squ);
  double tau = 1.0/this->gamma(ou_kmin);
  double dt = tau/1000.;

  for(int i=0; i<2000; i++)
196
  {
197
    this->step(dt);
198
199
200
  }
}

sniklas142's avatar
sniklas142 committed
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
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;
                  }
                }
              }

      }

  );

}

229
230
template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_replace(
231
  field<rnumber, be, THREE> *src, std::string uv)
232
233
{
  assert(src->real_space_representation==false);
234
235
236
237
238
239
240
  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;

241
242
243
244
245
246
247
  this->kk->CLOOP_K2(
          [&](ptrdiff_t cindex,
            ptrdiff_t xindex,
            ptrdiff_t yindex,
            ptrdiff_t zindex,
            double k2){

sniklas142's avatar
sniklas142 committed
248
              if (k2 <= this->ou_kmax_squ && k2 >= this->ou_kmin_squ){
249

sniklas142's avatar
sniklas142 committed
250
251
252
253
                rnumber tmp;

                for(int cc=0; cc < 3; cc++){
                  for(int imag=0; imag < 2; imag++){
254
                    tmp = field_to_replace->cval(cindex,cc,imag);
sniklas142's avatar
sniklas142 committed
255
256
257
                    src->cval(cindex,cc,imag) = tmp;
                  }
                }
258
259
260
261
262
263
264
265
              }

      }

  );
}


266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
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();
}
290
291
292

template class ornstein_uhlenbeck_process<float,FFTW>;
template class ornstein_uhlenbeck_process<double,FFTW>;