diff --git a/.gitignore b/.gitignore
index e18d19c5b68084e3bef12fa670095f1b73b7af45..ff58edaef9b1c0512d146c6ab80654b6531263ce 100644
--- a/.gitignore
+++ b/.gitignore
@@ -11,4 +11,5 @@ bfps/install_info.pickle
bfps/lib*.a
obj
build*
-venv*
\ No newline at end of file
+venv*
+cpp/turtle_config.hpp
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index bdd2cc323a4e31dfc968ae74286996d519ec3492..87cf869306997fa294d82b030ee420a2b6a1da6b 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -1447,6 +1447,7 @@ class DNS(_code):
c.set_ylim(bottom = 0)
c.set_xlabel('$t / \\tau_K$')
c.legend(loc = 'best')
+ fig.suptitle(self.simname)
fig.tight_layout()
fig.savefig(os.path.join(fig_dir, 'stats' + fig_suffix + '.' + fig_format))
return None
diff --git a/TurTLE/test/test_fftw.py b/TurTLE/test/test_fftw.py
index bb35935a4bad1660fd15893a162efd24ef7e3867..b0a787079608add7349ea8c6c37d1a377e4c2a1c 100644
--- a/TurTLE/test/test_fftw.py
+++ b/TurTLE/test/test_fftw.py
@@ -3,6 +3,7 @@
import numpy as np
import h5py
import sys
+import argparse
import TurTLE
from TurTLE import TEST
@@ -18,8 +19,12 @@ except:
cmocean = None
def main():
+ parser = argparse.ArgumentParser('Compare FFTW with `numpy.fft`.')
+ parser.add_argument('-p', '--plot', dest = 'plot_on', action = 'store_true')
+ opt = parser.parse_args()
+
niterations = 10
- nlist = [16, 32, 48, 24, 64, 12]
+ nlist = [8, 24, 32, 16, 12]
for ii in range(len(nlist)):
c = TEST()
c.launch(
@@ -54,9 +59,15 @@ def main():
print('found maximum difference {0} between field0 and field1 at {1}'.format(diff[ii[0], ii[1], ii[2]], ii))
assert(np.max(np.abs(diff)) < 1e-5)
+
+ def L2norm_Fourier(ff):
+ return np.sqrt(
+ np.sum(ff[:, :, 0 ] * np.conj(ff[:, :, 0 ]))
+ + 2*np.sum(ff[:, :, 1:] * np.conj(ff[:, :, 1:])))
+
# now compare numpy with fftw
np_field1_real = np.fft.irfftn(field1_complex, axes = (0, 1, 2)).transpose(1, 0, 2, 3)
- L2normr = np.sqrt(np.mean(np.sum(field1_real**2, axis = 3)))
+ L2normr = np.sqrt(np.mean(np.sum( field1_real**2, axis = 3)))
np_L2normr = np.sqrt(np.mean(np.sum(np_field1_real**2, axis = 3)))
err = np.max(np.abs(field1_real - np_field1_real*npoints)) / L2normr
print('found maximum difference {0} between field1_real and np_field1_real'.format(err))
@@ -64,15 +75,15 @@ def main():
np_field1_complex = np.fft.rfftn(field1_real.transpose(1, 0, 2, 3), axes = (0, 1, 2)) / npoints
- L2norm0 = np.sqrt(np.sum(np.abs(field1_complex[:, :, 0])**2) + 2*np.sum(np.abs(field1_complex[:, :, 1:])**2))
- L2norm1 = np.sqrt(np.sum(np.abs(np_field1_complex[:, :, 0])**2) + 2*np.sum(np.abs(np_field1_complex[:, :, 1:])**2))
+ L2norm0 = L2norm_Fourier(field1_complex)
+ L2norm1 = L2norm_Fourier(np_field1_complex)
err = np.max(np.abs(np_field1_complex - field1_complex)) / L2norm0
assert(err < 1e-5)
err = abs(L2normr - L2norm0) / L2norm0
assert(err < 1e-5)
- if not type(plt) == type(None):
+ if ((not type(plt) == type(None)) and (opt.plot_on)):
f = plt.figure()
a = f.add_subplot(221)
a.imshow(np.log(np.abs(np_field1_complex[:, :, 0, 0])), interpolation = 'nearest')
diff --git a/TurTLE/test/test_fluid_convergence.py b/TurTLE/test/test_fluid_convergence.py
new file mode 100644
index 0000000000000000000000000000000000000000..d79cd0039837232c33c36dc49c5151902780ef5d
--- /dev/null
+++ b/TurTLE/test/test_fluid_convergence.py
@@ -0,0 +1,151 @@
+#! /usr/bin/env python
+#######################################################################
+# #
+# Copyright 2021 Max Planck Institute #
+# for Dynamics and Self-Organization #
+# #
+# This file is part of TurTLE. #
+# #
+# 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, #
+# 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/> #
+# #
+# Contact: Cristian.Lalescu@mpcdf.mpg.de #
+# #
+#######################################################################
+
+
+
+import os
+import numpy as np
+import h5py
+import sys
+import argparse
+
+import TurTLE
+from TurTLE import DNS, PP
+
+
+try:
+ import matplotlib.pyplot as plt
+except:
+ plt = None
+
+def main_NSE():
+ niterations = 32
+ c0 = DNS()
+ c0.launch(
+ ['NSE',
+ '-n', '64',
+ '--simname', 'nse0',
+ '--np', '16',
+ '--ntpp', '1',
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(niterations//16),
+ '--niter_stat', '1',
+ '--checkpoints_per_file', '{0}'.format(3),
+ '--wd', './'])
+ c0.compute_statistics()
+ c0.plot_basic_stats()
+ #cc = PP()
+ #cc.launch(
+ # ['get_velocity',
+ # '--simname', 'nsve_for_long_comparison',
+ # '--np', '16',
+ # '--ntpp', '1',
+ # '--iter0', '0',
+ # '--iter1', '{0}'.format(niterations)])
+ #c1 = DNS()
+ #c1.launch(
+ # ['NSE',
+ # '-n', '64',
+ # '--src-simname', 'nsve_for_long_comparison',
+ # '--src-wd', './',
+ # '--src-iteration', '0',
+ # '--simname', 'nse_for_long_comparison',
+ # '--np', '16',
+ # '--ntpp', '1',
+ # '--fftw_plan_rigor', 'FFTW_PATIENT',
+ # '--niter_todo', '{0}'.format(niterations),
+ # '--niter_out', '{0}'.format(niterations),
+ # '--niter_stat', '1',
+ # '--forcing_type', 'linear',
+ # '--checkpoints_per_file', '{0}'.format(3),
+ # '--njobs', '{0}'.format(njobs),
+ # '--wd', './'])
+ #def get_err(a, b):
+ # return np.abs(a - b) / (0.5 * (a + b))
+ #if not type(plt) == type(None):
+ # c0 = DNS(simname = 'nsve_for_long_comparison')
+ # c1 = DNS(simname = 'nse_for_long_comparison')
+ # for cc in [c0, c1]:
+ # cc.compute_statistics()
+ # cc.plot_basic_stats()
+ # ##########
+ # fig = plt.figure()
+ # a = fig.add_subplot(211)
+ # a.plot(c0.statistics['energy(t)'],
+ # label = 'NSVE')
+ # a.plot(c1.statistics['energy(t)'],
+ # label = 'NSE')
+ # a.legend(loc = 'best')
+ # a = fig.add_subplot(212)
+ # a.plot(c0.statistics['kshell'],
+ # c0.statistics['energy(k)'],
+ # label = 'NSVE')
+ # a.plot(c1.statistics['kshell'],
+ # c1.statistics['energy(k)'],
+ # label = 'NSE')
+ # a.set_xscale('log')
+ # a.set_yscale('log')
+ # a.set_ylim(1e-7, 10)
+ # fig.tight_layout()
+ # fig.savefig('comparison_energy_long.pdf')
+ # plt.close(fig)
+ # ##########
+ # fig = plt.figure()
+ # a = fig.add_subplot(211)
+ # a.plot(c0.statistics['vel_max(t)'],
+ # label = 'NSVE')
+ # a.plot(c1.statistics['vel_max(t)'],
+ # label = 'NSE')
+ # a.legend(loc = 'best')
+ # a = fig.add_subplot(212)
+ # err = get_err(
+ # c0.statistics['vel_max(t)'],
+ # c1.statistics['vel_max(t)'])
+ # a.plot(c0.statistics['t'] / c0.statistics['tauK'],
+ # err,
+ # label = 'velocity')
+ # vort0_max = c0.get_data_file()['statistics/moments/vorticity'][:, 9, 3]
+ # vort1_max = c1.get_data_file()['statistics/moments/vorticity'][:, 9, 3]
+ # a.plot(c0.statistics['t'] / c0.statistics['tauK'],
+ # get_err(vort0_max, vort1_max),
+ # label = 'vorticity')
+ # a.legend(loc = 'best')
+ # a.set_yscale('log')
+ # a.set_ylim(bottom = 1e-6)
+ # fig.tight_layout()
+ # fig.savefig('comparison_velmax_long.pdf')
+ # plt.close(fig)
+ return None
+
+if __name__ == '__main__':
+ parser = argparse.ArgumentParser('')
+ parser.add_argument('--NSE', dest = 'NSE', action = 'store_true')
+ opt = parser.parse_args()
+ if opt.NSE:
+ main_NSE()
+ else:
+ pass
+
diff --git a/TurTLE/test/test_turtle_NSE.py b/TurTLE/test/test_turtle_NSE.py
index 6fe5633dfa35e2a5ed1acf4ebc5b4000bb8792d2..62ff60fcf398d7bd97de3843c9ee7ebb02732ae4 100644
--- a/TurTLE/test/test_turtle_NSE.py
+++ b/TurTLE/test/test_turtle_NSE.py
@@ -29,11 +29,11 @@ import os
import numpy as np
import h5py
import sys
+import argparse
import TurTLE
from TurTLE import DNS, PP
-from TurTLE_addons import NSReader
try:
import matplotlib.pyplot as plt
@@ -61,8 +61,7 @@ def main_basic():
'--niter_stat', '1',
'--checkpoints_per_file', '{0}'.format(3),
'--njobs', '{0}'.format(njobs),
- '--wd', './'] +
- sys.argv[1:])
+ '--wd', './'])
cc = PP()
cc.launch(
['get_velocity',
@@ -70,8 +69,7 @@ def main_basic():
'--np', '4',
'--ntpp', '2',
'--iter0', '0',
- '--iter1', '64'] +
- sys.argv[1:])
+ '--iter1', '64'])
c1 = DNS()
c1.launch(
['NSE',
@@ -89,8 +87,13 @@ def main_basic():
'--forcing_type', 'linear',
'--checkpoints_per_file', '{0}'.format(3),
'--njobs', '{0}'.format(njobs),
- '--wd', './'] +
- sys.argv[1:])
+ '--wd', './'])
+ if not type(plt) == type(None):
+ c0 = DNS(simname = 'nsve_for_long_comparison')
+ c1 = DNS(simname = 'nse_for_long_comparison')
+ for cc in [c0, c1]:
+ cc.compute_statistics()
+ cc.plot_basic_stats()
f0 = h5py.File(
'nsve_for_comparison.h5', 'r')
f1 = h5py.File(
@@ -133,8 +136,8 @@ def main_long():
'--src-wd', TurTLE.data_dir,
'--src-iteration', '0',
'--simname', 'nsve_for_long_comparison',
- '--np', '4',
- '--ntpp', '2',
+ '--np', '16',
+ '--ntpp', '1',
'--fftw_plan_rigor', 'FFTW_PATIENT',
'--niter_todo', '{0}'.format(niterations),
'--niter_out', '{0}'.format(niterations),
@@ -142,17 +145,15 @@ def main_long():
'--niter_stat', '1',
'--checkpoints_per_file', '{0}'.format(3),
'--njobs', '{0}'.format(njobs),
- '--wd', './'] +
- sys.argv[1:])
+ '--wd', './'])
cc = PP()
cc.launch(
['get_velocity',
'--simname', 'nsve_for_long_comparison',
- '--np', '4',
- '--ntpp', '2',
+ '--np', '16',
+ '--ntpp', '1',
'--iter0', '0',
- '--iter1', '64'] +
- sys.argv[1:])
+ '--iter1', '{0}'.format(niterations)])
c1 = DNS()
c1.launch(
['NSE',
@@ -161,8 +162,8 @@ def main_long():
'--src-wd', './',
'--src-iteration', '0',
'--simname', 'nse_for_long_comparison',
- '--np', '4',
- '--ntpp', '2',
+ '--np', '16',
+ '--ntpp', '1',
'--fftw_plan_rigor', 'FFTW_PATIENT',
'--niter_todo', '{0}'.format(niterations),
'--niter_out', '{0}'.format(niterations),
@@ -170,19 +171,22 @@ def main_long():
'--forcing_type', 'linear',
'--checkpoints_per_file', '{0}'.format(3),
'--njobs', '{0}'.format(njobs),
- '--wd', './'] +
- sys.argv[1:])
- c0 = NSReader(simname = 'nsve_for_long_comparison')
- c1 = NSReader(simname = 'nse_for_long_comparison')
- c0.do_plots()
- c1.do_plots()
+ '--wd', './'])
+ def get_err(a, b):
+ return np.abs(a - b) / (0.5 * (a + b))
if not type(plt) == type(None):
+ c0 = DNS(simname = 'nsve_for_long_comparison')
+ c1 = DNS(simname = 'nse_for_long_comparison')
+ for cc in [c0, c1]:
+ cc.compute_statistics()
+ cc.plot_basic_stats()
+ ##########
fig = plt.figure()
a = fig.add_subplot(211)
a.plot(c0.statistics['energy(t)'],
label = 'NSVE')
a.plot(c1.statistics['energy(t)'],
- label = 'NSVE')
+ label = 'NSE')
a.legend(loc = 'best')
a = fig.add_subplot(212)
a.plot(c0.statistics['kshell'],
@@ -197,9 +201,40 @@ def main_long():
fig.tight_layout()
fig.savefig('comparison_energy_long.pdf')
plt.close(fig)
+ ##########
+ fig = plt.figure()
+ a = fig.add_subplot(211)
+ a.plot(c0.statistics['vel_max(t)'],
+ label = 'NSVE')
+ a.plot(c1.statistics['vel_max(t)'],
+ label = 'NSE')
+ a.legend(loc = 'best')
+ a = fig.add_subplot(212)
+ err = get_err(
+ c0.statistics['vel_max(t)'],
+ c1.statistics['vel_max(t)'])
+ a.plot(c0.statistics['t'] / c0.statistics['tauK'],
+ err,
+ label = 'velocity')
+ vort0_max = c0.get_data_file()['statistics/moments/vorticity'][:, 9, 3]
+ vort1_max = c1.get_data_file()['statistics/moments/vorticity'][:, 9, 3]
+ a.plot(c0.statistics['t'] / c0.statistics['tauK'],
+ get_err(vort0_max, vort1_max),
+ label = 'vorticity')
+ a.legend(loc = 'best')
+ a.set_yscale('log')
+ a.set_ylim(bottom = 1e-6)
+ fig.tight_layout()
+ fig.savefig('comparison_velmax_long.pdf')
+ plt.close(fig)
return None
if __name__ == '__main__':
- main_basic()
- #main_long()
+ parser = argparse.ArgumentParser('compare NSE with NSVE')
+ parser.add_argument('--long', dest = 'long', action = 'store_true')
+ opt = parser.parse_args()
+ if opt.long:
+ main_long()
+ else:
+ main_basic()
diff --git a/cpp/full_code/NSE.cpp b/cpp/full_code/NSE.cpp
index 1ed1339fa0d089f8a05ace5693d23920e07b4c3f..52bee4186629bb564c39c895c5e466d90e8fcb28 100644
--- a/cpp/full_code/NSE.cpp
+++ b/cpp/full_code/NSE.cpp
@@ -34,8 +34,7 @@ int NSE<rnumber>::initialize(void)
TIMEZONE("NSE::initialize");
this->read_iteration();
this->read_parameters();
- if (this->myrank == 0)
- {
+ if (this->myrank == 0) {
// set caching parameters
hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
@@ -47,11 +46,11 @@ int NSE<rnumber>::initialize(void)
fapl);
}
int data_file_problem;
- if (this->myrank == 0)
+ if (this->myrank == 0) {
data_file_problem = this->grow_file_datasets();
+ }
MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
- if (data_file_problem > 0)
- {
+ if (data_file_problem > 0) {
std::cerr <<
data_file_problem <<
" problems growing file datasets.\ntrying to exit now." <<
@@ -61,53 +60,82 @@ int NSE<rnumber>::initialize(void)
/* allocate vector fields */
this->velocity = new field<rnumber, FFTW, THREE>(
- this->nx, this->ny, this->nz,
+ this->nx,
+ this->ny,
+ this->nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->tmp_velocity1 = new field<rnumber, FFTW, THREE>(
- this->nx, this->ny, this->nz,
+ this->nx,
+ this->ny,
+ this->nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->tmp_velocity2 = new field<rnumber, FFTW, THREE>(
- this->nx, this->ny, this->nz,
+ this->nx,
+ this->ny,
+ this->nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->tmp_vec_field0 = new field<rnumber, FFTW, THREE>(
- this->nx, this->ny, this->nz,
+ this->nx,
+ this->ny,
+ this->nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
this->tmp_vec_field1 = new field<rnumber, FFTW, THREE>(
- this->nx, this->ny, this->nz,
+ this->nx,
+ this->ny,
+ this->nz,
this->comm,
fftw_planner_string_to_flag[this->fftw_plan_rigor]);
/* initialize kspace */
this->kk = new kspace<FFTW, SMOOTH>(
this->velocity->clayout, this->dkx, this->dky, this->dkz);
+ if (this->myrank == 0 && this->iteration == 0) {
+ this->kk->store(stat_file);
+ }
- {
- this->velocity->real_space_representation = false;
- this->velocity->io(
- this->get_current_fname(),
- "velocity",
- this->iteration,
- true);
- TIMEZONE("NSE::initialize::read_initial_condition");
+ /* initialize field */
+ this->velocity->real_space_representation = false;
+ if ((this->iteration == 0)
+ && (this->field_random_seed != 0)) {
+ TIMEZONE("NSE::initialize::generate_initial_condition");
+ // generate initial condition
+ make_gaussian_random_field(
+ this->kk,
+ this->velocity,
+ this->field_random_seed,
+ this->injection_rate,
+ 1.0, // Lint
+ 1.5 / this->kk->kM, // etaK
+ 6.78,
+ 0.40,
+ 3./2.);
this->kk->template low_pass<rnumber, THREE>(
this->velocity->get_cdata(),
this->kk->kM);
- this->kk->template force_divfree<rnumber>(
+ this->kk->template project_divfree<rnumber>(
this->velocity->get_cdata());
this->velocity->symmetrize();
- MPI_Barrier(this->comm);
+ this->velocity->io(
+ this->get_current_fname(),
+ "velocity",
+ this->iteration,
+ false);
+ } else {
+ TIMEZONE("NSE::initialize::read_initial_condition");
+ this->velocity->io(
+ this->get_current_fname(),
+ "velocity",
+ this->iteration,
+ true);
}
-
- if (this->myrank == 0 && this->iteration == 0)
- this->kk->store(stat_file);
return EXIT_SUCCESS;
}
@@ -557,6 +585,15 @@ int NSE<rnumber>::do_stats()
this->iteration / niter_stat,
this->max_velocity_estimate/std::sqrt(3));
+ double CFL_velocity = tmp_vec_field0->compute_CFL_velocity();
+ if (this->myrank == 0)
+ hdf5_tools::update_time_series_with_single_rank<double>(
+ stat_group,
+ "CFL_velocity",
+ this->iteration / niter_stat,
+ 1,
+ &CFL_velocity);
+
compute_curl(this->kk,
this->velocity,
this->tmp_vec_field0);
@@ -582,6 +619,7 @@ int NSE<rnumber>::read_parameters(void)
this->nu = hdf5_tools::read_value<double>(parameter_file, "parameters/nu");
this->dt = hdf5_tools::read_value<double>(parameter_file, "parameters/dt");
this->fmode = hdf5_tools::read_value<int>(parameter_file, "parameters/fmode");
+ this->field_random_seed = hdf5_tools::read_value<int>(parameter_file, "parameters/field_random_seed");
this->famplitude = hdf5_tools::read_value<double>(parameter_file, "parameters/famplitude");
this->friction_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/friction_coefficient");
this->injection_rate = hdf5_tools::read_value<double>(parameter_file, "parameters/injection_rate");
diff --git a/cpp/full_code/NSE.hpp b/cpp/full_code/NSE.hpp
index ef0813e74aa3a68980468d0d0a377015eda0edb2..4ff447e1afe2fd5900e8a7630803b903a771ed46 100644
--- a/cpp/full_code/NSE.hpp
+++ b/cpp/full_code/NSE.hpp
@@ -55,6 +55,7 @@ class NSE: public direct_numerical_simulation
double energy;
double injection_rate;
int fmode;
+ int field_random_seed;
std::string forcing_type;
int histogram_bins;
double max_velocity_estimate;
diff --git a/cpp/full_code/NSVE.cpp b/cpp/full_code/NSVE.cpp
index 24fc890fef55d50d07d2b2e0439f276b86e54913..9f323c587579ed301a717c404b8503f27bcf52b4 100644
--- a/cpp/full_code/NSVE.cpp
+++ b/cpp/full_code/NSVE.cpp
@@ -112,6 +112,7 @@ int NSVE<rnumber>::initialize(void)
compute_curl(this->fs->kk, this->fs->cvelocity, this->fs->cvorticity);
this->fs->kk->template project_divfree<rnumber>(
this->fs->cvorticity->get_cdata());
+ this->fs->cvorticity->symmetrize();
this->fs->io_checkpoint(false);
} else {
this->fs->io_checkpoint();
diff --git a/cpp/full_code/code_base.hpp b/cpp/full_code/code_base.hpp
index 2886447dc7b3a5d8658d1a8f08ec8d994318cc59..3cce833bab8490ad365ed8dccdbf8dfdc35b39f4 100644
--- a/cpp/full_code/code_base.hpp
+++ b/cpp/full_code/code_base.hpp
@@ -109,6 +109,12 @@ class code_base
return EXIT_SUCCESS;
}
+ inline MPI_Comm get_communicator() const {return this->comm;}
+
+ inline int get_nx() const {return this->nx;}
+ inline int get_ny() const {return this->ny;}
+ inline int get_nz() const {return this->nz;}
+
/** Reads parameters
* \warning This method should ensure the parameter file is closed by
* all MPI processes when finished.
diff --git a/cpp/full_code/shared_array_test.cpp b/cpp/full_code/shared_array_test.cpp
index e1be24be04b0b1bf26beda67c45a7c80ee3122e0..1791adec63a9b07e2f4c1143d747779a9379a40f 100644
--- a/cpp/full_code/shared_array_test.cpp
+++ b/cpp/full_code/shared_array_test.cpp
@@ -57,9 +57,9 @@ template <typename rnumber>
int shared_array_test<rnumber>::do_work(void)
{
TIMEZONE("shared_array_test::do_work");
- const int array_size = 1000;
+ const int array_size = 500;
double value = 0.0;
- for (auto ii = 0; ii < 2000; ii++)
+ for (auto ii = 0; ii < 1000; ii++)
{
shared_array<double> local_value(
array_size,
diff --git a/examples/field_convergence/EDNS.py b/examples/field_convergence/EDNS.py
new file mode 100644
index 0000000000000000000000000000000000000000..de8e5056ac5ac62763665cf4e0d6d676c87ea8a3
--- /dev/null
+++ b/examples/field_convergence/EDNS.py
@@ -0,0 +1,296 @@
+import os
+import sys
+import getpass
+import subprocess
+
+import numpy as np
+import h5py
+import matplotlib.pyplot as plt
+
+import TurTLE
+import TurTLE._base
+import TurTLE.DNS
+from TurTLE._code import _code
+
+username = getpass.getuser()
+homefolder = os.path.expanduser('~')
+cpp_location = None
+
+import pathlib
+cpp_location = pathlib.Path(__file__).parent.resolve()
+
+
+class EDNS(TurTLE.DNS):
+ def __init__(
+ self,
+ **kwargs):
+ TurTLE.DNS.__init__(self, **kwargs)
+ self.list_of_particle_codes = []
+ self.extra_parameters = {
+ 'kraichnan_field' : {},
+ 'NSVE_Stokes_particles' : {},
+ 'NSVEp_extra_sampling' : {}}
+ self.check_current_vorticity_exists = False
+ return None
+ def write_src(
+ self):
+ self.version_message = (
+ '/***********************************************************************\n' +
+ '* this code automatically generated by TurTLE\n' +
+ '* version {0}\n'.format(TurTLE.__version__) +
+ '***********************************************************************/\n\n\n')
+ self.include_list = [
+ '"full_code/main_code.hpp"',
+ '<cmath>']
+ self.main = """
+ int main(int argc, char *argv[])
+ {{
+ bool fpe = (
+ (getenv("BFPS_FPE_OFF") == nullptr) ||
+ (getenv("BFPS_FPE_OFF") != std::string("TRUE")));
+ return main_code< {0} >(argc, argv, fpe);
+ }}
+ """.format(self.dns_type + '<{0}>'.format(self.C_field_dtype))
+ self.includes = '\n'.join(
+ ['#include ' + hh
+ for hh in self.include_list])
+ self.name = 'fluid_error'
+ self.dns_type = 'fluid_error'
+ self.definitions = ''
+ self.definitions += open(
+ os.path.join(
+ cpp_location, 'abstract_fluid_solver.hpp'), 'r').read()
+ self.definitions += open(
+ os.path.join(
+ cpp_location, 'fluid_solver.hpp'), 'r').read()
+ self.definitions += open(
+ os.path.join(
+ cpp_location, 'fluid_error.hpp'), 'r').read()
+ self.definitions += open(
+ os.path.join(
+ cpp_location, 'fluid_error.cpp'), 'r').read()
+ with open(self.name + '.cpp', 'w') as outfile:
+ outfile.write(self.version_message + '\n\n')
+ outfile.write(self.includes + '\n\n')
+ outfile.write(self.definitions + '\n\n')
+ outfile.write(self.main + '\n')
+ return None
+ def generate_default_parameters(self):
+ TurTLE.DNS.generate_default_parameters(self)
+ self.extra_parameters['fluid_error'] = {
+ 'tfs' : 1,
+ 'max_err_estimate' : 1e-5,
+ 'max_err_relative_estimate' : 1e-4,
+ }
+ return None
+
+ def write_par(
+ self,
+ iter0 = 0):
+ assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
+ assert (self.parameters['niter_todo'] % self.parameters['niter_out'] == 0)
+ assert (self.parameters['niter_out'] % self.parameters['niter_stat'] == 0)
+ #os.chdir(os.path.join(self.work_dir, self.simname))
+ _code.write_par(self, iter0 = iter0)
+ with h5py.File(self.get_data_file_name(), 'r+') as ofile:
+ ofile['code_info/exec_name'] = self.name
+ kspace = self.get_kspace()
+ for k in kspace.keys():
+ ofile['kspace/' + k] = kspace[k]
+ nshells = kspace['nshell'].shape[0]
+ kspace = self.get_kspace()
+ nshells = kspace['nshell'].shape[0]
+ vec_stat_datasets = [
+ 'velocity_err1',
+ 'velocity_err2',
+ 'velocity_err1_relative',
+ 'velocity_err2_relative',
+ 'vorticity_err1',
+ 'vorticity_err2',
+ 'vorticity_err1_relative',
+ 'vorticity_err2_relative']
+ scal_stat_datasets = []
+ for k in vec_stat_datasets:
+ time_chunk = 2**20//(8*3*3*nshells)
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+ (1, nshells, 3, 3),
+ chunks = (time_chunk, nshells, 3, 3),
+ maxshape = (None, nshells, 3, 3),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*4*10)
+ time_chunk = max(time_chunk, 1)
+ a = ofile.create_dataset('statistics/moments/' + k,
+ (1, 10, 4),
+ chunks = (time_chunk, 10, 4),
+ maxshape = (None, 10, 4),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/histograms/' + k,
+ (1,
+ self.parameters['histogram_bins'],
+ 4),
+ chunks = (time_chunk,
+ self.parameters['histogram_bins'],
+ 4),
+ maxshape = (None,
+ self.parameters['histogram_bins'],
+ 4),
+ dtype = np.int64)
+ ofile['checkpoint'] = int(0)
+ for member_id in [1, 2, 4]:
+ simname = self.simname + str(member_id)
+ cc = TurTLE.DNS()
+ fluid_type = 'NSVE'
+ if self.parameters['tfs'] == 1:
+ fluid_type = 'NSE'
+ cc.prepare_launch(
+ args = [
+ fluid_type,
+ '--simname', simname,
+ '--wd', self.work_dir])
+ for kk in cc.parameters.keys():
+ cc.parameters[kk] = self.parameters[kk]
+ # member-specific parameters
+ cc.parameters['dt'] = self.parameters['dt'] / member_id
+ cc.parameters['niter_todo'] = self.parameters['niter_todo'] * member_id
+ cc.parameters['niter_out'] = self.parameters['niter_out'] * member_id
+ cc.parameters['niter_stat'] = self.parameters['niter_stat'] * member_id
+ cc.write_par()
+ return None
+ def prepare_launch(
+ self,
+ **kwargs):
+ opt = TurTLE.DNS.prepare_launch(self, **kwargs)
+ return opt
+ def launch(
+ self,
+ args = [],
+ **kwargs):
+ # if no --simname argument is added to args with the same value as self.ensname (or self.simname),
+ # self.simname is overwritten with the DNS default inside of prepare_launch.
+ # to avoid multiple assignments of self.ensname, we enforce the correct naming scheme here.
+ assert('--simname' not in args)
+ args += ['--simname', self.simname]
+
+ opt = self.prepare_launch(args = args)
+ if not os.path.exists(self.get_data_file_name()):
+ self.write_par()
+ self.run(
+ nb_processes = opt.nb_processes,
+ nb_threads_per_process = opt.nb_threads_per_process,
+ njobs = opt.njobs,
+ hours = opt.minutes // 60,
+ minutes = opt.minutes % 60,
+ no_submit = opt.no_submit,
+ no_debug = opt.no_debug)
+ return None
+ def run_sanity_check(
+ self,
+ fname = None):
+ if type(fname) == type(None):
+ fname = self.simname
+ ofile = open(fname + '.tex', 'w')
+ ofile.write('\\documentclass{article}\n')
+ ofile.write('\\usepackage{graphicx}\n')
+ ofile.write('\\usepackage[a4paper,'
+ + 'bindingoffset = 0cm,'
+ + 'left = 1cm,'
+ + 'right = 1cm,'
+ + 'top = 1cm,'
+ + 'bottom = 2cm]{geometry}\n')
+ ofile.write('\\begin{document}\n')
+ members = self.get_member_DNS_list()
+ for cc in members:
+ cc.read_parameters()
+ print(cc.simname)
+ print(cc.parameters)
+ cc.compute_statistics()
+ cc.plot_basic_stats(fig_dir = './figs')
+ ofile.write('\\verb|' + cc.work_dir + '|\n\n')
+ ofile.write('\\verb|' + cc.simname + '|\n\n')
+ ofile.write('forcing type is \\verb|' + cc.parameters['forcing_type'] + '|\n\n')
+ ofile.write('\\includegraphics[width=\\textwidth]{'
+ + './figs/stats_' + cc.simname + '.pdf}\n')
+ ofile.write('\\newpage\n')
+ ofile.write('\\end{document}\n')
+ ofile.close()
+ subprocess.check_call(['pdflatex', fname + '.tex'])
+ return members
+ def get_member_DNS_list(self):
+ """
+ Method called "get_member_DNS_list" to make it obvious that
+ this returns a list of DNS objects.
+ """
+ member_list = [ TurTLE.DNS(
+ simname = '{0}{1}'.format(self.simname, i),
+ work_dir = self.work_dir)
+ for i in [1, 2, 4]]
+ return member_list
+ def read_error_statistics(
+ self):
+ df = self.get_data_file()
+ for ff in ['velocity', 'vorticity']:
+ for ii in [1, 2]:
+ self.statistics['error/moments/{0}_err{1}'.format(ff, ii)] = \
+ df['statistics/moments/{0}_err{1}'.format(ff, ii)][()]
+ self.statistics['error/moments/{0}_err{1}_relative'.format(ff, ii)] = \
+ df['statistics/moments/{0}_err{1}_relative'.format(ff, ii)][()]
+ self.statistics['error/histograms/{0}_err{1}'.format(ff, ii)] = \
+ df['statistics/histograms/{0}_err{1}'.format(ff, ii)][()]
+ self.statistics['error/histograms/{0}_err{1}_relative'.format(ff, ii)] = \
+ df['statistics/histograms/{0}_err{1}_relative'.format(ff, ii)][()]
+ df.close()
+ return None
+ def plot_error_statistics(
+ self,
+ fname = None,
+ err_norm = 1.0):
+ if type(fname) == type(None):
+ fname = self.simname + '_err_stats'
+ for ff in ['velocity', 'vorticity']:
+ f = plt.figure(figsize = (8, 8))
+ a = f.add_subplot(221)
+ a.plot(self.statistics['error/moments/{0}_err1'.format(ff)][:, 2, 3]**0.5)
+ a.plot(self.statistics['error/moments/{0}_err2'.format(ff)][:, 2, 3]**0.5)
+ a.set_title('absolute')
+ a.set_xlabel('iteration')
+ a = f.add_subplot(222)
+ a.plot(self.statistics['error/moments/{0}_err1_relative'.format(ff)][:, 2, 3]**0.5)
+ a.plot(self.statistics['error/moments/{0}_err2_relative'.format(ff)][:, 2, 3]**0.5)
+ a.set_title('relative')
+ a.set_xlabel('iteration')
+ a = f.add_subplot(223)
+ a.plot(self.statistics['error/histograms/{0}_err1'.format(ff)][-1, :, 3])
+ a.plot(self.statistics['error/histograms/{0}_err2'.format(ff)][-1, :, 3])
+ a.set_title('absolute')
+ a.set_yscale('log')
+ a = f.add_subplot(224)
+ a.plot(self.statistics['error/histograms/{0}_err1_relative'.format(ff)][-1, :, 3])
+ a.plot(self.statistics['error/histograms/{0}_err2_relative'.format(ff)][-1, :, 3])
+ a.set_title('relative')
+ a.set_yscale('log')
+ f.tight_layout()
+ f.suptitle(ff)
+ f.savefig(fname + '_' + ff + '.pdf')
+ plt.close(f)
+ return None
+
+def main():
+ bla = EDNS()
+ bla.launch(
+ ['fluid_error',
+ '--np', '4',
+ '--ntpp', '1',
+ '--tfs', '1',
+ '--field_random_seed', '2']
+ + sys.argv[1:])
+ bla.read_parameters()
+ mm = bla.run_sanity_check()
+ bla.read_error_statistics()
+ bla.plot_error_statistics(err_norm = mm[0].statistics['Uint'])
+ return None
+
+if __name__ == '__main__':
+ main()
diff --git a/examples/field_convergence/README.rst b/examples/field_convergence/README.rst
new file mode 100644
index 0000000000000000000000000000000000000000..530f46f2e2075e496ad0f9de81d334a13426add7
--- /dev/null
+++ b/examples/field_convergence/README.rst
@@ -0,0 +1,13 @@
+Error analysis for Navier-Stokes solvers
+========================================
+
+This is a code to analyze errors of the fluid solver in TurTLE.
+
+There are two main goals to achieve:
+
+* perform error analysis for Navier-Stokes solvers.
+
+* provide an example of TurTLE facilitating "exotic" studies, i.e. use TurTLE as a framework.
+
+This folder contains a basic C++ class to compute the errors, as well as
+a Python wrapper that will launch jobs and read the results.
diff --git a/examples/field_convergence/abstract_fluid_solver.hpp b/examples/field_convergence/abstract_fluid_solver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..37debf4bf1e1b9f74742349671d17131c4df7b56
--- /dev/null
+++ b/examples/field_convergence/abstract_fluid_solver.hpp
@@ -0,0 +1,80 @@
+/******************************************************************************
+* *
+* Copyright 2023 TurTLE team *
+* *
+* This file is part of TurTLE. *
+* *
+* 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, *
+* 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/> *
+* *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de *
+* *
+******************************************************************************/
+
+
+
+#ifndef ABSTRACT_FLUID_SOLVER_HPP
+#define ABSTRACT_FLUID_SOLVER_HPP
+
+#include "field.hpp"
+#include "kspace.hpp"
+
+/** \brief Abstract fluid solver class
+ */
+
+template <typename rnumber>
+class abstract_fluid_solver
+{
+ public:
+ abstract_fluid_solver(
+ const MPI_Comm COMM,
+ const std::string &simname)
+ {}
+ virtual ~abstract_fluid_solver() noexcept(false) {}
+
+ virtual MPI_Comm get_communicator() const = 0;
+ virtual int get_nx() const = 0;
+ virtual int get_ny() const = 0;
+ virtual int get_nz() const = 0;
+
+ virtual int write_checkpoint(void) = 0;
+ virtual int initialize(void) = 0;
+ virtual int step(void) = 0;
+ virtual int do_stats(void) = 0;
+ virtual int finalize(void) = 0;
+
+ virtual int get_cvelocity(field<rnumber, FFTW, THREE> *) = 0;
+ virtual int get_cvorticity(field<rnumber, FFTW, THREE> *) = 0;
+
+ virtual int copy_state_from(abstract_fluid_solver<rnumber> &) = 0;
+
+ virtual kspace<FFTW, SMOOTH>* get_kspace() = 0;
+
+ virtual unsigned get_fftw_plan_rigor() const = 0;
+
+ int get_rvelocity(field<rnumber, FFTW, THREE> *destination)
+ {
+ this->get_cvelocity(destination);
+ destination->ift();
+ return EXIT_SUCCESS;
+ }
+ int get_rvorticity(field<rnumber, FFTW, THREE> *destination)
+ {
+ this->get_cvorticity(destination);
+ destination->ift();
+ return EXIT_SUCCESS;
+ }
+};
+
+#endif//ABSTRACT_FLUID_SOLVER_HPP
+
diff --git a/examples/field_convergence/fluid_error.cpp b/examples/field_convergence/fluid_error.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..245085def49b5e9ad446ddac3cb96a08ccabb5b7
--- /dev/null
+++ b/examples/field_convergence/fluid_error.cpp
@@ -0,0 +1,319 @@
+/******************************************************************************
+* *
+* Copyright 2023 TurTLE team *
+* *
+* This file is part of TurTLE. *
+* *
+* 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, *
+* 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/> *
+* *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de *
+* *
+******************************************************************************/
+
+
+//#include "fluid_error.hpp"
+#include "scope_timer.hpp"
+
+#include <string>
+#include <cmath>
+#include <filesystem>
+
+template <typename rnumber>
+int fluid_error<rnumber>::initialize(void)
+{
+ TIMEZONE("fluid_error::intialize");
+ this->read_parameters();
+
+ this->fs1 = create_fluid_solver<rnumber>(
+ this->tfs,
+ this->get_communicator(),
+ this->simname + "1");
+
+ this->fs2 = create_fluid_solver<rnumber>(
+ this->tfs,
+ this->get_communicator(),
+ this->simname + "2");
+
+ this->fs4 = create_fluid_solver<rnumber>(
+ this->tfs,
+ this->get_communicator(),
+ this->simname + "4");
+
+ this->fs1->initialize();
+ this->fs2->initialize();
+ this->fs4->initialize();
+
+ this->fs2->copy_state_from(*(this->fs1));
+ this->fs4->copy_state_from(*(this->fs1));
+
+ // open stat file
+ if (this->myrank == 0) {
+ // set caching parameters
+ hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
+ herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
+ variable_used_only_in_assert(cache_err);
+ DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
+ this->stat_file = H5Fopen(
+ (this->simname + ".h5").c_str(),
+ H5F_ACC_RDWR,
+ fapl);
+ }
+ int data_file_problem;
+ if (this->myrank == 0) {
+ data_file_problem = this->grow_file_datasets();
+ }
+ MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->get_communicator());
+ if (data_file_problem > 0) {
+ std::cerr <<
+ data_file_problem <<
+ " problems growing file datasets.\ntrying to exit now." <<
+ std::endl;
+ return EXIT_FAILURE;
+ }
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int fluid_error<rnumber>::step(void)
+{
+ TIMEZONE("fluid_error::step");
+ this->fs1->step();
+
+ this->fs2->step();
+ this->fs2->step();
+
+ this->fs4->step();
+ this->fs4->step();
+ this->fs4->step();
+ this->fs4->step();
+
+ this->iteration++;
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int fluid_error<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("fluid_error::write_checkpoint");
+ this->fs1->write_checkpoint();
+ this->fs2->write_checkpoint();
+ this->fs4->write_checkpoint();
+ this->write_iteration();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int fluid_error<rnumber>::finalize(void)
+{
+ TIMEZONE("fluid_error::finalize");
+ this->fs1->finalize();
+ this->fs2->finalize();
+ this->fs4->finalize();
+ delete this->fs1;
+ delete this->fs2;
+ delete this->fs4;
+ if (this->myrank == 0)
+ H5Fclose(this->stat_file);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+inline rnumber relative_difference(
+ const rnumber a,
+ const rnumber b)
+{
+ const rnumber difference = a - b;
+ const rnumber denominator2 = a*a + b*b;
+ return ((denominator2 > std::numeric_limits<rnumber>::epsilon())
+ ? (difference / std::sqrt(denominator2))
+ : 0);
+}
+
+template <typename rnumber>
+int compare_real_fields(
+ kspace<FFTW, SMOOTH>* kk,
+ const field<rnumber, FFTW, THREE>* a,
+ const field<rnumber, FFTW, THREE>* b,
+ const std::string prefix,
+ const std::string suffix,
+ const hid_t stat_group,
+ const int ii,
+ const double max_err_estimate,
+ const double max_err_relative_estimate)
+{
+ TIMEZONE("compare_real_fields");
+ field<rnumber, FFTW, THREE> * err;
+ field<rnumber, FFTW, THREE> * err_relative;
+ err = new field<rnumber, FFTW, THREE>(
+ a->get_nx(), a->get_ny(), a->get_nz(),
+ a->comm,
+ a->fftw_plan_rigor);
+ err_relative = new field<rnumber, FFTW, THREE>(
+ a->get_nx(), a->get_ny(), a->get_nz(),
+ a->comm,
+ a->fftw_plan_rigor);
+ err->RLOOP(
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (auto ii = 0; ii < 3; ii++) {
+ err->rval(rindex, ii) = (
+ a->rval(rindex, ii)
+ - b->rval(rindex, ii));
+ err_relative->rval(rindex, ii) = relative_difference(
+ a->rval(rindex, ii),
+ b->rval(rindex, ii));
+ }
+ });
+ err->compute_stats(
+ kk,
+ stat_group,
+ prefix + "_err" + suffix,
+ ii,
+ max_err_estimate);
+ err_relative->compute_stats(
+ kk,
+ stat_group,
+ prefix + "_err" + suffix + "_relative",
+ ii,
+ max_err_relative_estimate);
+ delete err;
+ delete err_relative;
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int fluid_error<rnumber>::do_stats()
+{
+ TIMEZONE("fluid_error::do_stats");
+ this->fs1->do_stats();
+ this->fs2->do_stats();
+ this->fs4->do_stats();
+ // compute errors and error stats
+ if (!(this->iteration % this->niter_stat == 0))
+ return EXIT_SUCCESS;
+ hid_t stat_group;
+ if (this->myrank == 0)
+ stat_group = H5Gopen(
+ this->stat_file,
+ "statistics",
+ H5P_DEFAULT);
+ else
+ stat_group = 0;
+ field<rnumber, FFTW, THREE> * tmp_vec0;
+ field<rnumber, FFTW, THREE> * tmp_vec1;
+ tmp_vec0 = new field<rnumber, FFTW, THREE>(
+ this->fs1->get_nx(),
+ this->fs1->get_ny(),
+ this->fs1->get_nz(),
+ this->fs1->get_communicator(),
+ this->fs1->get_fftw_plan_rigor());
+ tmp_vec1 = new field<rnumber, FFTW, THREE>(
+ this->fs2->get_nx(),
+ this->fs2->get_ny(),
+ this->fs2->get_nz(),
+ this->fs2->get_communicator(),
+ this->fs2->get_fftw_plan_rigor());
+
+ this->fs1->get_rvelocity(tmp_vec0);
+ this->fs2->get_rvelocity(tmp_vec1);
+ compare_real_fields(
+ this->fs1->get_kspace(),
+ tmp_vec1,
+ tmp_vec0,
+ "velocity",
+ "1",
+ stat_group,
+ this->iteration / niter_stat,
+ this->max_err_estimate,
+ this->max_err_relative_estimate);
+ this->fs1->get_rvorticity(tmp_vec0);
+ this->fs2->get_rvorticity(tmp_vec1);
+ compare_real_fields(
+ this->fs1->get_kspace(),
+ tmp_vec1,
+ tmp_vec0,
+ "vorticity",
+ "1",
+ stat_group,
+ this->iteration / niter_stat,
+ this->max_err_estimate,
+ this->max_err_relative_estimate);
+
+ // clumsy delete/new because in the future we may want different
+ // problem sizes for the different individual DNS.
+ delete tmp_vec0;
+ tmp_vec0 = new field<rnumber, FFTW, THREE>(
+ this->fs4->get_nx(),
+ this->fs4->get_ny(),
+ this->fs4->get_nz(),
+ this->fs4->get_communicator(),
+ this->fs4->get_fftw_plan_rigor());
+
+ this->fs4->get_rvorticity(tmp_vec0);
+ compare_real_fields(
+ this->fs2->get_kspace(),
+ tmp_vec0,
+ tmp_vec1,
+ "vorticity",
+ "2",
+ stat_group,
+ this->iteration / niter_stat,
+ this->max_err_estimate,
+ this->max_err_relative_estimate);
+ this->fs2->get_rvelocity(tmp_vec1);
+ this->fs4->get_rvelocity(tmp_vec0);
+ compare_real_fields(
+ this->fs2->get_kspace(),
+ tmp_vec0,
+ tmp_vec1,
+ "velocity",
+ "2",
+ stat_group,
+ this->iteration / niter_stat,
+ this->max_err_estimate,
+ this->max_err_relative_estimate);
+
+ delete tmp_vec1;
+ delete tmp_vec0;
+
+ if (this->myrank == 0)
+ H5Gclose(stat_group);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int fluid_error<rnumber>::read_parameters(void)
+{
+ TIMEZONE("fluid_error::read_parameters");
+ this->direct_numerical_simulation::read_parameters();
+ hid_t parameter_file = H5Fopen(
+ (this->simname + ".h5").c_str(),
+ H5F_ACC_RDONLY,
+ H5P_DEFAULT);
+ this->tfs = type_of_fluid_solver(hdf5_tools::read_value<int>(
+ parameter_file, "parameters/tfs"));
+ this->max_err_estimate = hdf5_tools::read_value<double>(
+ parameter_file, "parameters/max_err_estimate");
+ this->max_err_relative_estimate = hdf5_tools::read_value<double>(
+ parameter_file, "parameters/max_err_relative_estimate");
+ H5Fclose(parameter_file);
+ MPI_Barrier(this->get_communicator()); // "lock" on parameter file
+ return EXIT_SUCCESS;
+}
+
+template class fluid_error<float>;
+template class fluid_error<double>;
+
diff --git a/examples/field_convergence/fluid_error.hpp b/examples/field_convergence/fluid_error.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..af4fd90924b66bf5e63f1589c3019907247f2407
--- /dev/null
+++ b/examples/field_convergence/fluid_error.hpp
@@ -0,0 +1,67 @@
+/******************************************************************************
+* *
+* Copyright 2023 TurTLE team *
+* *
+* This file is part of TurTLE. *
+* *
+* 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, *
+* 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/> *
+* *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de *
+* *
+******************************************************************************/
+
+
+
+#ifndef FLUID_ERROR_HPP
+#define FLUID_ERROR_HPP
+
+//#include "fluid_solver.hpp"
+
+/** \brief Error analysis for fluid solver.
+ */
+
+template <typename rnumber>
+class fluid_error: public direct_numerical_simulation
+{
+ public:
+
+ /* fundamental ensemble entities */
+ type_of_fluid_solver tfs;
+ abstract_fluid_solver<rnumber> * fs1;
+ abstract_fluid_solver<rnumber> * fs2;
+ abstract_fluid_solver<rnumber> * fs4;
+
+ double max_err_estimate;
+ double max_err_relative_estimate;
+
+ fluid_error(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simname):
+ direct_numerical_simulation(
+ COMMUNICATOR,
+ simname)
+ {}
+ ~fluid_error(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//FLUID_ERROR_HPP
+
diff --git a/examples/field_convergence/fluid_solver.hpp b/examples/field_convergence/fluid_solver.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a963745ef14c50f98ccbeb47ba12daa111028d87
--- /dev/null
+++ b/examples/field_convergence/fluid_solver.hpp
@@ -0,0 +1,182 @@
+/******************************************************************************
+* *
+* Copyright 2023 TurTLE team *
+* *
+* This file is part of TurTLE. *
+* *
+* 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, *
+* 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/> *
+* *
+* Contact: Cristian.Lalescu@mpcdf.mpg.de *
+* *
+******************************************************************************/
+
+
+
+#ifndef FLUID_SOLVER_HPP
+#define FLUID_SOLVER_HPP
+
+//#include "abstract_fluid_solver.hpp"
+#include "full_code/NSVE.hpp"
+#include "full_code/NSE.hpp"
+
+
+template <typename rnumber>
+class fluid_solver_NSVE:
+ public abstract_fluid_solver<rnumber>,
+ public NSVE<rnumber>
+{
+ public:
+ fluid_solver_NSVE(const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name)
+ : abstract_fluid_solver<rnumber>(COMMUNICATOR, simulation_name)
+ , NSVE<rnumber>(COMMUNICATOR, simulation_name){}
+
+ ~fluid_solver_NSVE() noexcept(false) {}
+
+ int get_cvelocity(field<rnumber, FFTW, THREE> *destination)
+ {
+ this->fs->compute_velocity(this->fs->cvorticity);
+ *destination = *(this->fs->cvelocity);
+ return EXIT_SUCCESS;
+ }
+ int get_cvorticity(field<rnumber, FFTW, THREE> *destination)
+ {
+ *destination = *(this->fs->cvorticity);
+ return EXIT_SUCCESS;
+ }
+ int copy_state_from(abstract_fluid_solver<rnumber> &src)
+ {
+ src.get_cvorticity(this->fs->cvorticity);
+ return EXIT_SUCCESS;
+ }
+ kspace<FFTW, SMOOTH>* get_kspace()
+ {
+ return this->fs->kk;
+ }
+ unsigned get_fftw_plan_rigor() const
+ {
+ return this->fs->cvorticity->fftw_plan_rigor;
+ }
+
+ inline MPI_Comm get_communicator() const
+ {
+ return this->NSVE<rnumber>::get_communicator();
+ }
+ inline int get_nx() const {return this->NSVE<rnumber>::get_nx();}
+ inline int get_ny() const {return this->NSVE<rnumber>::get_ny();}
+ inline int get_nz() const {return this->NSVE<rnumber>::get_nz();}
+
+ inline int write_checkpoint(void)
+ {return this->NSVE<rnumber>::write_checkpoint();}
+ inline int initialize(void)
+ {return this->NSVE<rnumber>::initialize();}
+ inline int step(void)
+ {return this->NSVE<rnumber>::step();}
+ inline int do_stats(void)
+ {return this->NSVE<rnumber>::do_stats();}
+ inline int finalize(void)
+ {return this->NSVE<rnumber>::finalize();}
+};
+
+
+template <typename rnumber>
+class fluid_solver_NSE:
+ public abstract_fluid_solver<rnumber>,
+ public NSE<rnumber>
+{
+ public:
+ fluid_solver_NSE(const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name)
+ : abstract_fluid_solver<rnumber>(COMMUNICATOR, simulation_name)
+ , NSE<rnumber>(COMMUNICATOR, simulation_name){}
+ ~fluid_solver_NSE() noexcept(false) {}
+
+ int get_cvelocity(field<rnumber, FFTW, THREE> *destination)
+ {
+ *destination = *(this->velocity);
+ return EXIT_SUCCESS;
+ }
+ int get_cvorticity(field<rnumber, FFTW, THREE> *destination)
+ {
+ compute_curl(this->kk,
+ this->velocity,
+ destination);
+ return EXIT_SUCCESS;
+ }
+ int copy_state_from(abstract_fluid_solver<rnumber> &src)
+ {
+ src.get_cvelocity(this->velocity);
+ return EXIT_SUCCESS;
+ }
+ kspace<FFTW, SMOOTH>* get_kspace()
+ {
+ return this->kk;
+ }
+ unsigned get_fftw_plan_rigor() const
+ {
+ return this->velocity->fftw_plan_rigor;
+ }
+
+ inline MPI_Comm get_communicator() const
+ {
+ return this->NSE<rnumber>::get_communicator();
+ }
+ inline int get_nx() const {return this->NSE<rnumber>::get_nx();}
+ inline int get_ny() const {return this->NSE<rnumber>::get_ny();}
+ inline int get_nz() const {return this->NSE<rnumber>::get_nz();}
+
+ inline int write_checkpoint(void)
+ {return this->NSE<rnumber>::write_checkpoint();}
+ inline int initialize(void)
+ {return this->NSE<rnumber>::initialize();}
+ inline int step(void)
+ {return this->NSE<rnumber>::step();}
+ inline int do_stats(void)
+ {return this->NSE<rnumber>::do_stats();}
+ inline int finalize(void)
+ {return this->NSE<rnumber>::finalize();}
+};
+
+enum type_of_fluid_solver {
+ tfs_NSVE = 0,
+ tfs_NSE = 1
+};
+
+template <typename rnumber>
+abstract_fluid_solver<rnumber> * create_fluid_solver(
+ const type_of_fluid_solver tt,
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name)
+{
+ switch (tt)
+ {
+ case tfs_NSVE:
+ return new fluid_solver_NSVE<rnumber>(
+ COMMUNICATOR,
+ simulation_name);
+ break;
+ case tfs_NSE:
+ return new fluid_solver_NSE<rnumber>(
+ COMMUNICATOR,
+ simulation_name);
+ break;
+ default:
+ DEBUG_MSG("wrong type of fluid solver passed to create_fluid_solver.\n");
+ assert(false);
+ break;
+ }
+}
+
+#endif//FLUID_SOLVER_HPP
+