ornstein_uhlenbeck_process.cpp 9.24 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
191
192
193
194
195
196
197

template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
{
  //add some logic here TODO
  for(int i=0; i<1000; i++)
  {
    this->step(2e-3);
  }
}

sniklas142's avatar
sniklas142 committed
198
199
200
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
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;
                  }
                }
              }

      }

  );

}

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

238
239
240
241
242
243
244
  this->kk->CLOOP_K2(
          [&](ptrdiff_t cindex,
            ptrdiff_t xindex,
            ptrdiff_t yindex,
            ptrdiff_t zindex,
            double k2){

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

sniklas142's avatar
sniklas142 committed
247
248
249
250
                rnumber tmp;

                for(int cc=0; cc < 3; cc++){
                  for(int imag=0; imag < 2; imag++){
251
                    tmp = field_to_replace->cval(cindex,cc,imag);
sniklas142's avatar
sniklas142 committed
252
253
254
                    src->cval(cindex,cc,imag) = tmp;
                  }
                }
255
256
257
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
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();
}
287
288
289

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