diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c0c1ef36964f482d9bb061bbf89afa3f12d2f99..97b263614e7383d6fc637d2ddda45092af14f27b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -211,7 +211,8 @@ set(cpp_for_lib ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.cpp ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.cpp - ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp) + ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp + ${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.cpp) set(hpp_for_lib ${PROJECT_SOURCE_DIR}/cpp/full_code/code_base.hpp ${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.hpp @@ -285,6 +286,7 @@ set(hpp_for_lib ${PROJECT_SOURCE_DIR}/cpp/omputils.hpp ${PROJECT_SOURCE_DIR}/cpp/shared_array.hpp ${PROJECT_SOURCE_DIR}/cpp/spline.hpp + ${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.hpp ) #file(GLOB_RECURSE hpp_for_lib ${PROJECT_SOURCE_DIR}/*.hpp) LIST(APPEND source_files ${hpp_for_lib} ${cpp_for_lib}) @@ -316,4 +318,3 @@ else() install(CODE "execute_process(COMMAND ${CMAKE_COMMAND} -E copy ${PROJECT_SOURCE_DIR}/pc_host_info.py ${PROJECT_BINARY_DIR}/python/TurTLE/host_info.py)") endif() install(CODE "execute_process(COMMAND python3 ${PROJECT_SOURCE_DIR}/setup.py install --force --prefix=${CMAKE_INSTALL_PREFIX} WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/python/)") - diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp new file mode 100644 index 0000000000000000000000000000000000000000..37ea946e4a61151a202ced69d8a33945b3253516 --- /dev/null +++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp @@ -0,0 +1,147 @@ +#include "ornstein_uhlenbeck_process.hpp" +#include <cmath> +#include <cstring> +#include <cassert> +#include "scope_timer.hpp" +#define NDEBUG + + +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, + 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; + this->ou_field->dft(); //have it complex all zero TODO:Ask Chichi + 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; + assert(this->kk->kM2 >= this->ou_kmax_squ); //I can avoid one if statement later in the step method + 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; + delete this->B; +} + +template <class rnumber,field_backend be> +void ornstein_uhlenbeck_process<rnumber,be>::step( + double dt) +{ + 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) + { + double g = this->gamma(k2); + double random[6] = + { + this->d(this->gen),this->d(this->gen),this->d(this->gen), + this->d(this->gen),this->d(this->gen),this->d(this->gen) + }; + 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] ); + + } + }); + // this->ou_field->symmetrize(); ??? +} + +// template <class rnumber,field_backend be> +// void fill_number_of_modes<rnumber,be>::fill_number_of_modes( +// double dk, double kmax + +template <class rnumber, field_backend be> +void ornstein_uhlenbeck_process<rnumber,be>::initialize_B() +{ + TIMEZONE("ornstein_uhlenbeck_process::initialize_B"); + + 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); + double kabs = sqrt(ksqu); + double sigma = + this->ou_energy_amplitude * + sqrt(4*this->gamma(kabs)*pow(kabs,-5./3.) + /this->kk->nshell[(int)(kabs/this->kk->dk)]); + + 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); + + if (i+j == 4) + this->B->cval(cindex,i,j,0) = + sigma/2. * (1-pow(this->kk->kz[zindex],2)/ksqu); + + if (i+j == 1) + this->B->cval(cindex,i,j,0) = + sigma/2. * (0-this->kk->kx[xindex]*this->kk->ky[yindex]/ksqu); + + if (i+j == 3) + this->B->cval(cindex,i,j,0) = + sigma/2. * (0-this->kk->kz[zindex]*this->kk->ky[yindex]/ksqu); + + if (i+j == 2) { + + if(i==j) + this->B->cval(cindex,i,j,0) = + sigma/2. * (1-pow(this->kk->ky[yindex],2)/ksqu); + + if(i!=j) + this->B->cval(cindex,i,j,0) = + sigma/2. * (0-this->kk->kx[xindex]*this->kk->kz[zindex]/ksqu); + } + } + } + }); +} + + +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 new file mode 100644 index 0000000000000000000000000000000000000000..07c7af594bf7e1445917c57b7e144ad555ffffdd --- /dev/null +++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp @@ -0,0 +1,60 @@ +#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; + + field<rnumber,be,THREE> *ou_field; + field<rnumber,be,THREExTHREE> *B; + // field<rnumber,be,ONE> *number_of_modes; + kspace<be,SMOOTH> *kk; + // std::unordered_map<double,double> map_n_kabs; + + 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); + inline double gamma(double ksqu) + { + return this->ou_energy_amplitude * pow(ksqu,1./3.); + } + // double sigma(double k_abs); + void initialize_B(void); + // void fill_number_of_modes(double dk, double kmax); + +}; + + +#endif