diff --git a/TurTLE/DNS_statistics.py b/TurTLE/DNS_statistics.py
index 00a016f87478f12fdfc50d50bdc992f4749d39af..90353155fb4e552b12b92b4c48cb4a123de294ec 100644
--- a/TurTLE/DNS_statistics.py
+++ b/TurTLE/DNS_statistics.py
@@ -67,13 +67,15 @@ def compute_statistics(
self.statistics['dk'] = pp_file['kspace/dk'][...]
self.statistics['kshell'] = pp_file['kspace/kshell'][...]
self.statistics['nshell'] = pp_file['kspace/nshell'][...]
+ else:
+ return None
else:
self.read_parameters()
with self.get_data_file() as data_file:
if 'moments' not in data_file['statistics'].keys():
return None
- iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
- self.parameters['niter_stat']-1),
+ iter0 = min(((data_file['statistics/moments/velocity'].shape[0]-1) *
+ self.parameters['niter_stat']),
iter0)
if type(iter1) == type(None):
iter1 = data_file['iteration'][...]
diff --git a/TurTLE/test/test_turtle_NSE.py b/TurTLE/test/test_turtle_NSE.py
index 9678b9bb205f07eedbb715d62203514b77d1e3b0..f3b632293b1b8c8f92fcb07906f18ea96c74800a 100644
--- a/TurTLE/test/test_turtle_NSE.py
+++ b/TurTLE/test/test_turtle_NSE.py
@@ -88,8 +88,8 @@ def main_basic():
'--njobs', '{0}'.format(njobs),
'--wd', './'])
if not type(plt) == type(None):
- c0 = DNS(simname = 'nsve_for_long_comparison')
- c1 = DNS(simname = 'nse_for_long_comparison')
+ c0 = DNS(simname = 'nsve_for_comparison')
+ c1 = DNS(simname = 'nse_for_comparison')
for cc in [c0, c1]:
cc.compute_statistics()
cc.plot_basic_stats()
diff --git a/cpp/full_code/NSE.cpp b/cpp/full_code/NSE.cpp
index 29acebf71436d04e97567428223ec7bb818a86b7..a16fad01a21f89fa477ade3f764ec41d677fd267 100644
--- a/cpp/full_code/NSE.cpp
+++ b/cpp/full_code/NSE.cpp
@@ -134,8 +134,8 @@ int NSE<rnumber>::step(void)
TIMEZONE("NSE::step");
this->iteration++;
//return this->Euler_step();
+ //return this->Heun_step();
return this->ShuOsher();
- return EXIT_SUCCESS;
}
template <typename rnumber>
@@ -267,7 +267,7 @@ int NSE<rnumber>::add_nonlinear_term(
* \f]
* The result is given in Fourier representation.
* */
- const int sign_array[2] = {1, -1};
+ const std::array<int, 2> sign_array = {1, -1};
/* copy velocity */
this->tmp_vec_field0->real_space_representation = false;
@@ -277,7 +277,7 @@ int NSE<rnumber>::add_nonlinear_term(
this->tmp_vec_field0->ift();
/* compute uu */
- /* store 11 22 33 components */
+ /* store 00 11 22 components */
this->tmp_vec_field1->real_space_representation = true;
this->tmp_vec_field1->RLOOP(
[&](const ptrdiff_t rindex,
@@ -289,7 +289,7 @@ int NSE<rnumber>::add_nonlinear_term(
this->tmp_vec_field0->rval(rindex,cc) *
this->tmp_vec_field0->rval(rindex,cc) ) / this->tmp_vec_field1->npoints;
});
- /* take 11 22 33 components to Fourier representation */
+ /* take 00 11 22 components to Fourier representation */
this->tmp_vec_field1->dft();
this->kk->template dealias<rnumber, THREE>(this->tmp_vec_field1->get_cdata());
@@ -309,7 +309,7 @@ int NSE<rnumber>::add_nonlinear_term(
}
});
- /* store 12 23 31 components */
+ /* store 01 12 20 components */
this->tmp_vec_field1->real_space_representation = true;
this->tmp_vec_field1->RLOOP(
[&](const ptrdiff_t rindex,
@@ -321,7 +321,7 @@ int NSE<rnumber>::add_nonlinear_term(
this->tmp_vec_field0->rval(rindex,cc) *
this->tmp_vec_field0->rval(rindex,(cc+1)%3) ) / this->tmp_vec_field1->npoints;
});
- /* take 12 23 31 components to Fourier representation */
+ /* take 01 12 20 components to Fourier representation */
this->tmp_vec_field1->dft();
this->kk->template dealias<rnumber, THREE>(this->tmp_vec_field1->get_cdata());
@@ -344,6 +344,7 @@ int NSE<rnumber>::add_nonlinear_term(
}
});
+ this->kk->template force_divfree<rnumber>(dst->get_cdata());
dst->symmetrize();
/* done */
@@ -396,7 +397,7 @@ int NSE<rnumber>::Euler_step(void)
const ptrdiff_t yindex,
const ptrdiff_t zindex,
const double k2){
- const double factor = exp(-this->nu*k2 * this->dt);
+ const double factor = std::exp(-this->nu*k2 * this->dt);
for (int cc = 0; cc < 3; cc++)
for (int ccc = 0; ccc < 2; ccc++)
{
@@ -412,6 +413,100 @@ int NSE<rnumber>::Euler_step(void)
return EXIT_SUCCESS;
}
+template <typename rnumber>
+int NSE<rnumber>::Heun_step(void)
+{
+ // technically this is not the Heun method, but rather a
+ // composition between the exact diffusion solution and the Heun method
+ // for the nonlinear term
+ TIMEZONE("NSE::Heun_step");
+
+ // temporary field
+ field<rnumber, FFTW, THREE> *acc1 = new field<rnumber, FFTW, THREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+ field<rnumber, FFTW, THREE> *acc2 = new field<rnumber, FFTW, THREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+
+ this->kk->CLOOP_K2(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ const double factor = std::exp(-0.5*this->nu*k2 * this->dt);
+ for (int cc = 0; cc < 3; cc++)
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ this->velocity->cval(cindex, cc, ccc) *= factor;
+ }
+ });
+ *acc1 = rnumber(0);
+ this->add_nonlinear_term(
+ this->velocity,
+ acc1);
+ this->add_forcing_term(this->velocity, acc1);
+ *(this->tmp_velocity1) = *(this->velocity);
+ this->kk->CLOOP(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (int cc = 0; cc < 3; cc++)
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ this->tmp_velocity1->cval(cindex, cc, ccc) +=
+ this->dt*acc1->cval(cindex, cc, ccc);
+ }
+ });
+ this->kk->template force_divfree<rnumber>(this->tmp_velocity1->get_cdata());
+ this->velocity->symmetrize();
+ this->tmp_velocity1->symmetrize();
+
+ *acc2 = rnumber(0);
+ this->add_nonlinear_term(
+ this->tmp_velocity1,
+ acc2);
+ this->add_forcing_term(this->tmp_velocity1, acc2);
+ this->kk->CLOOP(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (int cc = 0; cc < 3; cc++)
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ this->velocity->cval(cindex, cc, ccc) +=
+ 0.5*this->dt*(acc1->cval(cindex, cc, ccc)
+ + acc2->cval(cindex, cc, ccc));
+ }
+ });
+ this->kk->template force_divfree<rnumber>(this->velocity->get_cdata());
+ this->kk->CLOOP_K2(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ const double factor = std::exp(-0.5*this->nu*k2 * this->dt);
+ for (int cc = 0; cc < 3; cc++)
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ this->velocity->cval(cindex, cc, ccc) *= factor;
+ }
+ });
+ this->velocity->symmetrize();
+
+ // delete temporary field
+ delete acc1;
+ delete acc2;
+ return EXIT_SUCCESS;
+}
+
+
/** \brief Time step with Shu Osher method
*
* The Navier Stokes equations are integrated with a
@@ -438,71 +533,88 @@ int NSE<rnumber>::Euler_step(void)
template <typename rnumber>
int NSE<rnumber>::ShuOsher(void)
{
- *this->tmp_velocity1 = rnumber(0);
- this->add_nonlinear_term(this->velocity, this->tmp_velocity1);
- this->add_forcing_term(this->velocity, this->tmp_velocity1);
+ TIMEZONE("NSE::ShuOsher");
+
+ // temporary field
+ field<rnumber, FFTW, THREE> *acc = new field<rnumber, FFTW, THREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+
+ *acc = rnumber(0);
+ this->add_nonlinear_term(
+ this->velocity,
+ acc);
+ this->add_forcing_term(this->velocity, acc);
+ *(this->tmp_velocity1) = rnumber(0.0);
this->kk->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
const ptrdiff_t yindex,
const ptrdiff_t zindex,
const double k2){
- const double factor = exp(-this->nu*k2*this->dt);
+ const double factor1 = std::exp(-this->nu * k2 * this->dt);
for (int cc = 0; cc < 3; cc++)
for (int ccc = 0; ccc < 2; ccc++)
{
- this->tmp_velocity1->cval(cindex, cc, ccc) = factor*(
- this->velocity->cval(cindex, cc, ccc) +
- this->dt*this->tmp_velocity1->cval(cindex, cc, ccc));
- }});
- this->impose_forcing(this->tmp_velocity1);
- this->kk->template force_divfree<rnumber>(this->tmp_velocity1->get_cdata());
+ this->tmp_velocity1->cval(cindex, cc, ccc) =
+ (this->velocity->cval(cindex, cc, ccc)
+ + this->dt*acc->cval(cindex, cc, ccc))*factor1;
+ }
+ });
+ this->tmp_velocity1->symmetrize();
- *this->tmp_velocity2 = rnumber(0);
- this->add_nonlinear_term(this->tmp_velocity1, this->tmp_velocity2);
- this->add_forcing_term(this->tmp_velocity1, this->tmp_velocity2);
+ *acc = rnumber(0);
+ this->add_nonlinear_term(
+ this->tmp_velocity1,
+ acc);
+ this->add_forcing_term(this->tmp_velocity1, acc);
+ *(this->tmp_velocity2) = rnumber(0.0);
this->kk->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
const ptrdiff_t yindex,
const ptrdiff_t zindex,
const double k2){
- const double factor0 = exp(-0.5*this->nu*k2*this->dt)*0.75;
- const double factor1 = exp( 0.5*this->nu*k2*this->dt)*0.25;
+ const double factor0 = double(3)*std::exp(-0.5*this->nu * k2 * this->dt)/4;
+ const double factor1 = std::exp(0.5*this->nu * k2 * this->dt)/4;
for (int cc = 0; cc < 3; cc++)
for (int ccc = 0; ccc < 2; ccc++)
{
this->tmp_velocity2->cval(cindex, cc, ccc) = (
- this->velocity->cval(cindex, cc, ccc)*factor0 +
- (this->tmp_velocity1->cval(cindex, cc, ccc) +
- this->dt*this->tmp_velocity2->cval(cindex, cc, ccc))*factor1);
- }});
- this->impose_forcing(this->tmp_velocity2);
- this->kk->template force_divfree<rnumber>(this->tmp_velocity2->get_cdata());
-
- *this->tmp_velocity1 = rnumber(0);
- this->add_nonlinear_term(this->tmp_velocity2, this->tmp_velocity1);
- this->add_forcing_term(this->tmp_velocity2, this->tmp_velocity1);
+ this->velocity->cval(cindex, cc, ccc)*factor0
+ + (this->tmp_velocity1->cval(cindex, cc, ccc)
+ + this->dt*acc->cval(cindex, cc, ccc))*factor1);
+ }
+ });
+ this->tmp_velocity2->symmetrize();
+
+ *acc = rnumber(0);
+ this->add_nonlinear_term(
+ this->tmp_velocity2,
+ acc);
+ this->add_forcing_term(this->tmp_velocity2, acc);
this->kk->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
const ptrdiff_t yindex,
const ptrdiff_t zindex,
const double k2){
- const double factor0 = exp(- this->nu*k2*this->dt) / 3;
- const double factor1 = 2*exp(-0.5*this->nu*k2*this->dt) / 3;
+ const double factor0 = std::exp(-this->nu * k2 * this->dt)/3;
+ const double factor1 = double(2)*std::exp(-0.5*this->nu * k2 * this->dt)/3;
for (int cc = 0; cc < 3; cc++)
for (int ccc = 0; ccc < 2; ccc++)
{
this->velocity->cval(cindex, cc, ccc) = (
- this->velocity->cval(cindex, cc, ccc)*factor0 +
- (this->tmp_velocity2->cval(cindex, cc, ccc) +
- this->dt*this->tmp_velocity1->cval(cindex, cc, ccc))*factor1);
- }});
- this->impose_forcing(this->velocity);
- this->kk->template force_divfree<rnumber>(this->velocity->get_cdata());
-
+ this->velocity->cval(cindex, cc, ccc)*factor0
+ + (this->tmp_velocity2->cval(cindex, cc, ccc)
+ + this->dt*acc->cval(cindex, cc, ccc))*factor1);
+ }
+ });
this->velocity->symmetrize();
+
+ // delete temporary field
+ delete acc;
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/NSE.hpp b/cpp/full_code/NSE.hpp
index 8473425ee79380df04a004fe3a5d8dad4f403aaa..d0d674d29a90f1d999756d75071066a412ab772f 100644
--- a/cpp/full_code/NSE.hpp
+++ b/cpp/full_code/NSE.hpp
@@ -83,6 +83,7 @@ class NSE: public direct_numerical_simulation
int finalize(void);
int Euler_step(void);
+ int Heun_step(void);
int ShuOsher(void);
int add_field_band(
diff --git a/documentation/cpp/cpp_config b/documentation/cpp/cpp_config
index e5fa2a6d5b6ab0e54cba96d3cb8e1f94805342af..9b7c7915b6c503a1cb40d66bbb5bae4f8d824227 100644
--- a/documentation/cpp/cpp_config
+++ b/documentation/cpp/cpp_config
@@ -1540,7 +1540,7 @@ FORMULA_TRANSPARENT = YES
# The default value is: NO.
# This tag requires that the tag GENERATE_HTML is set to YES.
-USE_MATHJAX = NO
+USE_MATHJAX = YES
# When MathJax is enabled you can set the default output format to be used for
# the MathJax output. See the MathJax site (see: