From 61f3bcc30816c82152ab9c08b7bd29ea730f89d8 Mon Sep 17 00:00:00 2001
From: sniklas142 <niklas.schnierstein@googlemail.com>
Date: Fri, 4 Oct 2019 09:20:10 +0200
Subject: [PATCH] added vorticity field to ou process
---
cpp/full_code/ornstein_uhlenbeck_process.cpp | 34 ++++++++++++++++++++
cpp/full_code/ornstein_uhlenbeck_process.hpp | 3 ++
2 files changed, 37 insertions(+)
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index 6e6893fe..00479dd2 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -28,6 +28,12 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
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->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);
@@ -48,6 +54,7 @@ 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;
}
@@ -81,6 +88,9 @@ void ornstein_uhlenbeck_process<rnumber,be>::step(
}
});
+
+ this->calc_ou_vorticity();
+
// this->ou_field->symmetrize(); ???
}
@@ -244,6 +254,30 @@ void ornstein_uhlenbeck_process<rnumber,be>::add_to_field_gaussian(
);
}
+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>;
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index 6883c965..768025e2 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -22,6 +22,7 @@ class ornstein_uhlenbeck_process{
double ou_energy_amplitude;
field<rnumber,be,THREE> *ou_field;
+ field<rnumber,be,THREE> *ou_field_vort;
field<rnumber,be,THREExTHREE> *B;
// field<rnumber,be,ONE> *number_of_modes;
kspace<be,SMOOTH> *kk;
@@ -61,6 +62,8 @@ class ornstein_uhlenbeck_process{
void add_to_field_replace(
field<rnumber,be,THREE> *src);
+ void calc_ou_vorticity(void);
+
};
--
GitLab