Merge branch 'feature/ou_with_vort' into develop

CMakeLists.txt

@@ -225,7 +225,10 @@ 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
+    ${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.cpp)
+
 set(hpp_for_lib
     ${PROJECT_SOURCE_DIR}/cpp/full_code/code_base.hpp
     ${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.hpp
@@ -305,6 +308,8 @@ 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
+    ${PROJECT_SOURCE_DIR}/cpp/full_code/ou_vorticity_equation.hpp
     )
 #file(GLOB_RECURSE hpp_for_lib ${PROJECT_SOURCE_DIR}/*.hpp)
 LIST(APPEND source_files ${hpp_for_lib} ${cpp_for_lib})
@@ -336,4 +341,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/)")
-

TurTLE/DNS.py

@@ -2,20 +2,20 @@
 #                                                                              #
 #  Copyright 2015-2019 Max Planck Institute for Dynamics and Self-Organization #
 #                                                                              #
-#  This file is part of TurTLE.                                                  #
+#  This file is part of TurTLE.                                                #
 #                                                                              #
-#  TurTLE is free software: you can redistribute it and/or modify                #
+#  TurTLE is free software: you can redistribute it and/or modify              #
 #  it under the terms of the GNU General Public License as published           #
 #  by the Free Software Foundation, either version 3 of the License,           #
 #  or (at your option) any later version.                                      #
 #                                                                              #
-#  TurTLE is distributed in the hope that it will be useful,                     #
+#  TurTLE is distributed in the hope that it will be useful,                   #
 #  but WITHOUT ANY WARRANTY; without even the implied warranty of              #
 #  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               #
 #  GNU General Public License for more details.                                #
 #                                                                              #
 #  You should have received a copy of the GNU General Public License           #
-#  along with TurTLE.  If not, see <http://www.gnu.org/licenses/>                #
+#  along with TurTLE.  If not, see <http://www.gnu.org/licenses/>              #
 #                                                                              #
 # Contact: Cristian.Lalescu@ds.mpg.de                                          #
 #                                                                              #
@@ -657,7 +657,9 @@ class DNS(_code):
         parser_NSVEp_extra = subparsers.add_parser(
                 'NSVEp_extra_sampling',
                 help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers, that sample velocity gradient, as well as pressure and its derivatives.')
-
+        parser_NSVE_ou = subparsers.add_parser(
+                'NSVE_ou_forcing',
+                help = 'plain Navier-Stokes vorticity formulation, with ornstein-uhlenbeck forcing')
         for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra', 'static_field']:
             eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
             eval('self.job_parser_arguments({0})'.format('parser_' + parser))
@@ -1060,4 +1062,3 @@ class DNS(_code):
                 no_submit = opt.no_submit,
                 no_debug = opt.no_debug)
         return None
-

TurTLE/TEST.py

@@ -275,12 +275,18 @@ class TEST(_code):
         parser_test_interpolation = subparsers.add_parser(
                 'test_interpolation',
                 help = 'test velocity gradient interpolation')
+
+        parser_ornstein_uhlenbeck_test = subparsers.add_parser(
+                'ornstein_uhlenbeck_test',
+                help = 'test ornstein uhlenbeck process')
+
         for parser in ['parser_filter_test',
                        'parser_write_filtered_test',
                        'parser_field_test',
                        'parser_symmetrize_test',
                        'parser_field_output_test',
-                       'parser_test_interpolation']:
+                       'parser_test_interpolation',
+                       'parser_ornstein_uhlenbeck_test']:
             eval('self.simulation_parser_arguments(' + parser + ')')
             eval('self.job_parser_arguments(' + parser + ')')
             eval('self.parameters_to_parser_arguments(' + parser + ')')
@@ -368,6 +374,8 @@ class TEST(_code):
                     ofile.require_group('tracers0')
                     for kk in ['position', 'velocity', 'vorticity', 'velocity_gradient']:
                         ofile['tracers0'].require_group(kk)
+
+
         self.run(
                 nb_processes = opt.nb_processes,
                 nb_threads_per_process = opt.nb_threads_per_process,
@@ -376,4 +384,3 @@ class TEST(_code):
                 minutes = opt.minutes % 60,
                 no_submit = opt.no_submit)
         return None
-

cpp/full_code/NSVE.cpp

@@ -62,11 +62,11 @@ int NSVE<rnumber>::initialize(void)
     }
     this->fs = new vorticity_equation<rnumber, FFTW>(
             simname.c_str(),
-            nx, ny, nz,
-            dkx, dky, dkz,
+            this->nx, this->ny, this->nz,
+            this->dkx, this->dky, this->dkz,
             fftw_planner_string_to_flag[this->fftw_plan_rigor]);
     this->tmp_vec_field = new field<rnumber, FFTW, THREE>(
-            nx, ny, nz,
+            this->nx, this->ny, this->nz,
             this->comm,
             fftw_planner_string_to_flag[this->fftw_plan_rigor]);
 

cpp/full_code/ornstein_uhlenbeck_process.cpp

+#include "ornstein_uhlenbeck_process.hpp"
+#include <cmath>
+#include <cstring>
+#include <cassert>
+#include "scope_timer.hpp"
+#include <algorithm>
+#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(); 
+
+    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();
+
+    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;
+    this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);
+
+    assert(this->kk->kM2 >= this->ou_kmax_squ); 
+
+    gen.resize(omp_get_max_threads());
+
+    for(int thread=0;thread<omp_get_max_threads();thread++)
+    {
+            int current_seed = 
+                    this->ou_field->clayout->myrank*omp_get_max_threads() +
+                    thread;
+            gen[thread].seed(current_seed);
+    }
+
+    
+    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->ou_field_vort;
+    delete this->B;
+}
+
+template <class rnumber,field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::step(
+    double dt)
+{
+    // put in "CFL"-criterium TODO!!!
+    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(sqrt(k2));
+            int tr = omp_get_thread_num();
+            double random[6] =
+            {
+                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])
+            };
+            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();
+    this->calc_ou_vorticity();
+
+}
+
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
+{
+    TIMEZONE("ornstein_uhlenbeck_process::initialize_B");
+    *this->B = 0.0;
+    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);
+            if (ksqu <= this->ou_kmax_squ && ksqu >= this->ou_kmin_squ)
+            {
+                double kabs = sqrt(ksqu);
+                double sigma;
+                if (ksqu == 0 || ksqu == this->ou_kmax_squ)
+                {
+                    ksqu = 1; kabs = 1; sigma = 0;
+                }
+                else
+                {
+                    sigma =
+                        sqrt(4*this->gamma(kabs)*this->energy(kabs)
+                        /(4*M_PI*ksqu));
+                }
+
+                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 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);
+  }
+}
+
+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;
+                  }
+                }
+              }
+
+      }
+
+  );
+
+}
+
+template <class rnumber, field_backend be>
+void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_replace(
+  field<rnumber, be, THREE> *src, std::string uv)
+{
+  assert(src->real_space_representation==false);
+  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;
+
+  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){
+
+                rnumber tmp;
+
+                for(int cc=0; cc < 3; cc++){
+                  for(int imag=0; imag < 2; imag++){
+                    tmp = field_to_replace->cval(cindex,cc,imag);
+                    src->cval(cindex,cc,imag) = tmp;
+                  }
+                }
+              }
+
+      }
+
+  );
+}
+
+
+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();
+}
+
+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 <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;
+        double kolmogorov_constant = 2;
+        double epsilon;
+
+        field<rnumber,be,THREE> *ou_field;
+        field<rnumber,be,THREE> *ou_field_vort;
+        field<rnumber,be,THREExTHREE> *B;
+        kspace<be,SMOOTH> *kk;
+
+        std::vector<std::mt19937_64> gen; 
+        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 energy(double kabs)
+        {
+            return this->ou_energy_amplitude * pow(kabs,-5./3);
+        }
+
+        inline double gamma(double kabs)
+        {
+            return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
+                            this->ou_energy_amplitude);
+        }
+        
+
+        void initialize_B(void);
+
+        void let_converge(void);
+
+        void add_to_field_replace(
+          field<rnumber,be,THREE> *src, std::string uv);
+
+        void strip_from_field(
+          field<rnumber, be, THREE> *src);
+
+        void calc_ou_vorticity(void);
+
+
+};
+
+
+#endif

cpp/full_code/ou_vorticity_equation.cpp

+#include "ou_vorticity_equation.hpp"
+#include <cmath>
+#include <cstring>
+#include <cassert>
+#include "scope_timer.hpp"
+#define NDEBUG
+
+
+template <class rnumber, field_backend be>
+ou_vorticity_equation<rnumber,be>::~ou_vorticity_equation()
+{
+  delete this->ou;
+}
+
+template <class rnumber,
+          field_backend be>
+void ou_vorticity_equation<rnumber, be>::step(double dt)
+{
+    DEBUG_MSG("vorticity_equation::step\n");
+    TIMEZONE("vorticity_equation::step");
+    this->ou->add_to_field_replace(this->cvorticity, "vorticity");
+    *this->v[1] = 0.0;
+    this->omega_nonlin(0);
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2)
+        {
+            double factor0;
+            factor0 = exp(-this->nu * k2 * dt);
+            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+                this->v[1]->cval(cindex,cc,i) = (
+                        this->v[0]->cval(cindex,cc,i) +
+                        dt*this->u->cval(cindex,cc,i))*factor0;
+        }
+    }
+    );
+
+    this->omega_nonlin(1);
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2)
+        {
+            double factor0, factor1;
+            factor0 = exp(-this->nu * k2 * dt/2);
+            factor1 = exp( this->nu * k2 * dt/2);
+            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+                this->v[2]->cval(cindex, cc, i) = (
+                        3*this->v[0]->cval(cindex,cc,i)*factor0 +
+                        ( this->v[1]->cval(cindex,cc,i) +
+                         dt*this->u->cval(cindex,cc,i))*factor1)*0.25;
+        }
+    }
+    );
+
+    this->omega_nonlin(2);
+    // store old vorticity
+    *this->v[1] = *this->v[0];
+    this->kk->CLOOP_K2(
+                [&](ptrdiff_t cindex,
+                    ptrdiff_t xindex,
+                    ptrdiff_t yindex,
+                    ptrdiff_t zindex,
+                    double k2){
+        if (k2 <= this->kk->kM2)
+        {
+            double factor0;
+            factor0 = exp(-this->nu * k2 * dt * 0.5);
+            for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
+                this->v[3]->cval(cindex,cc,i) = (
+                        this->v[0]->cval(cindex,cc,i)*factor0 +
+                        2*(this->v[2]->cval(cindex,cc,i) +
+                           dt*this->u->cval(cindex,cc,i)))*factor0/3;
+        }
+    }
+    );
+
+    this->kk->template force_divfree<rnumber>(this->cvorticity->get_cdata());
+    this->cvorticity->symmetrize();
+    this->iteration++;
+}
+
+
+template <class rnumber, field_backend be>
+void ou_vorticity_equation<rnumber, be>::omega_nonlin(int src)
+{
+  DEBUG_MSG("ou_vorticity_equation::omega_nonlin(%d)\n", src);
+  assert(src >= 0 && src < 3);
+  this->compute_velocity(this->v[src]);
+
+  this->add_ou_forcing_velocity();
+  // TIME STEP TODO
+
+  /* get fields from Fourier space to real space */
+  this->u->ift();
+  this->rvorticity->real_space_representation = false;
+  *this->rvorticity = this->v[src]->get_cdata();
+  this->rvorticity->ift();
+  /* compute cross product $u \times \omega$, and normalize */
+  this->u->RLOOP(
+              [&](ptrdiff_t rindex,
+                  ptrdiff_t xindex,
+                  ptrdiff_t yindex,
+                  ptrdiff_t zindex){
+      rnumber tmp[3];
+      for (int cc=0; cc<3; cc++)
+          tmp[cc] = (this->u->rval(rindex,(cc+1)%3)*this->rvorticity->rval(rindex,(cc+2)%3) -
+                     this->u->rval(rindex,(cc+2)%3)*this->rvorticity->rval(rindex,(cc+1)%3));