first implementation of OU forcing (not tested but compiles)

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/)")
-

cpp/full_code/ornstein_uhlenbeck_process.cpp

+#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>;

cpp/full_code/ornstein_uhlenbeck_process.hpp

+#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