diff --git a/CMakeLists.txt b/CMakeLists.txt
index 4a05ebf148c3d180f53c744b6c4c4b20e732ac4b..2ff1a6330838916e0b36f5471ff867910934bed5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -230,7 +230,7 @@ if(NOT DEFINED ENV{HDF5_ROOT})
endif()
set(HDF5_PREFER_PARALLEL TRUE)
-set(HDF5_NO_FIND_PACKAGE_CONFIG_FILE TRUE)
+#set(HDF5_NO_FIND_PACKAGE_CONFIG_FILE TRUE)
find_package(HDF5 REQUIRED COMPONENTS C)
target_link_libraries(TurTLE
PRIVATE
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 6c3d3409d7479993f14c1e5e639b7db2fe3801d7..530c0149f72c341f95dd668397bc7201a98c0a52 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -286,7 +286,7 @@ class DNS(_code):
ofile['checkpoint'] = int(0)
if self.dns_type in ['static_field_with_ghost_collisions']:
ofile.create_group('statistics/collisions')
- if (self.dns_type in ['NSE', 'NSVE', 'NSVE_no_output']):
+ if (self.dns_type in ['NSE', 'NSE_alt_dealias', 'NSVE', 'NSVE_no_output']):
return None
return None
def simulation_parser_arguments(
@@ -394,6 +394,10 @@ class DNS(_code):
'NSE',
help = 'plain Navier-Stokes equation')
+ parser_list['NSE_alt_dealias'] = subparsers.add_parser(
+ 'NSE_alt_dealias',
+ help = 'plain Navier-Stokes equation alternate dealiasing')
+
parser_list['NSVE'] = subparsers.add_parser(
'NSVE',
help = 'plain Navier-Stokes vorticity equation')
@@ -833,12 +837,13 @@ class DNS(_code):
need_field = False
if self.check_current_vorticity_exists:
need_field = True
- if self.dns_type == 'NSE':
+ if self.dns_type in ['NSE', 'NSE_alt_dealias']:
checkpoint_field = 'velocity'
else:
checkpoint_field = 'vorticity'
if self.dns_type in [
'NSE',
+ 'NSE_alt_dealias',
'NSVE',
'NSVE_no_output',
'static_field',
@@ -960,7 +965,7 @@ class DNS(_code):
if (('test' in self.dns_type) or
(self.dns_type in ['kraichnan_field', 'NSVE'])):
self.check_current_vorticity_exists = False
- if self.dns_type == 'NSE':
+ if self.dns_type in ['NSE', 'NSE_alt_dealias']:
self.check_current_vorticity_exists = False
self.check_current_velocity_exists = True
self.run(
diff --git a/TurTLE/PP.py b/TurTLE/PP.py
index 2056cf3296cd8b76f00a903050e75f76b0af3684..e68f321483fa53698fc3fd4754355e4912ba519f 100644
--- a/TurTLE/PP.py
+++ b/TurTLE/PP.py
@@ -76,16 +76,27 @@ class PP(_code):
def generate_default_parameters(self):
# these parameters are relevant for all PP classes
self.parameters['fftw_plan_rigor'] = 'FFTW_ESTIMATE'
+ self.parameter_description['fftw_plan_rigor'] = 'FFTW plan rigor to use. One of `FFTW_ESTIMATE`, `FFTW_MEASURE`, `FFTW_PATIENT`. Please see FFTW documentation.'
self.parameters['dealias_type'] = int(1)
+ self.parameter_description['dealias_type'] = 'Dealiasing mehtod to use, integer. Options are: two-thirds (0) or smooth (1).'
self.parameters['dkx'] = float(1.0)
+ self.parameter_description['dkx'] = 'Smallest wavenumber in the x direction for a pseudo-spectral run.'
self.parameters['dky'] = float(1.0)
+ self.parameter_description['dky'] = 'Smallest wavenumber in the y direction for a pseudo-spectral run.'
self.parameters['dkz'] = float(1.0)
+ self.parameter_description['dkz'] = 'Smallest wavenumber in the z direction for a pseudo-spectral run.'
self.parameters['nu'] = float(0.1)
+ self.parameter_description['nu'] = 'Viscosity value used in the equations, given in code units.'
self.parameters['fmode'] = int(1)
+ self.parameter_description['fmode'] = 'Forcing parameter: mode to use for the Kolmogorov forcing.'
self.parameters['famplitude'] = float(0.5)
+ self.parameter_description['famplitude'] = 'Forcing parameter: amplitude of Kolmogorov forcing, in code units.'
self.parameters['fk0'] = float(2.0)
+ self.parameter_description['fk0'] = 'Forcing parameter: if forcing acts on wavenumber band, this is the smallest wavenumber where it acts (in code units).'
self.parameters['fk1'] = float(4.0)
+ self.parameter_description['fk1'] = 'Forcing parameter: if forcing acts on wavenumber band, this is the largest wavenumber where it acts (in code units).'
self.parameters['forcing_type'] = 'linear'
+ self.parameter_description['forcing_type'] = 'Forcing parameter: what type of force to use.'
self.pp_parameters = {}
self.pp_parameters['iteration_list'] = np.zeros(1).astype(int)
return None
@@ -396,7 +407,7 @@ class PP(_code):
', parameters = self.extra_postprocessing_parameters("' +
pp_type +
'"))')
- return None
+ return subparsers
def prepare_launch(
self,
args = []):
@@ -756,6 +767,8 @@ class PP(_code):
with h5py.File(self.get_fields_fname(), 'a') as ff:
ff.require_group('vorticity')
ff.require_group('vorticity/complex')
+ ff.require_group('velocity')
+ ff.require_group('velocity/complex')
checkpoint_file_list = [self.get_checkpoint_fname(checkpoint = cp)
for cp in range(df['checkpoint'][()]+1)]
for cpf_name in checkpoint_file_list:
@@ -763,13 +776,21 @@ class PP(_code):
cpf = h5py.File(cpf_name, 'r')
if 'vorticity' not in cpf.keys():
print('file ', cpf_name, ' does not have vorticity group')
- continue
else:
for iter_name in cpf['vorticity/complex'].keys():
if iter_name not in ff['vorticity/complex'].keys():
ff['vorticity/complex/' + iter_name] = h5py.ExternalLink(
cpf_name,
'vorticity/complex/' + iter_name)
+ if 'velocity' not in cpf.keys():
+ print('file ', cpf_name, ' does not have velocity group')
+ else:
+ if 'complex' in cpf['velocity'].keys():
+ for iter_name in cpf['velocity/complex'].keys():
+ if iter_name not in ff['velocity/complex'].keys():
+ ff['velocity/complex/' + iter_name] = h5py.ExternalLink(
+ cpf_name,
+ 'velocity/complex/' + iter_name)
cpf.close()
if self.dns_type == 'resize':
with h5py.File(os.path.join(self.work_dir, self.pp_parameters['new_simname'] + '.h5'), 'a') as ff:
diff --git a/TurTLE/_code.py b/TurTLE/_code.py
index 0490ad8e6e88788ee05f15d0784e5b57ac1bccfd..06ea006f9c68a5486f752e41ce50cc9ad8a297ff 100644
--- a/TurTLE/_code.py
+++ b/TurTLE/_code.py
@@ -96,7 +96,7 @@ class _code(_base):
outfile.write('find_package(MPI REQUIRED COMPONENTS C)\n')
outfile.write('set(HDF5_STATIC ON)\n')
outfile.write('set(HDF5_PREFER_PARALLEL TRUE)\n')
- outfile.write('set(HDF5_NO_FIND_PACKAGE_CONFIG_FILE TRUE)\n')
+ #outfile.write('set(HDF5_NO_FIND_PACKAGE_CONFIG_FILE TRUE)\n')
outfile.write('find_package(HDF5 REQUIRED COMPONENTS C)\n')
if not no_debug:
outfile.write('message(VERBOSE "found HDF5 include dir at ${HDF5_C_INCLUDE_DIR}")\n')
diff --git a/TurTLE/test/test_fftw.py b/TurTLE/test/test_fftw.py
index a1034eccdc9a36a830de4b73e0718b81086decff..4f8cd4fb2fa71f0a31a2fcf4dde3f52040871d34 100644
--- a/TurTLE/test/test_fftw.py
+++ b/TurTLE/test/test_fftw.py
@@ -39,7 +39,6 @@ def main():
'--ntpp', '1',
'--wd', './'] +
sys.argv[1:])
- df = h5py.File(c.simname + '.h5', 'r')
df = h5py.File(c.simname + '_fields.h5', 'r')
field0_complex = df['field0/complex/0'][...]
field1_complex = df['field1/complex/0'][...]
diff --git a/TurTLE/test/test_turtle_NSE_dealias.py b/TurTLE/test/test_turtle_NSE_dealias.py
new file mode 100644
index 0000000000000000000000000000000000000000..0cb00592e571ca3d6800dbd23ddfef483dd3a432
--- /dev/null
+++ b/TurTLE/test/test_turtle_NSE_dealias.py
@@ -0,0 +1,254 @@
+#! /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 TurTLE
+from TurTLE import DNS, PP
+
+from TurTLE_addons import NSReader
+
+try:
+ import matplotlib.pyplot as plt
+except:
+ plt = None
+
+
+def main_basic(
+ wd = './'
+ ):
+ niterations = 256
+ njobs = 2
+ c0 = DNS()
+ c0.launch(
+ ['NSVE',
+ '-n', '32',
+ '--src-simname', 'B32p1e4',
+ '--src-wd', TurTLE.data_dir,
+ '--src-iteration', '0',
+ '--simname', 'nsve_for_comparison',
+ '--np', '4',
+ '--ntpp', '1',
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(niterations),
+ '--forcing_type', 'linear',
+ '--niter_stat', '1',
+ '--checkpoints_per_file', '{0}'.format(3),
+ '--njobs', '{0}'.format(njobs),
+ '--wd', wd] +
+ sys.argv[1:])
+ cc = PP()
+ cc.launch(
+ ['get_velocity',
+ '--simname', 'nsve_for_comparison',
+ '--np', '4',
+ '--ntpp', '1',
+ '--iter0', '0',
+ '--iter1', '{0}'.format(niterations),
+ '--wd', wd] +
+ sys.argv[1:])
+ c1 = DNS()
+ c1.launch(
+ ['NSE',
+ '-n', '32',
+ '--src-simname', 'nsve_for_comparison',
+ '--src-wd', wd,
+ '--src-iteration', '{0}'.format(niterations),
+ '--simname', 'nse1',
+ '--np', '4',
+ '--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', wd] +
+ sys.argv[1:])
+ c2 = DNS()
+ c2.launch(
+ ['NSE_alt_dealias',
+ '-n', '32',
+ '--src-simname', 'nsve_for_comparison',
+ '--src-wd', wd,
+ '--src-iteration', '{0}'.format(niterations),
+ '--simname', 'nse2',
+ '--np', '4',
+ '--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', wd] +
+ sys.argv[1:])
+ f1 = h5py.File(
+ os.path.join(wd, 'nse1.h5'), 'r')
+ f2 = h5py.File(
+ os.path.join(wd, 'nse2.h5'), 'r')
+ if not type(plt) == type(None):
+ e1 = 0.5*f1['statistics/moments/velocity'][:, 2, 3]
+ e2 = 0.5*f2['statistics/moments/velocity'][:, 2, 3]
+ fig = plt.figure()
+ a = fig.add_subplot(211)
+ a.set_title('energy')
+ a.plot(e1)
+ a.plot(e2)
+ e1 = 0.5*f1['statistics/moments/vorticity'][:, 2, 3]
+ e2 = 0.5*f2['statistics/moments/vorticity'][:, 2, 3]
+ a = fig.add_subplot(212)
+ a.set_title('enstrophy')
+ a.plot(e1)
+ a.plot(e2)
+ fig.tight_layout()
+ fig.savefig('comparison.pdf')
+ plt.close(fig)
+ f1.close()
+ f2.close()
+ print('SUCCESS! Basic test passed.')
+ return None
+
+
+def main_long():
+ niterations = 2048
+ N = 64
+ njobs = 4
+ c0 = DNS()
+ c0.launch(
+ ['NSVE',
+ '-n', '{0}'.format(N),
+ '--src-simname', 'B32p1e4',
+ '--src-wd', TurTLE.data_dir,
+ '--src-iteration', '0',
+ '--simname', 'nsve_for_comparison',
+ '--np', '4',
+ '--ntpp', '1',
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(niterations),
+ '--forcing_type', 'linear',
+ '--niter_stat', '1',
+ '--checkpoints_per_file', '{0}'.format(3),
+ '--wd', './'] +
+ sys.argv[1:])
+ cc = PP()
+ cc.launch(
+ ['get_velocity',
+ '--simname', 'nsve_for_comparison',
+ '--np', '4',
+ '--ntpp', '1',
+ '--iter0', '0',
+ '--iter1', '{0}'.format(niterations)] +
+ sys.argv[1:])
+ c1 = DNS()
+ c1.launch(
+ ['NSE',
+ '-n', '{0}'.format(N),
+ '--src-simname', 'nsve_for_comparison',
+ '--src-wd', './',
+ '--src-iteration', '{0}'.format(niterations),
+ '--simname', 'nse1',
+ '--np', '4',
+ '--ntpp', '1',
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(niterations//64),
+ '--niter_stat', '1',
+ '--forcing_type', 'linear',
+ '--njobs', '{0}'.format(njobs),
+ '--wd', './'] +
+ sys.argv[1:])
+ c2 = DNS()
+ c2.launch(
+ ['NSE_alt_dealias',
+ '-n', '{0}'.format(N),
+ '--src-simname', 'nsve_for_comparison',
+ '--src-wd', './',
+ '--src-iteration', '{0}'.format(niterations),
+ '--simname', 'nse2',
+ '--np', '4',
+ '--ntpp', '1',
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(niterations//64),
+ '--niter_stat', '1',
+ '--forcing_type', 'linear',
+ '--njobs', '{0}'.format(njobs),
+ '--wd', './'] +
+ sys.argv[1:])
+ cc1 = NSReader(simname = 'nse1')
+ cc2 = NSReader(simname = 'nse2')
+ cc1.read_full_tk_spectrum()
+ cc2.read_full_tk_spectrum()
+ if not type(plt) == type(None):
+ fig = plt.figure()
+ a = fig.add_subplot(211)
+ a.set_ylabel('energy')
+ a.set_xlabel('iteration')
+ a.plot(cc1.statistics['energy(t)'])
+ a.plot(cc2.statistics['energy(t)'])
+ a = fig.add_subplot(212)
+ a.set_ylabel('enstrophy')
+ a.set_xlabel('iteration')
+ a.plot(cc1.statistics['enstrophy(t)'])
+ a.plot(cc2.statistics['enstrophy(t)'])
+ fig.suptitle('total physical times are {0:.1f} and {1:.1f} Tint'.format(
+ cc1.statistics['t'][-1] / cc1.statistics['Tint'],
+ cc2.statistics['t'][-1] / cc2.statistics['Tint']
+ ))
+ fig.tight_layout()
+ fig.savefig('comparison_averages.pdf')
+ plt.close(fig)
+ fig = plt.figure()
+ a = fig.add_subplot(111)
+ a.plot(cc1.statistics['kshell'], cc1.statistics['energy(t, k)'][0], color = 'C0', label = 'dealias after')
+ a.plot(cc2.statistics['kshell'], cc2.statistics['energy(t, k)'][0], color = 'C1', label = 'dealias before')
+ a.plot(cc1.statistics['kshell'], cc1.statistics['energy(t, k)'][2048], color = 'C0', dashes = (1, 1))
+ a.plot(cc2.statistics['kshell'], cc2.statistics['energy(t, k)'][2048], color = 'C1', dashes = (1, 1))
+ a.plot(cc1.statistics['kshell'], cc1.statistics['energy(t, k)'][4096], color = 'C0', dashes = (2, 3))
+ a.plot(cc2.statistics['kshell'], cc2.statistics['energy(t, k)'][4096], color = 'C1', dashes = (2, 3))
+ a.plot(cc1.statistics['kshell'], cc1.statistics['energy(t, k)'][8192], color = 'C0', dashes = (3, 5))
+ a.plot(cc2.statistics['kshell'], cc2.statistics['energy(t, k)'][8192], color = 'C1', dashes = (3, 5))
+ a.set_xscale('log')
+ a.set_yscale('log')
+ a.legend(loc = 'best')
+ fig.tight_layout()
+ fig.savefig('comparison_spectra.pdf')
+ plt.close(fig)
+ return None
+
+if __name__ == '__main__':
+ #main_basic()
+ main_long()
+
diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt
index 6fbbc57cf33f1df54112d4181af51d2b6c6c5f87..25d39ce899f89bc5bb5eee575914a01d2e82f510 100644
--- a/cpp/CMakeLists.txt
+++ b/cpp/CMakeLists.txt
@@ -28,6 +28,7 @@ set(cpp_for_lib
${CMAKE_CURRENT_LIST_DIR}/full_code/code_base.cpp
${CMAKE_CURRENT_LIST_DIR}/full_code/direct_numerical_simulation.cpp
${CMAKE_CURRENT_LIST_DIR}/full_code/NSE.cpp
+ ${CMAKE_CURRENT_LIST_DIR}/full_code/NSE_alt_dealias.cpp
${CMAKE_CURRENT_LIST_DIR}/full_code/NSVE.cpp
${CMAKE_CURRENT_LIST_DIR}/full_code/static_field.cpp
${CMAKE_CURRENT_LIST_DIR}/full_code/kraichnan_field.cpp
@@ -117,6 +118,7 @@ set(hpp_for_lib
${CMAKE_CURRENT_LIST_DIR}/full_code/main_code.hpp
${CMAKE_CURRENT_LIST_DIR}/full_code/native_binary_to_hdf5.hpp
${CMAKE_CURRENT_LIST_DIR}/full_code/NSE.hpp
+ ${CMAKE_CURRENT_LIST_DIR}/full_code/NSE_alt_dealias.hpp
${CMAKE_CURRENT_LIST_DIR}/full_code/NSVEcomplex_particles.hpp
${CMAKE_CURRENT_LIST_DIR}/full_code/NSVE_field_stats.hpp
${CMAKE_CURRENT_LIST_DIR}/full_code/NSVE.hpp
diff --git a/cpp/field.cpp b/cpp/field.cpp
index b2644571c659fe4ce12c28a400d0e896194070ab..fe03fa374867304e80c5acdd3b6fe4910abada31 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -273,7 +273,9 @@ int field<rnumber, be, fc>::io(
/* file dataset has same dimensions as field */
TIMEZONE("field::io");
hid_t file_id, dset_id, plist_id;
+ file_id = H5I_BADID;
dset_id = H5I_BADID;
+ plist_id = H5I_BADID;
std::string representation = std::string(
this->real_space_representation ?
"real" : "complex");
@@ -918,8 +920,6 @@ int field<rnumber, be, fc>::write_filtered(
}
}
}
- //DEBUG_MSG("count[0] = %ld, offset[0] = %ld\n",
- // count[0], offset[0]);
if (nz>=2)
{
assert(nz%2==0);
@@ -931,8 +931,6 @@ int field<rnumber, be, fc>::write_filtered(
offset[1] = cz*nz/2;
memshape [1] = this->clayout->sizes[1];
memoffset[1] = cz*(this->clayout->sizes[1] - nz/2);
- //DEBUG_MSG("cz = %d, count[1] + offset[1] = %ld\n",
- // cz, count[1] + offset[1]);
//now write data
mspace = H5Screate_simple(ndim(fc), memshape, NULL);
@@ -949,7 +947,6 @@ int field<rnumber, be, fc>::write_filtered(
offset[1] = 0;
memshape [1] = this->clayout->sizes[1];
memoffset[1] = 0;
- //DEBUG_MSG("filtered write nz=1\n");
//now write data
mspace = H5Screate_simple(ndim(fc), memshape, NULL);
@@ -1099,6 +1096,10 @@ void field<rnumber, be, fc>::compute_rspace_stats(
hsize_t dims[ndim(fc)-1];
int ndims;
dset = H5Dopen(group, ("moments/" + dset_name).c_str(), H5P_DEFAULT);
+ if(dset < 0)
+ {
+ DEBUG_MSG("couldn't open dataset 'moments/%s'\n", dset_name.c_str());
+ }
wspace = H5Dget_space(dset);
ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
// the following -1 comes from -3+2
@@ -1205,7 +1206,7 @@ void field<rnumber, be, fc>::compute_rspace_stats(
local_moments[0*nvals+3] = val_tmp[3];
if (val_tmp[3] > local_moments[9*nvals+3])
local_moments[9*nvals+3] = val_tmp[3];
- const int bin = int(floor(val_tmp[3]*2/binsize[3]));
+ const int bin = int(std::floor(val_tmp[3]*2/binsize[3]));
if (bin >= 0 && bin < nbins){
local_hist[bin*nvals+3]++;
}
@@ -1216,7 +1217,7 @@ void field<rnumber, be, fc>::compute_rspace_stats(
local_moments[0*nvals+i] = val_tmp[i];
if (val_tmp[i] > local_moments[(nmoments-1)*nvals+i])
local_moments[(nmoments-1)*nvals+i] = val_tmp[i];
- const int bin = int(floor((val_tmp[i] + max_estimate[i]) / binsize[i]));
+ const int bin = int(std::floor((val_tmp[i] + max_estimate[i]) / binsize[i]));
if (bin >= 0 && bin < nbins)
local_hist[bin*nvals+i]++;
}
@@ -1460,21 +1461,6 @@ void field<rnumber, be, fc>::compute_rspace_zaverage(
delete[] zaverage;
}
-template <typename rnumber,
- field_backend be,
- field_components fc>
-void field<rnumber, be, fc>::normalize()
-{
- this->RLOOP(
- [&](const ptrdiff_t rindex,
- const ptrdiff_t xindex,
- const ptrdiff_t yindex,
- const ptrdiff_t zindex){
- for (unsigned int i=0; i<ncomp(fc); i++)
- this->data[rindex*ncomp(fc) + i] /= this->npoints;
- });
-}
-
template <typename rnumber,
field_backend be,
field_components fc>
@@ -2274,15 +2260,12 @@ void field<rnumber, be, fc>::compute_stats(
max_estimate_vector);
did_rspace = true;
this->dft();
- // normalize
- TIMEZONE("field::normalize");
- for (hsize_t tmp_index=0; tmp_index<this->rmemlayout->local_size; tmp_index++)
- this->data[tmp_index] /= this->npoints;
+ this->normalize();
}
// what follows gave me a headache until I found this link:
// http://stackoverflow.com/questions/8256636/expected-primary-expression-error-on-template-method-using
kk->template cospectrum<rnumber, fc>(
- (typename fftw_interface<rnumber>::complex*)this->data,
+ (typename fftw_interface<rnumber>::complex*)(this->data),
group,
dset_name + "_" + dset_name,
toffset);
@@ -2511,6 +2494,169 @@ int invert_curl(
return EXIT_SUCCESS;
}
+inline int bin_id(
+ const double value,
+ const double max_estimate,
+ const double binsize,
+ const int total_number_of_bins)
+{
+ assert(max_estimate > 0);
+ assert(binsize > 0);
+ assert(total_number_of_bins > 0);
+ return std::min(
+ std::max(int(0), int(std::floor((value + max_estimate) / binsize))),
+ total_number_of_bins-1);
+}
+
+template <typename rnumber,
+ field_backend be,
+ field_components fc>
+int conditional_rspace_PDF(
+ field<rnumber, be, fc>* ff,
+ field<rnumber, be, ONE>* cc,
+ const hid_t group,
+ const std::string dset_name,
+ const hsize_t toffset,
+ const std::vector<double> max_ff_estimate,
+ const double max_cc_estimate)
+{
+ TIMEZONE("conditional_rspace_PDF");
+ assert(ff->real_space_representation);
+ assert(cc->real_space_representation);
+ assert(max_ff_estimate.size() >= ncomp(fc));
+ int nbins_ff, nbins_cc;
+
+ if (cc->myrank == 0)
+ {
+ hid_t dset, wspace;
+ hsize_t dims[5];
+ int ndims;
+ variable_used_only_in_assert(ndims);
+ dset = H5Dopen(
+ group,
+ ("histograms/" + dset_name).c_str(),
+ H5P_DEFAULT);
+ if (dset <= 0) {
+ throw std::runtime_error("Couldn't open " + dset_name + " for reading dimensions.\n");
+ }
+ wspace = H5Dget_space(dset);
+ ndims = H5Sget_simple_extent_dims(wspace, dims, NULL);
+ if (fc == ONE)
+ assert(ndims == 3);
+ else if (fc == THREE) {
+ assert(ndims == 4);
+ assert(dims[3] == 3);
+ } else if (fc == THREExTHREE) {
+ assert(ndims == 5);
+ assert(dims[3] == 3);
+ assert(dims[4] == 3);
+ }
+ H5Sclose(wspace);
+ H5Dclose(dset);
+ nbins_cc = dims[1];
+ nbins_ff = dims[2];
+ }
+ {
+ TIMEZONE("MPI_Bcast");
+ MPI_Bcast(&nbins_cc, 1, MPI_INT, 0, cc->comm);
+ MPI_Bcast(&nbins_ff, 1, MPI_INT, 0, cc->comm);
+ }
+
+ // allocate histogram
+ shared_array<ptrdiff_t> local_hist_threaded(
+ nbins_cc*nbins_ff*ncomp(fc),
+ [&](ptrdiff_t* local_hist){
+ std::fill_n(local_hist, nbins_cc*nbins_ff*ncomp(fc), 0);
+ });
+
+ /// set up bin sizes
+ const double bin_cc_size = 2*max_cc_estimate / nbins_cc;
+ std::vector<double> bin_ff_size;
+ bin_ff_size.resize(ncomp(fc));
+ for (unsigned int i=0; i<ncomp(fc); i++)
+ {
+ bin_ff_size[i] = 2*max_ff_estimate[i] / nbins_ff;
+ }
+
+ {
+ TIMEZONE("field::RLOOP");
+ cc->RLOOP(
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ ptrdiff_t *local_hist = local_hist_threaded.getMine();
+ const int bin_cc = bin_id(
+ cc->rval(rindex),
+ max_cc_estimate,
+ bin_cc_size,
+ nbins_cc);
+ for (unsigned int component = 0;
+ component < ncomp(fc);
+ component++) {
+ const int bin_ff = bin_id(
+ ff->rval(rindex, component),
+ max_ff_estimate[component],
+ bin_ff_size[component],
+ nbins_ff);
+ local_hist[
+ (nbins_ff * bin_cc + bin_ff)*ncomp(fc)+component]++;
+ }
+ });
+ }
+ // OpenMP merge
+ local_hist_threaded.mergeParallel();
+ ptrdiff_t *hist = new ptrdiff_t[nbins_cc*nbins_ff*ncomp(fc)];
+ // MPI merge
+ MPI_Allreduce(
+ (void*)local_hist_threaded.getMasterData(),
+ (void*)hist,
+ nbins_cc*nbins_ff*ncomp(fc),
+ MPI_INT64_T, MPI_SUM, cc->comm);
+
+ // output
+ if (cc->myrank == 0)
+ {
+ TIMEZONE("root-work");
+ hid_t dset, wspace, mspace;
+ hsize_t count[5], offset[5];
+ offset[0] = toffset;
+ offset[1] = 0;
+ offset[2] = 0;
+ offset[3] = 0;
+ offset[4] = 0;
+ count[0] = 1;
+ count[1] = nbins_cc;
+ count[2] = nbins_ff;
+ count[3] = 3;
+ count[4] = 3;
+ dset = H5Dopen(
+ group,
+ ("histograms/" + dset_name).c_str(),
+ H5P_DEFAULT);
+ if (dset <= 0) {
+ throw std::runtime_error(
+ "Couldn't open " + dset_name + " for writing.\n");
+ }
+ wspace = H5Dget_space(dset);
+ if (fc == ONE) {
+ mspace = H5Screate_simple(3, count, NULL);
+ } else if (fc == THREE) {
+ mspace = H5Screate_simple(4, count, NULL);
+ } else if (fc == THREExTHREE) {
+ mspace = H5Screate_simple(5, count, NULL);
+ }
+ H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
+ H5Dwrite(dset, H5T_NATIVE_INT64, mspace, wspace, H5P_DEFAULT, hist);
+ H5Sclose(wspace);
+ H5Sclose(mspace);
+ H5Dclose(dset);
+ }
+
+ delete[] hist;
+ return EXIT_SUCCESS;
+}
+
template <typename rnumber,
field_backend be,
field_components fc>
@@ -2638,25 +2784,25 @@ int joint_rspace_PDF(
{
double val1 = f1->rval(rindex, i);
mag1 += val1*val1;
- int bin1 = int(floor((val1 + max_f1_estimate[i])/bin1size[i]));
+ int bin1 = int(std::floor((val1 + max_f1_estimate[i])/bin1size[i]));
mag2 = 0.0;
for (unsigned int j=0; j<3; j++)
{
double val2 = f2->rval(rindex, j);
mag2 += val2*val2;
- int bin2 = int(floor((val2 + max_f2_estimate[j])/bin2size[j]));
+ int bin2 = int(std::floor((val2 + max_f2_estimate[j])/bin2size[j]));
if ((bin1 >= 0 && bin1 < nbins) &&
(bin2 >= 0 && bin2 < nbins))
local_histc[(bin1*nbins + bin2)*9 + i*3 + j]++;
}
}
- bin1 = int(floor(std::sqrt(mag1)/bin1size[3]));
- bin2 = int(floor(std::sqrt(mag2)/bin2size[3]));
+ bin1 = int(std::floor(std::sqrt(mag1)/bin1size[3]));
+ bin2 = int(std::floor(std::sqrt(mag2)/bin2size[3]));
}
else if (fc == ONE)
{
- bin1 = int(floor((f1->rval(rindex) + max_f1_estimate[0])/bin1size[3]));
- bin2 = int(floor((f2->rval(rindex) + max_f2_estimate[0])/bin2size[3]));
+ bin1 = int(std::floor((f1->rval(rindex) + max_f1_estimate[0])/bin1size[3]));
+ bin2 = int(std::floor((f2->rval(rindex) + max_f2_estimate[0])/bin2size[3]));
}
if ((bin1 >= 0 && bin1 < nbins) &&
(bin2 >= 0 && bin2 < nbins))
@@ -2690,19 +2836,19 @@ int joint_rspace_PDF(
TIMEZONE("root-work");
hid_t dset, wspace, mspace;
hsize_t count[5], offset[5];
+ offset[0] = toffset;
+ offset[1] = 0;
+ offset[2] = 0;
+ count[0] = 1;
+ count[1] = nbins;
+ count[2] = nbins;
if (fc == THREE)
{
dset = H5Dopen(group, dsetc.c_str(), H5P_DEFAULT);
assert(dset > 0);
wspace = H5Dget_space(dset);
- offset[0] = toffset;
- offset[1] = 0;
- offset[2] = 0;
offset[3] = 0;
offset[4] = 0;
- count[0] = 1;
- count[1] = nbins;
- count[2] = nbins;
count[3] = 3;
count[4] = 3;
mspace = H5Screate_simple(5, count, NULL);
@@ -2714,12 +2860,6 @@ int joint_rspace_PDF(
}
dset = H5Dopen(group, dsetm.c_str(), H5P_DEFAULT);
assert(dset > 0);
- offset[0] = toffset;
- offset[1] = 0;
- offset[2] = 0;
- count[0] = 1;
- count[1] = nbins;
- count[2] = nbins;
mspace = H5Screate_simple(3, count, NULL);
wspace = H5Dget_space(dset);
H5Sselect_hyperslab(wspace, H5S_SELECT_SET, offset, NULL, count, NULL);
@@ -2817,9 +2957,9 @@ int joint_rspace_3PDF(
int bin2 = 0;
int bin3 = 0;
- bin1 = int(floor((f1->rval(rindex) + max_f1_estimate[0])/bin1size[0]));
- bin2 = int(floor((f2->rval(rindex) + max_f2_estimate[0])/bin2size[0]));
- bin3 = int(floor((f3->rval(rindex) + max_f3_estimate[0])/bin3size[0]));
+ bin1 = int(std::floor((f1->rval(rindex) + max_f1_estimate[0])/bin1size[0]));
+ bin2 = int(std::floor((f2->rval(rindex) + max_f2_estimate[0])/bin2size[0]));
+ bin3 = int(std::floor((f3->rval(rindex) + max_f3_estimate[0])/bin3size[0]));
if ((bin1 >= 0 && bin1 < nbins) &&
(bin2 >= 0 && bin2 < nbins) &&
(bin3 >= 0 && bin3 < nbins))
@@ -3263,15 +3403,6 @@ int generate_random_phase_field(
output_field->cval(cindex, 0) = 0;
output_field->cval(cindex, 1) = 0;
}
- //if ((xindex == output_field->clayout->sizes[2]-1) &&
- // (zindex == output_field->clayout->sizes[1]-1))
- //{
- // DEBUG_MSG("yindex = %ld, phase_field = %g + j%g\n",
- // yindex,
- // output_field->cval(cindex, 0),
- // output_field->cval(cindex, 1)
- // );
- //}
});
}
@@ -3496,6 +3627,57 @@ template int joint_rspace_PDF<double, FFTW, ONE>(
const std::vector<double>,
const std::vector<double>);
+template int conditional_rspace_PDF<float, FFTW, THREExTHREE>(
+ field<float, FFTW, THREExTHREE> *,
+ field<float, FFTW, ONE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const double);
+template int conditional_rspace_PDF<double, FFTW, THREExTHREE>(
+ field<double, FFTW, THREExTHREE> *,
+ field<double, FFTW, ONE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const double);
+
+template int conditional_rspace_PDF<float, FFTW, THREE>(
+ field<float, FFTW, THREE> *,
+ field<float, FFTW, ONE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const double);
+template int conditional_rspace_PDF<double, FFTW, THREE>(
+ field<double, FFTW, THREE> *,
+ field<double, FFTW, ONE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const double);
+
+template int conditional_rspace_PDF<float, FFTW, ONE>(
+ field<float, FFTW, ONE> *,
+ field<float, FFTW, ONE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const double);
+template int conditional_rspace_PDF<double, FFTW, ONE>(
+ field<double, FFTW, ONE> *,
+ field<double, FFTW, ONE> *,
+ const hid_t,
+ const std::string,
+ const hsize_t,
+ const std::vector<double>,
+ const double);
+
template int joint_rspace_3PDF<float, FFTW>(
field<float, FFTW, ONE> *,
field<float, FFTW, ONE> *,
diff --git a/cpp/field.hpp b/cpp/field.hpp
index b98eb74fe2f05bc09aa6ba975f20ac8445073e05..6759af2683592c92a15e77d773a5d998efdd27d8 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -115,7 +115,6 @@ class field
/* essential FFT stuff */
void dft();
void ift();
- void normalize();
void symmetrize();
void symmetrize_FFT();
void symmetrize_alternate();
@@ -306,6 +305,36 @@ class field
}
finish_mpi_profiling_zone(turtle_mpi_pcontrol::FIELD);
}
+ void normalize()
+ {
+ // duplicates RLOOP for loops, BUT x loop is over rmemlayout, not rlayout.
+ start_mpi_profiling_zone(turtle_mpi_pcontrol::FIELD);
+ switch(be)
+ {
+ case FFTW:
+ #pragma omp parallel
+ {
+ const hsize_t start = OmpUtils::ForIntervalStart(this->rlayout->subsizes[1]);
+ const hsize_t end = OmpUtils::ForIntervalEnd(this->rlayout->subsizes[1]);
+
+ for (hsize_t zindex = 0; zindex < this->rlayout->subsizes[0]; zindex++)
+ //#pragma omp simd
+ for (hsize_t yindex = start; yindex < end; yindex++)
+ {
+ const ptrdiff_t rindex = (
+ zindex * this->rlayout->subsizes[1] + yindex)*(
+ this->rmemlayout->subsizes[2])*ncomp(fc);
+ #pragma omp simd
+ for (hsize_t xindex = 0;
+ xindex < this->rmemlayout->subsizes[2]*ncomp(fc);
+ xindex++)
+ this->data[rindex+xindex] /= this->npoints;
+ }
+ }
+ break;
+ }
+ finish_mpi_profiling_zone(turtle_mpi_pcontrol::FIELD);
+ }
ptrdiff_t get_cindex(
ptrdiff_t xindex,
ptrdiff_t yindex,
@@ -382,6 +411,18 @@ int compute_curl(
field<rnumber, be, THREE> *source,
field<rnumber, be, THREE> *destination);
+template <typename rnumber,
+ field_backend be,
+ field_components fc>
+int conditional_rspace_PDF(
+ field<rnumber, be, fc>* ff,
+ field<rnumber, be, ONE>* cc,
+ const hid_t group,
+ const std::string dset_name,
+ const hsize_t toffset,
+ const std::vector<double> max_ff_estimate,
+ const double max_cc_estimate);
+
template <typename rnumber,
field_backend be,
field_components fc>
diff --git a/cpp/full_code/NSE_alt_dealias.cpp b/cpp/full_code/NSE_alt_dealias.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..44a325fb9e6306e179fe15903771321a2f80bb01
--- /dev/null
+++ b/cpp/full_code/NSE_alt_dealias.cpp
@@ -0,0 +1,604 @@
+/******************************************************************************
+* *
+* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
+* *
+* This file is part of bfps. *
+* *
+* bfps 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. *
+* *
+* bfps 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@ds.mpg.de *
+* *
+******************************************************************************/
+
+
+
+#include <string>
+#include <cmath>
+#include "NSE_alt_dealias.hpp"
+#include "scope_timer.hpp"
+#include "fftw_tools.hpp"
+#include "hdf5_tools.hpp"
+#include "shared_array.hpp"
+
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::initialize(void)
+{
+ TIMEZONE("NSE_alt_dealias::initialize");
+ this->read_iteration();
+ this->read_parameters();
+ 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->comm);
+ if (data_file_problem > 0)
+ {
+ std::cerr <<
+ data_file_problem <<
+ " problems growing file datasets.\ntrying to exit now." <<
+ std::endl;
+ return EXIT_FAILURE;
+ }
+
+ /* allocate vector fields */
+ this->velocity = new field<rnumber, FFTW, THREE>(
+ 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->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->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->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->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);
+
+ {
+ this->velocity->real_space_representation = false;
+ this->velocity->io(
+ this->get_current_fname(),
+ "velocity",
+ this->iteration,
+ true);
+ TIMEZONE("NSE_alt_dealias::initialize::read_initial_condition");
+ this->kk->template low_pass<rnumber, THREE>(
+ this->velocity->get_cdata(),
+ this->kk->kM);
+ this->kk->template force_divfree<rnumber>(
+ this->velocity->get_cdata());
+ this->velocity->symmetrize();
+ MPI_Barrier(this->comm);
+ }
+
+ if (this->myrank == 0 && this->iteration == 0)
+ this->kk->store(stat_file);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::step(void)
+{
+ TIMEZONE("NSE_alt_dealias::step");
+ this->iteration++;
+ //return this->Euler_step();
+ return this->ShuOsher();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::add_field_band(
+ field<rnumber, FFTW, THREE> *dst,
+ field<rnumber, FFTW, THREE> *src,
+ const double k0,
+ const double k1,
+ const double prefactor)
+{
+ TIMEZONE("NSE_alt_dealias::add_field_band");
+ this->kk->CLOOP(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ const double knorm = sqrt(this->kk->kx[xindex]*this->kk->kx[xindex] +
+ this->kk->ky[yindex]*this->kk->ky[yindex] +
+ this->kk->kz[zindex]*this->kk->kz[zindex]);
+ if ((k0 <= knorm) &&
+ (k1 >= knorm))
+ for (int c=0; c<3; c++)
+ for (int i=0; i<2; i++)
+ dst->cval(cindex,c,i) += prefactor*src->cval(cindex,c,i);
+ }
+ );
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::add_forcing_term(
+ field<rnumber, FFTW, THREE> *vel,
+ field<rnumber, FFTW, THREE> *dst)
+{
+ TIMEZONE("NSE_alt_dealias::add_forcing_term");
+ if (this->forcing_type == std::string("linear"))
+ {
+ return this->add_field_band(
+ dst, this->velocity,
+ this->fk0, this->fk1,
+ this->famplitude);
+ }
+ if ((this->forcing_type == std::string("fixed_energy_injection_rate")) ||
+ (this->forcing_type == std::string("fixed_energy_injection_rate_and_drag")))
+ {
+ // first, compute energy in shell
+ shared_array<double> local_energy_in_shell(1);
+ double energy_in_shell = 0;
+ this->kk->CLOOP_K2_NXMODES(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2,
+ const int nxmodes){
+ const double knorm = sqrt(k2);
+ if ((k2 > 0) &&
+ (this->fk0 <= knorm) &&
+ (this->fk1 >= knorm))
+ *local_energy_in_shell.getMine() += nxmodes*(
+ vel->cval(cindex, 0, 0)*vel->cval(cindex, 0, 0) + vel->cval(cindex, 0, 1)*vel->cval(cindex, 0, 1) +
+ vel->cval(cindex, 1, 0)*vel->cval(cindex, 1, 0) + vel->cval(cindex, 1, 1)*vel->cval(cindex, 1, 1) +
+ vel->cval(cindex, 2, 0)*vel->cval(cindex, 2, 0) + vel->cval(cindex, 2, 1)*vel->cval(cindex, 2, 1)
+ );
+ }
+ );
+ local_energy_in_shell.mergeParallel();
+ MPI_Allreduce(
+ local_energy_in_shell.getMasterData(),
+ &energy_in_shell,
+ 1,
+ MPI_DOUBLE,
+ MPI_SUM,
+ vel->comm);
+ // we should divide by 2, if we wanted energy;
+ // but then we would need to multiply the amplitude by 2 anyway,
+ // because what we really care about is force dotted into velocity,
+ // without the division by 2.
+
+ // now, modify amplitudes
+ if (energy_in_shell < 10*std::numeric_limits<rnumber>::epsilon())
+ energy_in_shell = 1;
+ double temp_famplitude = this->injection_rate / energy_in_shell;
+ this->add_field_band(
+ dst, vel,
+ this->fk0, this->fk1,
+ temp_famplitude);
+ // and add drag if desired
+ if (this->forcing_type == std::string("fixed_energy_injection_rate_and_drag"))
+ this->add_field_band(
+ dst, vel,
+ this->fmode, this->fmode + (this->fk1 - this->fk0),
+ -this->friction_coefficient);
+ return EXIT_SUCCESS;
+ }
+ if (this->forcing_type == std::string("fixed_energy"))
+ return EXIT_SUCCESS;
+ else
+ {
+ DEBUG_MSG("unknown forcing type printed on next line\n%s", this->forcing_type.c_str());
+ throw std::invalid_argument("unknown forcing type");
+ }
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::impose_forcing(
+ field<rnumber, FFTW, THREE> *dst)
+{
+ TIMEZONE("NSE_alt_dealias::impose_forcing");
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::add_nonlinear_term(
+ const field<rnumber, FFTW, THREE> *vel,
+ field<rnumber, FFTW, THREE> *dst)
+{
+ TIMEZONE("NSE_alt_dealias::add_nonlinear_term");
+ dst->real_space_representation = false;
+
+ /* This computes the nonlinear term coming from vel, for the Navier Stokes
+ * equations. Using vel=\f$u\f$ and dst=\f$a^{nl}\f$, we compute
+ * \f[
+ * a^{nl}_i = - u_j \partial_j u_i
+ * \f]
+ * The result is given in Fourier representation.
+ * */
+ const int sign_array[2] = {1, -1};
+
+ /* copy velocity */
+ this->tmp_vec_field0->real_space_representation = false;
+ *this->tmp_vec_field0 = *vel;
+
+ //dealias
+ this->kk->template dealias<rnumber, THREE>(this->tmp_vec_field0->get_cdata());
+ // put velocity in real space
+ this->tmp_vec_field0->ift();
+
+ /* compute uu */
+ /* store 11 22 33 components */
+ this->tmp_vec_field1->real_space_representation = true;
+ this->tmp_vec_field1->RLOOP(
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (int cc=0; cc<3; cc++)
+ this->tmp_vec_field1->rval(rindex,cc) = (
+ 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 */
+ this->tmp_vec_field1->dft();
+
+ this->kk->CLOOP(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ dst->cval(cindex, 0, ccc) += (
+ sign_array[ccc]*this->kk->kx[xindex] * this->tmp_vec_field1->cval(cindex, 0, (ccc+1)%2));
+ dst->cval(cindex, 1, ccc) += (
+ sign_array[ccc]*this->kk->ky[yindex] * this->tmp_vec_field1->cval(cindex, 1, (ccc+1)%2));
+ dst->cval(cindex, 2, ccc) += (
+ sign_array[ccc]*this->kk->kz[zindex] * this->tmp_vec_field1->cval(cindex, 2, (ccc+1)%2));
+ }
+ });
+
+ /* store 12 23 31 components */
+ this->tmp_vec_field1->real_space_representation = true;
+ this->tmp_vec_field1->RLOOP(
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (int cc=0; cc<3; cc++)
+ this->tmp_vec_field1->rval(rindex,cc) = (
+ 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 */
+ this->tmp_vec_field1->dft();
+
+ this->kk->CLOOP(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ dst->cval(cindex, 0, ccc) += (
+ sign_array[ccc]*this->kk->ky[yindex] * this->tmp_vec_field1->cval(cindex, 0, (ccc+1)%2) +
+ sign_array[ccc]*this->kk->kz[zindex] * this->tmp_vec_field1->cval(cindex, 2, (ccc+1)%2));
+ dst->cval(cindex, 1, ccc) += (
+ sign_array[ccc]*this->kk->kz[zindex] * this->tmp_vec_field1->cval(cindex, 1, (ccc+1)%2) +
+ sign_array[ccc]*this->kk->kx[xindex] * this->tmp_vec_field1->cval(cindex, 0, (ccc+1)%2));
+ dst->cval(cindex, 2, ccc) += (
+ sign_array[ccc]*this->kk->kx[xindex] * this->tmp_vec_field1->cval(cindex, 2, (ccc+1)%2) +
+ sign_array[ccc]*this->kk->ky[yindex] * this->tmp_vec_field1->cval(cindex, 1, (ccc+1)%2));
+ }
+ });
+
+ dst->symmetrize();
+ /* done */
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::Euler_step(void)
+{
+ // technically this is not a pure Euler method,
+ // because I use an exponential factor for the diffusion term
+ TIMEZONE("NSE_alt_dealias::Euler_step");
+
+ // 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);
+
+ // add nonlinear term
+ this->add_nonlinear_term(
+ this->velocity,
+ acc);
+
+ // add forcing
+ this->add_forcing_term(this->velocity, acc);
+
+ // perform Euler step
+ 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) += this->dt*acc->cval(cindex, cc, ccc);
+ }
+ });
+
+ // enforce incompressibility
+ this->kk->template force_divfree<rnumber>(this->velocity->get_cdata());
+
+ // apply diffusion exponential
+ 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);
+ for (int cc = 0; cc < 3; cc++)
+ for (int ccc = 0; ccc < 2; ccc++)
+ {
+ this->velocity->cval(cindex, cc, ccc) *= factor;
+ }
+ });
+
+ // enforce Hermitian symmetry
+ this->velocity->symmetrize();
+
+ // delete temporary field
+ delete acc;
+ return EXIT_SUCCESS;
+}
+
+/** \brief Time step with Shu Osher method
+ *
+ * The Navier Stokes equations are integrated with a
+ * third-order Runge-Kutta method [shu1988jcp]_, which is an explicit
+ * Runge-Kutta method with the Butcher tableau
+ *
+ * | | | | |
+ * |:--|:--|:--|:--|
+ * | 0 | | | |
+ * | 1 | 1 | | |
+ * | 1/2 | 1/4 | 1/4 | |
+ * | | 1/6 | 1/6 | 2/3 |
+ *
+ * In addition to the stability properties described in [shu1988jcp]_, this method has the advantage that it is memory-efficient, requiring only two additional field allocations, as can be seen from
+ * \f{eqnarray*}{
+ * \hat{\boldsymbol w}_1(\boldsymbol k) &= {\hat{\boldsymbol u}}(\boldsymbol k, t)e^{-\nu k^2 h} + h \textrm{NL}[\hat{\boldsymbol u}(\boldsymbol k, t)]e^{-\nu k^2 h}, \\
+ * \hat{\boldsymbol w}_2(\boldsymbol k) &= \frac{3}{4}\hat{\boldsymbol u}( \boldsymbol k, t)e^{-\nu k^2 h/2}
+ * +
+ * \frac{1}{4}(\hat{\boldsymbol w}_1(\boldsymbol k) + h \textrm{NL}[\hat{\boldsymbol w}_1(\boldsymbol k)])e^{\nu k^2 h/2}, \\
+ * \hat{\boldsymbol u}(\boldsymbol k, t+h) &= \frac{1}{3}\hat{\boldsymbol u}(\boldsymbol k, t)e^{-\nu k^2 h} +
+ * \frac{2}{3}(\hat{\boldsymbol w}_2(\boldsymbol k) + h \textrm{NL}[\hat{\boldsymbol w}_2(\boldsymbol k)])e^{-\nu k^2 h/2},
+ * \f}
+ */
+template <typename rnumber>
+int NSE_alt_dealias<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);
+ 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);
+ 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_velocity2 = rnumber(0);
+ this->add_nonlinear_term(this->tmp_velocity1, this->tmp_velocity2);
+ this->add_forcing_term(this->tmp_velocity1, this->tmp_velocity2);
+ 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;
+ 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->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;
+ 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->symmetrize();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("NSE_alt_dealias::write_checkpoint");
+ this->update_checkpoint();
+
+ assert(!this->velocity->real_space_representation);
+ std::string fname = this->get_current_fname();
+ this->velocity->io(
+ fname,
+ "velocity",
+ this->iteration,
+ false);
+
+ this->write_iteration();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::finalize(void)
+{
+ TIMEZONE("NSE_alt_dealias::finalize");
+ if (this->myrank == 0)
+ H5Fclose(this->stat_file);
+ delete this->kk;
+ delete this->tmp_velocity1;
+ delete this->tmp_velocity2;
+ delete this->tmp_vec_field1;
+ delete this->tmp_vec_field0;
+ delete this->velocity;
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute standard statistics for velocity and vorticity fields.
+ *
+ * IMPORTANT: at the end of this subroutine, `this->fs->cvelocity` contains
+ * the Fourier space representation of the velocity field, and
+ * `this->tmp_vec_field` contains the real space representation of the
+ * velocity field.
+ * This behavior is relied upon in the `NSEparticles` class, so please
+ * don't break it.
+ */
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::do_stats()
+{
+ TIMEZONE("NSE_alt_dealias::do_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;
+
+ this->tmp_vec_field0->real_space_representation = false;
+ *this->tmp_vec_field0 = *this->velocity;
+
+ this->tmp_vec_field0->compute_stats(
+ this->kk,
+ stat_group,
+ "velocity",
+ this->iteration / niter_stat,
+ this->max_velocity_estimate/sqrt(3));
+
+ compute_curl(this->kk,
+ this->velocity,
+ this->tmp_vec_field0);
+
+ this->tmp_vec_field0->compute_stats(
+ this->kk,
+ stat_group,
+ "vorticity",
+ this->iteration / niter_stat,
+ this->max_vorticity_estimate/sqrt(3));
+
+ if (this->myrank == 0)
+ H5Gclose(stat_group);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSE_alt_dealias<rnumber>::read_parameters(void)
+{
+ TIMEZONE("NSE_alt_dealias::read_parameters");
+ this->direct_numerical_simulation::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ 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->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");
+ this->fk0 = hdf5_tools::read_value<double>(parameter_file, "parameters/fk0");
+ this->fk1 = hdf5_tools::read_value<double>(parameter_file, "parameters/fk1");
+ this->energy = hdf5_tools::read_value<double>(parameter_file, "parameters/energy");
+ this->histogram_bins = hdf5_tools::read_value<int>(parameter_file, "parameters/histogram_bins");
+ this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_velocity_estimate");
+ this->max_vorticity_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_vorticity_estimate");
+ this->forcing_type = hdf5_tools::read_string(parameter_file, "parameters/forcing_type");
+ this->fftw_plan_rigor = hdf5_tools::read_string(parameter_file, "parameters/fftw_plan_rigor");
+ H5Fclose(parameter_file);
+ MPI_Barrier(this->comm);
+ return EXIT_SUCCESS;
+}
+
+template class NSE_alt_dealias<float>;
+template class NSE_alt_dealias<double>;
+
diff --git a/cpp/full_code/NSE_alt_dealias.hpp b/cpp/full_code/NSE_alt_dealias.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..eb81f48f614772e7cdd4c70679ea6aaefdd23987
--- /dev/null
+++ b/cpp/full_code/NSE_alt_dealias.hpp
@@ -0,0 +1,111 @@
+/**********************************************************************
+* *
+* Copyright 2017 Max Planck Institute *
+* for Dynamics and Self-Organization *
+* *
+* This file is part of bfps. *
+* *
+* bfps 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. *
+* *
+* bfps 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 bfps. If not, see <http://www.gnu.org/licenses/> *
+* *
+* Contact: Cristian.Lalescu@ds.mpg.de *
+* *
+**********************************************************************/
+
+
+
+#ifndef NSE_ALT_DEALIAS_HPP
+#define NSE_ALT_DEALIAS_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/direct_numerical_simulation.hpp"
+
+/** \brief Navier-Stokes solver.
+ *
+ * Solves the Navier Stokes equation.
+ *
+ * Allocates fields and a kspace object.
+ */
+
+
+template <typename rnumber>
+class NSE_alt_dealias: public direct_numerical_simulation
+{
+ public:
+ std::string name;
+
+ /* parameters that are read in read_parameters */
+ double dt;
+ double famplitude;
+ double friction_coefficient;
+ double fk0;
+ double fk1;
+ double energy;
+ double injection_rate;
+ int fmode;
+ std::string forcing_type;
+ int histogram_bins;
+ double max_velocity_estimate;
+ double max_vorticity_estimate;
+ double nu;
+ std::string fftw_plan_rigor;
+
+ /* other stuff */
+ kspace<FFTW, SMOOTH> *kk;
+ field<rnumber, FFTW, THREE> *velocity;
+ field<rnumber, FFTW, THREE> *tmp_velocity1; // for timestep
+ field<rnumber, FFTW, THREE> *tmp_velocity2; // for timestep
+ field<rnumber, FFTW, THREE> *tmp_vec_field0; // can be changed during methods
+ field<rnumber, FFTW, THREE> *tmp_vec_field1; // can be changed during methods
+
+ NSE_alt_dealias(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ direct_numerical_simulation(
+ COMMUNICATOR,
+ simulation_name){}
+ ~NSE_alt_dealias(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int Euler_step(void);
+ int ShuOsher(void);
+
+ int add_field_band(
+ field<rnumber, FFTW, THREE> *dst,
+ field<rnumber, FFTW, THREE> *src,
+ const double k0,
+ const double k1,
+ const double prefactor);
+ virtual int add_nonlinear_term(
+ const field<rnumber, FFTW, THREE> *vel,
+ field<rnumber, FFTW, THREE> *dst);
+ int add_forcing_term(
+ field<rnumber, FFTW, THREE> *vel,
+ field<rnumber, FFTW, THREE> *dst);
+ int impose_forcing(
+ field<rnumber, FFTW, THREE> *vel);
+
+ virtual int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//NSE_ALT_DEALIAS_HPP
+
diff --git a/cpp/full_code/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index b2af1679fd4cb8d3a8c9cfb53651d218dcdc9a6c..9fbb3a4165a5b95cb964bca0b46010f58b820432 100644
--- a/cpp/full_code/NSVE_field_stats.cpp
+++ b/cpp/full_code/NSVE_field_stats.cpp
@@ -28,11 +28,18 @@
#include "scope_timer.hpp"
#include "hdf5_tools.hpp"
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::read_parameters(void)
+{
+ this->postprocess::read_parameters();
+ return EXIT_SUCCESS;
+}
+
template <typename rnumber>
int NSVE_field_stats<rnumber>::initialize(void)
{
TIMEZONE("NSVE_field_stats::initialize");
- this->postprocess::read_parameters();
+ this->read_parameters();
this->vorticity = new field<rnumber, FFTW, THREE>(
nx, ny, nz,
this->comm,
diff --git a/cpp/full_code/NSVE_field_stats.hpp b/cpp/full_code/NSVE_field_stats.hpp
index 5e65e80661f8710392c1a38e33bccf1739015c0e..34292ef3bf44e6179c321833e4e48442c79bd17c 100644
--- a/cpp/full_code/NSVE_field_stats.hpp
+++ b/cpp/full_code/NSVE_field_stats.hpp
@@ -47,6 +47,7 @@ class NSVE_field_stats: public postprocess
simulation_name){}
virtual ~NSVE_field_stats() noexcept(false){}
+ virtual int read_parameters(void);
virtual int initialize(void);
virtual int work_on_current_iteration(void);
virtual int finalize(void);
diff --git a/cpp/hdf5_tools.cpp b/cpp/hdf5_tools.cpp
index 8aa91ec49c9159e725a34415717f83a3b112cf0e..42f626fcf3c0891263dc6fe3b0e2c33d0d599550 100644
--- a/cpp/hdf5_tools.cpp
+++ b/cpp/hdf5_tools.cpp
@@ -245,17 +245,26 @@ std::vector<number> hdf5_tools::read_vector(
const std::string dset_name)
{
TIMEZONE("hdf5_tools::read_vector");
+ DEBUG_MSG("hdf5_tools::read_vector dset_name = %s\n", dset_name.c_str());
std::vector<number> result;
hsize_t vector_length;
// first, read size of array
hid_t dset, dspace;
hid_t mem_dtype = H5Tcopy(hdf5_tools::hdf5_type_id<number>());
dset = H5Dopen(group, dset_name.c_str(), H5P_DEFAULT);
+ if (dset < 0)
+ {
+ DEBUG_MSG("had an error trying to open %s.\n H5Dopen returned dset = %ld\n", dset_name.c_str(), dset);
+ }
dspace = H5Dget_space(dset);
assert(H5Sget_simple_extent_ndims(dspace) == 1);
H5Sget_simple_extent_dims(dspace, &vector_length, NULL);
result.resize(vector_length);
- H5Dread(dset, mem_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &result.front());
+ herr_t status = H5Dread(dset, mem_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &result.front());
+ if (status < 0)
+ {
+ DEBUG_MSG("had an error trying to read %s.\n H5Dread returned status %ld\n", dset_name.c_str(), status);
+ }
H5Sclose(dspace);
H5Dclose(dset);
H5Tclose(mem_dtype);
diff --git a/cpp/kspace.cpp b/cpp/kspace.cpp
index fd5b0f9b6c1f654fc9c6ce814d27dd6979b9dd6e..d199d721a6fcd145327064a750a420693879e923 100644
--- a/cpp/kspace.cpp
+++ b/cpp/kspace.cpp
@@ -348,6 +348,7 @@ void kspace<be, dt>::Gauss_filter(
typename fftw_interface<rnumber>::complex *__restrict__ a,
const double sigma)
{
+ DEBUG_MSG("kspace::Gauss_filter called with sigma = %g\n", sigma);
this->template Gauss_n_filter<rnumber, fc, 2>(a, sigma);
}
@@ -403,7 +404,7 @@ void kspace<be, dt>::Gauss_n_filter(
typename fftw_interface<rnumber>::complex *__restrict__ a,
const double sigma)
{
- const double prefactor = -pow(sigma, n)/2;
+ const double prefactor = -std::pow(sigma, n)/2;
this->CLOOP_K2(
[&](const ptrdiff_t cindex,
const ptrdiff_t xindex,
@@ -411,9 +412,12 @@ void kspace<be, dt>::Gauss_n_filter(
const ptrdiff_t zindex,
const double k2){
{
- for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
- ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= std::exp(
+ const double factor = std::exp(
prefactor*hack_second_power<n>(k2));
+ for (unsigned int tcounter=0;
+ tcounter<2*ncomp(fc);
+ tcounter++)
+ ((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= factor;
}
});
}
@@ -473,43 +477,31 @@ int kspace<be, dt>::filter(
const double wavenumber,
std::string filter_type)
{
- if (filter_type == std::string("sharp_Fourier_sphere"))
- {
+ if (filter_type == std::string("sharp_Fourier_sphere")) {
this->template low_pass<rnumber, fc>(
a,
wavenumber);
- }
- else if (filter_type == std::string("Gauss"))
- {
+ } else if (filter_type == std::string("Gauss")) {
this->template Gauss_filter<rnumber, fc>(
a,
2*acos(0.)/wavenumber);
- }
- else if (filter_type == std::string("Gauss4"))
- {
+ } else if (filter_type == std::string("Gauss4")) {
this->template Gauss_n_filter<rnumber, fc, 4>(
a,
2*acos(0.)/wavenumber);
- }
- else if (filter_type == std::string("Gauss6"))
- {
+ } else if (filter_type == std::string("Gauss6")) {
this->template Gauss_n_filter<rnumber, fc, 6>(
a,
2*acos(0.)/wavenumber);
- }
- else if (filter_type == std::string("Gauss8"))
- {
+ } else if (filter_type == std::string("Gauss8")) {
this->template Gauss_n_filter<rnumber, fc, 8>(
a,
2*acos(0.)/wavenumber);
- }
- else if (filter_type == std::string("ball"))
- {
+ } else if (filter_type == std::string("ball")) {
this->template ball_filter<rnumber, fc>(
a,
2*acos(0.)/wavenumber);
- }
- else{
+ } else {
throw std::runtime_error(
"Filter type not available.");
}
@@ -869,6 +861,10 @@ void kspace<be, dt>::cospectrum(
hid_t dset, wspace, mspace;
hsize_t count[(ndim(fc)-2)*2], offset[(ndim(fc)-2)*2], dims[(ndim(fc)-2)*2];
dset = H5Dopen(group, ("spectra/" + dset_name).c_str(), H5P_DEFAULT);
+ if (dset < 0)
+ {
+ DEBUG_MSG("Error opening dataset %s\n", dset_name.c_str());
+ }
wspace = H5Dget_space(dset);
H5Sget_simple_extent_dims(wspace, dims, NULL);
switch (fc)