diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp index 653c2540f01ee18661346b04945c6cade09450d2..95a9cf9b385be1214db9639c95df347485a00a34 100644 --- a/cpp/full_code/ornstein_uhlenbeck_process.cpp +++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp @@ -186,18 +186,54 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B() template <class rnumber, field_backend be> -void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void) +void ornstein_uhlenbeck_process<rnumber, be>::setup_field(int iteration,field<rnumber, be, THREE> *src) { - 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++) - { - this->step(dt); + if (iteration == 0){ + 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++) + { + this->step(dt); + } + + } + else { + this->set_from_field(src); } } +template <class rnumber, field_backend be> +void ornstein_uhlenbeck_process<rnumber, be>::set_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++){ + this->ou_field_vort->cval(cindex,cc,imag) = src->cval(cindex,cc,imag); + } + } + } + + } + + ); + + this->calc_ou_velocity(); + +} + template <class rnumber, field_backend be> void ornstein_uhlenbeck_process<rnumber, be>::strip_from_field( field<rnumber, be, THREE> *src) @@ -288,5 +324,41 @@ void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_vorticity(void) this->ou_field_vort->symmetrize(); } +template <class rnumber, field_backend be> +void ornstein_uhlenbeck_process<rnumber,be>::calc_ou_velocity(void) +{ + 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) + { + this->ou_field->cval(cindex,0,0) = + -(this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,2,1) - this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,1,1)) / k2; + + this->ou_field->cval(cindex,0,1) = + (this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,2,0) - this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,1,0)) / k2; + + this->ou_field->cval(cindex,1,0) = + -(this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,0,1) - this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,2,1)) / k2; + + this->ou_field->cval(cindex,1,1) = + (this->kk->kz[zindex]*this->ou_field_vort->cval(cindex,0,0) - this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,2,0)) / k2; + + this->ou_field->cval(cindex,2,0) = + -(this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,1,1) - this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,0,1)) / k2; + + this->ou_field->cval(cindex,2,1) = + (this->kk->kx[xindex]*this->ou_field_vort->cval(cindex,1,0) - this->kk->ky[yindex]*this->ou_field_vort->cval(cindex,0,0)) / k2; + } + else + std::fill_n((rnumber*)(this->ou_field_vort->get_cdata()+3*cindex), 6, 0.0); + } + ); + this->ou_field->symmetrize(); +} + template class ornstein_uhlenbeck_process<float,FFTW>; template class ornstein_uhlenbeck_process<double,FFTW>; diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp index 5a1288c929afac4cad4a3a6f36ff5b2a0436b1ca..33417c976ed238a77b2406e6e03bb8c3aa4be5a7 100644 --- a/cpp/full_code/ornstein_uhlenbeck_process.hpp +++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp @@ -60,7 +60,7 @@ class ornstein_uhlenbeck_process{ void initialize_B(void); - void let_converge(void); + void setup_field(int iteration,field<rnumber, be, THREE> *src); void add_to_field_replace( field<rnumber,be,THREE> *src, std::string uv); @@ -68,7 +68,11 @@ class ornstein_uhlenbeck_process{ void strip_from_field( field<rnumber, be, THREE> *src); + void set_from_field( + field<rnumber, be, THREE> *src); + void calc_ou_vorticity(void); + void calc_ou_velocity(void); };