diff --git a/CMakeLists.txt b/CMakeLists.txt
index 414c0d6ab760e593c4480105d01a10cb13e36222..9bebb964949d1e275e05e2936f6caebb1aa1db2b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -3,20 +3,20 @@
# Copyright 2019 Max Planck Institute #
# for Dynamics and Self-Organization #
# #
-# This file is part of bfps. #
+# This file is part of TurTLE. #
# #
-# bfps is free software: you can redistribute it and/or modify #
+# TurTLE is free software: you can redistribute it and/or modify #
# it under the terms of the GNU General Public License as published #
# by the Free Software Foundation, either version 3 of the License, #
# or (at your option) any later version. #
# #
-# bfps is distributed in the hope that it will be useful, #
+# TurTLE is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
-# along with bfps. If not, see <http://www.gnu.org/licenses/> #
+# along with TurTLE. If not, see <http://www.gnu.org/licenses/> #
# #
# Contact: Cristian.Lalescu@ds.mpg.de #
# #
@@ -272,6 +272,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/write_rpressure.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/bandpass_stats.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_single_to_double.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/resize.cpp
@@ -300,6 +301,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_Stokes_particles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.cpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -320,6 +322,7 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_output_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/get_rfields.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/write_rpressure.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/bandpass_stats.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_single_to_double.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/resize.hpp
@@ -348,14 +351,16 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_Stokes_particles.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEp_extra_sampling.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer_2nd_order.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer_empty.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particles_input.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particles_output.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particles_system.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particles_system_with_p2p.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/alltoall_exchanger.hpp
- ${PROJECT_SOURCE_DIR}/cpp/particles/env_utils.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/lock_free_bool_array.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/p2p_computer_empty.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/p2p_computer.hpp
@@ -368,7 +373,6 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/particles/particles_distr_mpi.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_field_computer.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_generic_interp.hpp
- ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer_empty.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_input_hdf5.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_output_hdf5.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_output_mpiio.hpp
@@ -382,6 +386,7 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE_no_output.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles_no_output.hpp
${PROJECT_SOURCE_DIR}/cpp/base.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/env_utils.hpp
${PROJECT_SOURCE_DIR}/cpp/fftw_interface.hpp
${PROJECT_SOURCE_DIR}/cpp/bfps_timer.hpp
${PROJECT_SOURCE_DIR}/cpp/omputils.hpp
@@ -408,6 +413,13 @@ configure_file(${PROJECT_SOURCE_DIR}/cmake/TurTLEConfig.cmake.in ${PROJECT_BINAR
install(FILES "${PROJECT_BINARY_DIR}/TurTLEConfig.cmake" DESTINATION lib/)
export(TARGETS TurTLE FILE "${PROJECT_BINARY_DIR}/TurTLELibraryDepends.cmake")
install(EXPORT TURTLE_EXPORT DESTINATION lib/)
+if(EXISTS "${PROJECT_BINARY_DIR}/bash_setup_for_TurTLE.sh")
+ install(
+ FILES "${PROJECT_BINARY_DIR}/bash_setup_for_TurTLE.sh"
+ PERMISSIONS OWNER_READ GROUP_READ WORLD_READ
+ DESTINATION "lib/"
+ )
+endif()
#####################################################################################
@@ -420,3 +432,70 @@ else()
install(CODE "execute_process(COMMAND ${CMAKE_COMMAND} -E copy ${PROJECT_SOURCE_DIR}/pc_host_info.py ${PROJECT_BINARY_DIR}/python/TurTLE/host_info.py)")
endif()
install(CODE "execute_process(COMMAND python3 ${PROJECT_SOURCE_DIR}/setup.py install --force --prefix=${CMAKE_INSTALL_PREFIX} WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/python/)")
+
+
+#####################################################################################
+## Add tests
+include(CTest)
+set(TEST_OUTPUT_DIRECTORY "${PROJECT_BINARY_DIR}/test_runs/")
+enable_testing()
+if (BUILD_TESTING)
+ file(MAKE_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ ### basic functionality
+ add_test(
+ NAME test_fftw
+ COMMAND turtle.test_fftw
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_Parseval
+ COMMAND turtle.test_Parseval
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ ### compare DNS output to stored results
+ add_test(
+ NAME test_NSVEparticles
+ COMMAND turtle.test_NSVEparticles --ntpp 2
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ ### simple runs of post-processing tools
+ add_test(
+ NAME test_pp_single_to_double
+ COMMAND turtle PP field_single_to_double --simname dns_nsveparticles --iter0 32 --iter1 32
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_get_rfields
+ COMMAND turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_write_rpressure
+ COMMAND turtle PP write_rpressure --simname dns_nsveparticles --iter0 0 --iter1 64
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_joint_acc_vel_stats
+ COMMAND turtle PP joint_acc_vel_stats --simname dns_nsveparticles --iter0 0 --iter1 64
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_pp_resize
+ COMMAND turtle PP resize --simname dns_nsveparticles --new_nx 96 --new_ny 96 --new_nz 96 --new_simname dns_nsveparticles_resized
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ ### simple runs of different DNS
+ add_test(
+ NAME test_NSVEp_extra_sampling
+ COMMAND turtle DNS NSVEp_extra_sampling -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsvep_extra_sampling --nparticles 1000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_NSVEcomplex_particles
+ COMMAND turtle DNS NSVEcomplex_particles -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsvecomplex_particles --nparticles 1000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_static_field
+ COMMAND turtle DNS static_field --simname dns_static_field --nparticles 10000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_kraichnan_field
+ COMMAND turtle DNS kraichnan_field --simname dns_kraichnan_field --dtfactor 0.05 --nparticles 10000 --ntpp 2
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+ add_test(
+ NAME test_NSVE_Stokes_particles
+ COMMAND turtle DNS NSVE_Stokes_particles -n 32 --src-simname dns_nsveparticles --src-iteration 32 --simname dns_nsve_tokes_particles --nparticles 10000
+ WORKING_DIRECTORY ${TEST_OUTPUT_DIRECTORY})
+endif(BUILD_TESTING)
+
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index d77ffef50133dab6348408f5f0d75637b3c359e2..c0ac8328553e9c8a5e6785fef24a2150de288cb4 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -439,7 +439,14 @@ class DNS(_code):
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)
- if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles', 'static_field', 'static_field_with_ghost_collisions', 'kraichnan_field']:
+ if self.dns_type in [
+ 'NSVEparticles_no_output',
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
+ 'NSVEparticles',
+ 'static_field',
+ 'static_field_with_ghost_collisions',
+ 'kraichnan_field']:
assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
assert (self.parameters['niter_out'] % self.parameters['niter_part'] == 0)
_code.write_par(self, iter0 = iter0)
@@ -654,40 +661,45 @@ class DNS(_code):
'NSVEparticles',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
+ parser_NSVE_Stokes_particles = subparsers.add_parser(
+ 'NSVE_Stokes_particles',
+ help = 'plain Navier-Stokes vorticity formulation, with passive Stokes drag particles')
+
parser_NSVEp2p = subparsers.add_parser(
'NSVEcomplex_particles',
help = 'plain Navier-Stokes vorticity formulation, with oriented active particles')
parser_NSVEp_extra = subparsers.add_parser(
'NSVEp_extra_sampling',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers, that sample velocity gradient, as well as pressure and its derivatives.')
-
- for parser in [
+ for pp in [
'NSVE',
'NSVE_no_output',
'NSVEparticles_no_output',
'NSVEp2',
'NSVEp2p',
+ 'NSVE_Stokes_particles',
'NSVEp_extra',
'static_field',
'static_field_with_ghost_collisions',
'kraichnan_field']:
- eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
- eval('self.job_parser_arguments({0})'.format('parser_' + parser))
- eval('self.parameters_to_parser_arguments({0})'.format('parser_' + parser))
+ eval('self.simulation_parser_arguments({0})'.format('parser_' + pp))
+ eval('self.job_parser_arguments({0})'.format('parser_' + pp))
+ eval('self.parameters_to_parser_arguments({0})'.format('parser_' + pp))
eval('self.parameters_to_parser_arguments('
'parser_{0},'
- 'self.generate_extra_parameters(\'{0}\'))'.format(parser))
- for parser in [
+ 'self.generate_extra_parameters(\'{0}\'))'.format(pp))
+ for pp in [
'NSVEparticles_no_output',
'NSVEp2',
'NSVEp2p',
+ 'NSVE_Stokes_particles',
'NSVEp_extra',
'static_field',
'kraichnan_field']:
- eval('self.particle_parser_arguments({0})'.format('parser_' + parser))
+ eval('self.particle_parser_arguments({0})'.format('parser_' + pp))
eval('self.parameters_to_parser_arguments('
'parser_{0},'
- 'self.NSVEp_extra_parameters)'.format(parser))
+ 'self.NSVEp_extra_parameters)'.format(pp))
return None
def generate_extra_parameters(
self,
@@ -698,7 +710,11 @@ class DNS(_code):
pars['field_random_seed'] = int(1)
pars['spectrum_slope'] = float(-5./3)
pars['spectrum_k_cutoff'] = float(16)
- pars['spectrum_coefficient'] = float(1)
+ pars['spectrum_coefficient'] = float(0.1)
+ if dns_type == 'NSVE_Stokes_particles':
+ pars['initial_field_amplitude'] = float(0.0)
+ pars['initial_particle_vel'] = float(0.05)
+ pars['drag_coefficient'] = float(0.1)
return pars
def prepare_launch(
self,
@@ -737,6 +753,7 @@ class DNS(_code):
if self.dns_type in [
'NSVEparticles',
'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
'NSVEparticles_no_output',
'NSVEp_extra_sampling',
'static_field',
@@ -812,7 +829,15 @@ class DNS(_code):
# hardcoded FFTW complex representation size
field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
checkpoint_size = field_size
- if self.dns_type in ['static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling', 'static_field_with_ghost_collisions', 'kraichnan_field']:
+ if self.dns_type in [
+ 'kraichnan_field',
+ 'static_field',
+ 'static_field_with_ghost_collisions',
+ 'NSVEparticles',
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
+ 'NSVEparticles_no_output',
+ 'NSVEp_extra_sampling']:
rhs_size = self.parameters['tracers0_integration_steps']
if type(opt.tracers0_integration_steps) != type(None):
rhs_size = opt.tracers0_integration_steps
@@ -852,7 +877,9 @@ class DNS(_code):
integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps']
if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys():
integration_steps = self.parameters['tracers{0}_integration_steps'.format(species)]
- if self.dns_type == 'NSVEcomplex_particles' and species == 0:
+ if self.dns_type in [
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles'] and species == 0:
ncomponents = 6
with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
nn = self.parameters['nparticles']
@@ -884,12 +911,18 @@ class DNS(_code):
if nn > batch_size:
dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
if dset.shape[1] == 6:
- dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
+ if self.dns_type == 'NSVE_Stokes_particles':
+ dset[cc*batch_size:(cc+1)*batch_size, 3:] = self.parameters['initial_particle_vel']*get_random_versors(batch_size)
+ else:
+ dset[cc*batch_size:(cc+1)*batch_size, 3:] = get_random_versors(batch_size)
nn -= batch_size
else:
dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
if dset.shape[1] == 6:
- dset[cc*batch_size:cc*batch_size+nn, 3:] = get_random_versors(nn)
+ if self.dns_type == 'NSVE_Stokes_particles':
+ dset[cc*batch_size:cc*batch_size+nn, 3:] = self.parameters['initial_particle_vel']*get_random_versors(nn)
+ else:
+ dset[cc*batch_size:cc*batch_size+nn, 3:] = get_random_versors(nn)
nn = 0
cc += 1
except Exception as e:
@@ -1036,6 +1069,8 @@ class DNS(_code):
if self.dns_type in ['NSVEcomplex_particles']:
particle_file.create_group('tracers0/orientation')
particle_file.create_group('tracers0/velocity_gradient')
+ if self.dns_type in ['NSVE_Stokes_particles']:
+ particle_file.create_group('tracers0/momentum')
if self.dns_type in ['NSVEp_extra_sampling']:
particle_file.create_group('tracers0/velocity_gradient')
particle_file.create_group('tracers0/pressure')
@@ -1049,6 +1084,16 @@ class DNS(_code):
# first, check if initial field exists
need_field = False
if self.check_current_vorticity_exists:
+ need_field = True
+ if self.dns_type in [
+ 'NSVE',
+ 'NSVE_no_output',
+ 'static_field',
+ 'NSVEparticles',
+ 'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
+ 'NSVEparticles_no_output',
+ 'NSVEp_extra_sampling']:
if not os.path.exists(self.get_checkpoint_0_fname()):
need_field = True
else:
@@ -1082,10 +1127,16 @@ class DNS(_code):
f,
'vorticity/complex/{0}'.format(0))
else:
- data = self.generate_vector_field(
+ if self.dns_type == 'NSVE_Stokes_particles':
+ data = self.generate_vector_field(
+ write_to_file = False,
+ spectra_slope = 2.0,
+ amplitude = self.parameters['initial_field_amplitude'])
+ else:
+ data = self.generate_vector_field(
write_to_file = False,
spectra_slope = 2.0,
- amplitude = 0.05)
+ amplitude = 0.05)
f['vorticity/complex/{0}'.format(0)] = data
f.close()
if self.dns_type == 'kraichnan_field':
@@ -1100,6 +1151,7 @@ class DNS(_code):
'static_field_with_ghost_collisions',
'NSVEparticles',
'NSVEcomplex_particles',
+ 'NSVE_Stokes_particles',
'NSVEparticles_no_output',
'NSVEp_extra_sampling']:
self.generate_particle_data(opt = opt)
diff --git a/TurTLE/PP.py b/TurTLE/PP.py
index 236b4963ab05f69d49fb6367d598ac83f02d51e2..acf3b003c0881f89fbed97fffb67c897a6b00786 100644
--- a/TurTLE/PP.py
+++ b/TurTLE/PP.py
@@ -145,6 +145,8 @@ class PP(_code):
pars['filter_type'] = 'Gauss'
pars['max_velocity_estimate'] = float(10)
pars['histogram_bins'] = int(129)
+ elif dns_type == 'get_rfields':
+ pars['TrS2_on'] = int(0)
return pars
def get_data_file_name(self):
return os.path.join(self.work_dir, self.simname + '.h5')
@@ -447,13 +449,17 @@ class PP(_code):
parser_resize = subparsers.add_parser(
'resize',
help = 'get joint acceleration and velocity statistics')
+ parser_write_rpressure = subparsers.add_parser(
+ 'write_rpressure',
+ help = 'write real pressure field to binary')
for pp_type in [
'resize',
'joint_acc_vel_stats',
'bandpass_stats',
'get_rfields',
'field_single_to_double',
- 'native_binary_to_hdf5']:
+ 'native_binary_to_hdf5',
+ 'write_rpressure']:
eval('self.simulation_parser_arguments(parser_' + pp_type + ')')
eval('self.job_parser_arguments(parser_' + pp_type + ')')
eval('self.parameters_to_parser_arguments(parser_' + pp_type + ')')
@@ -805,15 +811,21 @@ class PP(_code):
with h5py.File(os.path.join(self.work_dir, self.simname + '_fields.h5'), 'a') as ff:
ff.require_group('vorticity')
ff.require_group('vorticity/complex')
- checkpoint_file_list = glob.glob(self.simname + '_checkpoint_*.h5')
+ checkpoint_file_list = [self.simname + '_checkpoint_{0}.h5'.format(cp)
+ for cp in range(df['checkpoint'][()]+1)]
for cpf_name in checkpoint_file_list:
- cpf = h5py.File(cpf_name, 'r')
- 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)
- cpf.close()
+ if os.path.exists(cpf_name):
+ 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)
+ cpf.close()
return None
def launch_jobs(
self,
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 4bf762e7ab9ff1159929fcee46a9247a3c2dff44..59757f6461416a31a4afe5936ecf21ab571a7f26 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -188,6 +188,7 @@ class TEST(_code):
scal_rspace_stats = []
if self.dns_type in ['Gauss_field_test']:
vec_spectra_stats.append('velocity')
+ vec_spectra_stats.append('k*velocity')
vec4_rspace_stats.append('velocity')
tens_rspace_stats.append('velocity_gradient')
scal_rspace_stats.append('velocity_divergence')
@@ -377,6 +378,7 @@ class TEST(_code):
eval('self.simulation_parser_arguments(parser_' + parser + ')')
eval('self.job_parser_arguments(parser_' + parser + ')')
eval('self.parameters_to_parser_arguments(parser_' + parser + ')')
+ eval('self.parameters_to_parser_arguments(parser_' + parser + ', self.generate_extra_parameters(dns_type = \'' + parser + '\'))')
return None
def prepare_launch(
self,
diff --git a/TurTLE/_code.py b/TurTLE/_code.py
index be898972439e485f3ad3f6b3e4a947c28bc14951..b49c8f2edc5e256912a4e357f17bfeb627a90c67 100644
--- a/TurTLE/_code.py
+++ b/TurTLE/_code.py
@@ -4,18 +4,18 @@
# #
# This file is part of bfps. #
# #
-# TurTLE is free software: you can redistribute it and/or modify #
+# TurTLE is free software: you can redistribute it and/or modify #
# it under the terms of the GNU General Public License as published #
# by the Free Software Foundation, either version 3 of the License, #
# or (at your option) any later version. #
# #
-# TurTLE is distributed in the hope that it will be useful, #
+# TurTLE is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
-# along with TurTLE. If not, see <http://www.gnu.org/licenses/> #
+# along with TurTLE. If not, see <http://www.gnu.org/licenses/> #
# #
# Contact: Cristian.Lalescu@ds.mpg.de #
# #
@@ -678,9 +678,9 @@ class _code(_base):
script_file.write('#SBATCH -o ' + out_file + '\n')
# set up environment
- script_file.write('#SBATCH --get-user-env\n')
- script_file.write('#SBATCH --partition={0}\n'.format(
- self.host_info['environment']))
+ if self.host_info['explicit_slurm_environment']:
+ script_file.write('#SBATCH --partition={0}\n'.format(
+ self.host_info['environment']))
if 'account' in self.host_info.keys():
script_file.write('#SBATCH --account={0}\n'.format(
self.host_info['account']))
@@ -717,14 +717,19 @@ class _code(_base):
script_file.write('#SBATCH --mail-type=none\n')
script_file.write('#SBATCH --time={0}:{1:0>2d}:00\n'.format(hours, minutes))
+ if 'extra_slurm_lines' in self.host_info.keys():
+ for line in self.host_info['extra_slurm_lines']:
+ script_file.write(line + '\n')
## following cleans up environment for job.
+ ## put these in the 'extra_slurm_lines' list if you need them
## make sure that "~/.config/TurTLE/bashrc" exists and builds desired job environment
- script_file.write('#SBATCH --export=NONE\n')
- script_file.write('#SBATCH --get-user-env\n')
- script_file.write('source ~/.config/TurTLE/bashrc\n')
+ #script_file.write('#SBATCH --export=NONE\n')
+ #script_file.write('#SBATCH --get-user-env\n')
+ #script_file.write('export OMP_PLACES=cores\n') # or threads, as appropriate
+ # also look up OMP_PROC_BIND and SLURM_HINT
+ #script_file.write('source ~/.config/TurTLE/bashrc\n')
if nb_threads_per_process > 1:
- script_file.write('export OMP_NUM_THREADS={0}\n'.format(nb_threads_per_process))
- script_file.write('export OMP_PLACES=cores\n')
+ script_file.write('export OMP_NUM_THREADS=${SLURM_CPUS_PER_TASK}\n')
script_file.write('echo "Start time is `date`"\n')
script_file.write('cd ' + self.work_dir + '\n')
diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
index 149073e7c1abe931368df1a026ef8f02ac4bc070..00b1692e084cbb346a3a19d3f046cdba5cb1d01b 100644
--- a/TurTLE/test/test_Gaussian_field.py
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -5,7 +5,7 @@ from scipy import trapz
from scipy.stats import norm
from scipy.integrate import quad
import h5py
-import sys
+import sys, os
import time
import TurTLE
@@ -17,24 +17,28 @@ except:
plt = None
def main():
- c = TEST()
# size of grid
- n = 256
+ n = 1024
slope = -5./3.
- k_cutoff = 30.
+ k_cutoff = 64.
func = lambda k, k_c=k_cutoff, s=slope : k**s*np.exp(-k/k_c)
- total_energy = quad(func, 1, k_cutoff*4)[0]
+ total_energy = quad(func, 0.6, k_cutoff*8)[0]
coeff = 1./total_energy
bin_no = 100
rseed = int(time.time())
+ simname = 'Gaussianity_test'
- c.launch(
+ if not os.path.exists(simname+'.h5'):
+ c = TEST()
+ opt = c.launch(
['Gauss_field_test',
'--nx', str(n),
'--ny', str(n),
'--nz', str(n),
- '--simname', 'Gaussianity_test',
+ '--simname', simname,
'--np', '4',
+ '--environment', 'short',
+ '--minutes', '60',
'--ntpp', '1',
'--wd', './',
'--histogram_bins', str(bin_no),
@@ -44,7 +48,7 @@ def main():
'--spectrum_coefficient', str(coeff),
'--field_random_seed', str(rseed)] +
sys.argv[1:])
- plot_stuff(c.simname, total_energy = total_energy)
+ plot_stuff(simname, total_energy = total_energy)
return None
def plot_stuff(simname, total_energy = 1.):
@@ -82,6 +86,7 @@ def plot_stuff(simname, total_energy = 1.):
f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
print('Energy analytically: {}'.format(total_energy))
+ print(np.sum(energy*np.arange(len(energy))**2))
print('Energy sum: {}'.format(np.sum(energy*df['kspace/dk'][()])))
print('Moment sum: {}'.format(df['statistics/moments/velocity'][0,2,3]/2))
print('Velocity variances: {}'.format(trapz(vel[:,None]**2*f_vel, vel[:,None], axis=0)))
@@ -101,6 +106,17 @@ def plot_stuff(simname, total_energy = 1.):
df['statistics/moments/velocity_divergence'][0, 2]))
print('Gradient second moment is: {0}'.format(
df['statistics/moments/velocity_gradient'][0, 2].mean()))
+
+ print('----------- k2-premultiplied spectrum -----------')
+ k2func = lambda k, k_c=k_cutoff, s=slope : k**(2+s)*np.exp(-k/k_c)
+ k2sum_analytic = quad(k2func, 0, k_cutoff*20)[0]
+ print('Analytically: {}'.format(k2sum_analytic))
+ k2spec_trace = (
+ df['statistics/spectra/k*velocity_k*velocity'][..., 0, 0]
+ + df['statistics/spectra/k*velocity_k*velocity'][..., 1, 1]
+ + df['statistics/spectra/k*velocity_k*velocity'][..., 2, 2])
+ print('Energy sum: {}'.format(np.sum(k2spec_trace*df['kspace/dk'][()])/2./coeff))
+
df.close()
return None
diff --git a/TurTLE/test/test_particle_clouds.py b/TurTLE/test/test_particle_clouds.py
index f72071143778bb79cbe2ea1128769c9c32004f5e..b6c401d18039d2b5b02991bd5a1cf1aae8152a50 100644
--- a/TurTLE/test/test_particle_clouds.py
+++ b/TurTLE/test/test_particle_clouds.py
@@ -33,22 +33,22 @@ import sys
import TurTLE
from TurTLE import DNS
-
-def main():
+def basic_test():
nclouds = 10
nparticles_per_cloud = 1000
nparticles = nclouds*nparticles_per_cloud
niterations = 32
c = DNS()
c.dns_type = 'NSVEparticles'
- c.parameters['nparticles'] = nparticles
- c.parameters['tracers1_integration_steps'] = 4
- c.generate_tracer_state(rseed = 2, species = 1)
- del c.parameters['nparticles']
- del c.parameters['tracers1_integration_steps']
+ c.simname = 'basic_cloud_test'
+ f0 = h5py.File(
+ os.path.join(
+ os.path.join(TurTLE.lib_dir, 'test'),
+ 'B32p1e4_checkpoint_0.h5'),
+ 'r')
ic_file = h5py.File(c.get_checkpoint_0_fname(), 'a')
- ic_file['tracers0/state/0'] = ic_file['tracers1/state/0'][...].reshape(nclouds, nparticles_per_cloud, 3)
- ic_file['tracers0/rhs/0'] = ic_file['tracers1/rhs/0'][...].reshape(4, nclouds, nparticles_per_cloud, 3)
+ ic_file['tracers0/state/0'] = f0['tracers0/state/0'][...].reshape(nclouds, nparticles_per_cloud, 3)
+ ic_file['tracers0/rhs/0'] = f0['tracers0/rhs/0'][...].reshape(4, nclouds, nparticles_per_cloud, 3)
ic_file.close()
c.launch(
['NSVEparticles',
@@ -57,12 +57,14 @@ def main():
'--forcing_type', 'linear',
'--src-wd', TurTLE.lib_dir + '/test',
'--src-iteration', '0',
+ '--simname', c.simname,
'--np', '4',
'--ntpp', '1',
'--fftw_plan_rigor', 'FFTW_PATIENT',
'--niter_todo', '{0}'.format(niterations),
'--niter_out', '{0}'.format(niterations),
'--niter_stat', '1',
+ '--checkpoints_per_file', '{0}'.format(3),
'--nparticles', '{0}'.format(nparticles),
'--njobs', '2',
'--wd', './'])
@@ -79,6 +81,7 @@ def main():
x0 = f0['tracers0/state/{0}'.format(iteration)][...]
x1 = f1['tracers0/state/{0}'.format(iteration)][...].reshape(x0.shape)
traj_error = np.max(np.abs(x0 - x1))
+ print(traj_error)
y0 = f0['tracers0/rhs/{0}'.format(iteration)][...]
y1 = f1['tracers0/rhs/{0}'.format(iteration)][...].reshape(y0.shape)
rhs_error = np.max(np.abs(y0 - y1))
@@ -88,6 +91,53 @@ def main():
print('SUCCESS! Basic test passed.')
return None
+def nasty_test():
+ nclouds = 10
+ nparticles_per_cloud = 1000000
+ nparticles = nclouds*nparticles_per_cloud
+ niterations = 8
+ c = DNS()
+ c.dns_type = 'NSVEparticles'
+ c.simname = 'nasty_cloud_test'
+ c.parameters['nparticles'] = nparticles
+ c.parameters['tracers1_integration_steps'] = 4
+ c.generate_tracer_state(rseed = 2, species = 1)
+ del c.parameters['nparticles']
+ del c.parameters['tracers1_integration_steps']
+ ic_file = h5py.File(c.get_checkpoint_0_fname(), 'a')
+ ic_file['tracers0/state/0'] = ic_file['tracers1/state/0'][...].reshape(nclouds, nparticles_per_cloud, 3)
+ ic_file['tracers0/rhs/0'] = ic_file['tracers1/rhs/0'][...].reshape(4, nclouds, nparticles_per_cloud, 3)
+ # put all particles in single z cell
+ ic_file['tracers0/state/0'][..., 2] = 0.0001
+ # put one cloud on another process
+ ic_file['tracers0/state/0'][1, :, 2] = np.pi + 0.0001
+ ic_file.close()
+ c.launch(
+ ['NSVEparticles',
+ '-n', '8',
+ '--src-simname', 'B32p1e4',
+ '--simname', c.simname,
+ '--forcing_type', 'linear',
+ '--src-wd', TurTLE.lib_dir + '/test',
+ '--src-iteration', '0',
+ '--np', '4',
+ '--ntpp', '2',
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(niterations),
+ '--niter_stat', '1',
+ '--nparticles', '{0}'.format(nparticles),
+ '--njobs', '1',
+ '--wd', './'])
+ print('SUCCESS! Nasty test passed.')
+ return None
+
+
+def main():
+ basic_test()
+ nasty_test()
+ return None
+
if __name__ == '__main__':
main()
diff --git a/TurTLE/test/test_turtle_NSVE_Stokes_particles.py b/TurTLE/test/test_turtle_NSVE_Stokes_particles.py
new file mode 100644
index 0000000000000000000000000000000000000000..400a0cdb12ba0790343489a8570f9383a16daa80
--- /dev/null
+++ b/TurTLE/test/test_turtle_NSVE_Stokes_particles.py
@@ -0,0 +1,100 @@
+#! /usr/bin/env python
+#######################################################################
+# #
+# Copyright 2019 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@ds.mpg.de #
+# #
+#######################################################################
+
+
+
+import os
+import numpy as np
+import h5py
+import sys
+import matplotlib.pyplot as plt
+
+
+import TurTLE
+from TurTLE import DNS
+
+def theory_exponential(t,initial_value,drag_coeff):
+ return initial_value*np.exp(-t*drag_coeff)
+
+def main():
+ """
+ Run test where the trajectory of initially moving particles with varying drag coefficient in a
+ quiescent flow is compared to the analytically expected exponential. The generated plot shows the accuracy of the corresponding particles for a certain fixed timestep.
+ """
+ drag_coefficients = ['30','20','10','1']
+ niterations = 100
+ nparticles = 1
+ njobs = 1
+ dt=0.01
+ fig = plt.figure()
+ plt.ylabel('$v$')
+ plt.xlabel('$t$')
+ plt.yscale('log')
+ j=0
+ for drag_coeff in drag_coefficients:
+ c = DNS()
+ c.launch(
+ ['NSVE_Stokes_particles',
+ '-n', '32',
+ '--simname', 'quiescent_nsve_stokes_particles_drag'+str(j),
+ '--np', '4',
+ '--ntpp', '1',
+ '--dt', '{0}'.format(dt),
+ '--fftw_plan_rigor', 'FFTW_PATIENT',
+ '--niter_todo', '{0}'.format(niterations),
+ '--niter_out', '{0}'.format(1),
+ '--niter_stat', '1',
+ '--nparticles', '{0}'.format(nparticles),
+ '--initial_field_amplitude', '0',
+ '--initial_particle_vel', '0.05',
+ '--drag_coefficient', drag_coeff,
+ '--famplitude', '0',
+ '--njobs', '{0}'.format(njobs),
+ '--wd', './'] +
+ sys.argv[1:])
+ f = h5py.File('quiescent_nsve_stokes_particles_drag'+str(j)+'_particles.h5', 'r')
+ f2 = h5py.File('quiescent_nsve_stokes_particles_drag'+str(j)+'.h5', 'r')
+ particle_momentum = [np.sqrt(f['tracers0/momentum/'+'{}'.format(i)][0][0]**2
+ +f['tracers0/momentum/'+'{}'.format(i)][0][1]**2
+ +f['tracers0/momentum/'+'{}'.format(i)][0][2]**2)
+ for i in range(niterations)]
+ time_values = np.arange(0, niterations*dt, dt)
+ numerics, = plt.plot(time_values, particle_momentum, alpha = 1, lw=2)
+ plt.plot(time_values,
+ theory_exponential(time_values, f2['parameters/initial_particle_vel'],np.float(drag_coeff)),
+ color='black', ls = '--')
+ plt.plot(0,0,color=numerics.get_color(),label=r'drag coefficient '+drag_coeff)
+ j+=1
+ plt.plot(0,
+ 0,
+ color='black', ls = '--',
+ label='analytic sol.')
+ plt.legend(loc='lower left')
+ plt.savefig('Stokes_quiescent_test.pdf')
+ return None
+
+if __name__ == '__main__':
+ main()
+
diff --git a/TurTLE/test/test_turtle_NSVEparticles.py b/TurTLE/test/test_turtle_NSVEparticles.py
index 0e600e7847fe6c01adc4578682f6aebfff80ac6f..8329a0b5ceebdb0f8477fbb45dc05645891a774e 100644
--- a/TurTLE/test/test_turtle_NSVEparticles.py
+++ b/TurTLE/test/test_turtle_NSVEparticles.py
@@ -4,20 +4,20 @@
# Copyright 2019 Max Planck Institute #
# for Dynamics and Self-Organization #
# #
-# This file is part of TurTLE. #
+# This file is part of TurTLE. #
# #
-# TurTLE is free software: you can redistribute it and/or modify #
+# TurTLE is free software: you can redistribute it and/or modify #
# it under the terms of the GNU General Public License as published #
# by the Free Software Foundation, either version 3 of the License, #
# or (at your option) any later version. #
# #
-# TurTLE is distributed in the hope that it will be useful, #
+# TurTLE is distributed in the hope that it will be useful, #
# but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. #
# #
# You should have received a copy of the GNU General Public License #
-# along with TurTLE. If not, see <http://www.gnu.org/licenses/> #
+# along with TurTLE. If not, see <http://www.gnu.org/licenses/> #
# #
# Contact: Cristian.Lalescu@ds.mpg.de #
# #
diff --git a/cpp/Lagrange_polys.cpp b/cpp/Lagrange_polys.cpp
index bbc8997f896aebba486e4fa601c18eaf2aacfcc1..eff66cd40cb5c8dd70a6813a4b5882e52893a467 100644
--- a/cpp/Lagrange_polys.cpp
+++ b/cpp/Lagrange_polys.cpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/Lagrange_polys.hpp b/cpp/Lagrange_polys.hpp
index 9f4742a3303d9586cc7832e568162465767a4922..7d4d27f8ac2868c0acab48e07e9af23d09b2423d 100644
--- a/cpp/Lagrange_polys.hpp
+++ b/cpp/Lagrange_polys.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/base.hpp b/cpp/base.hpp
index 5e53f4451689e3bf3bcb1f4902e83e3b5276d818..60fbb4168b5a091ceac476d41737c4d56543fc8a 100644
--- a/cpp/base.hpp
+++ b/cpp/base.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/bfps_timer.hpp b/cpp/bfps_timer.hpp
index b299dc115555bd3952e22ce5984f27eb9bfb1914..d5b1147ba472ddd7ba3c61cadad0f1d20188e6e1 100644
--- a/cpp/bfps_timer.hpp
+++ b/cpp/bfps_timer.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/env_utils.hpp b/cpp/env_utils.hpp
similarity index 93%
rename from cpp/particles/env_utils.hpp
rename to cpp/env_utils.hpp
index 829fd5b46f879c4485276d3f3866b8ae3d81e8d5..7530c2a6a11fc6292de57de3d7411fb5d8d9a687 100644
--- a/cpp/particles/env_utils.hpp
+++ b/cpp/env_utils.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/fftw_interface.hpp b/cpp/fftw_interface.hpp
index 0a840dd5ba3d864b36271515faa7cb81f3042c01..fd88ed4aff49285ff907a2f5ddcd01af42d42612 100644
--- a/cpp/fftw_interface.hpp
+++ b/cpp/fftw_interface.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/fftw_tools.cpp b/cpp/fftw_tools.cpp
index e28d86c0a4f8eba08faf979bcb52eced4efcd972..ae03966d686baab544d865cf1eab253faea6584d 100644
--- a/cpp/fftw_tools.cpp
+++ b/cpp/fftw_tools.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/fftw_tools.hpp b/cpp/fftw_tools.hpp
index b41cd2a453c2c0aa34f56febb17f2a650a2a9685..18cca206a3c2ffa32329c7e7ba3889ff7459569e 100644
--- a/cpp/fftw_tools.hpp
+++ b/cpp/fftw_tools.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 81f71e10d330ea0886091d13348584378e01011c..3ee9f2995dd1cbc2dea28abdc3d1b675d2cc6031 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -799,8 +799,8 @@ int field<rnumber, be, fc>::write_filtered(
}
}
}
- DEBUG_MSG("count[0] = %ld, offset[0] = %ld\n",
- count[0], offset[0]);
+ //DEBUG_MSG("count[0] = %ld, offset[0] = %ld\n",
+ // count[0], offset[0]);
if (nz>=2)
{
assert(nz%2==0);
@@ -812,8 +812,8 @@ 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]);
+ //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);
@@ -830,7 +830,7 @@ 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");
+ //DEBUG_MSG("filtered write nz=1\n");
//now write data
mspace = H5Screate_simple(ndim(fc), memshape, NULL);
@@ -872,12 +872,12 @@ void field<rnumber, be, fc>::compute_rspace_xincrement_stats(
this->rlayout->sizes[0],
this->rlayout->comm);
tmp_field->real_space_representation = true;
- this->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
- hsize_t rrindex = (xindex + xcells)%this->rlayout->sizes[2] + (
+ this->RLOOP_simd(
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ const hsize_t rrindex = (xindex + xcells)%this->rlayout->sizes[2] + (
zindex * this->rlayout->subsizes[1] + yindex)*(
this->rmemlayout->subsizes[2]);
for (unsigned int component=0; component < ncomp(fc); component++)
@@ -921,9 +921,9 @@ void field<rnumber, be, fc>::compute_rspace_stats(
// and possible component indices that are NOT averaged over
// HDF5 dataset has 2 extra indices for time and order of moment,
// and possible component indices
- DEBUG_MSG("ndims = %d, ndim(fc) = %d, dsetname = %s\n",
- ndims, ndim(fc),
- dset_name.c_str());
+ //DEBUG_MSG("ndims = %d, ndim(fc) = %d, dsetname = %s\n",
+ // ndims, ndim(fc),
+ // dset_name.c_str());
assert(ndims == int(ndim(fc))-1);
assert(dims[1] == nmoments);
switch(ndims)
@@ -981,10 +981,10 @@ void field<rnumber, be, fc>::compute_rspace_stats(
{
TIMEZONE("field::RLOOP");
this->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
double* pow_tmp = local_pow_tmp.getMine();
std::fill_n(pow_tmp, nvals, 1);
@@ -1005,7 +1005,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];
- int bin = int(floor(val_tmp[3]*2/binsize[3]));
+ const int bin = int(floor(val_tmp[3]*2/binsize[3]));
if (bin >= 0 && bin < nbins){
local_hist[bin*nvals+3]++;
}
@@ -1016,7 +1016,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];
- int bin = int(floor((val_tmp[i] + max_estimate[i]) / binsize[i]));
+ const int bin = int(floor((val_tmp[i] + max_estimate[i]) / binsize[i]));
if (bin >= 0 && bin < nbins)
local_hist[bin*nvals+i]++;
}
@@ -1146,10 +1146,10 @@ void field<rnumber, be, fc>::compute_rspace_zaverage(
{
TIMEZONE("field::RLOOP");
this->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
double *local_zaverage = local_zaverage_threaded.getMine();
ptrdiff_t zaverage_index = (yindex*this->rlayout->subsizes[2]+xindex)*ncomp(fc);
@@ -1467,10 +1467,10 @@ double field<rnumber, be, fc>::L2norm(
std::fill_n(local_moment, 1, 0);});
this->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
double *local_m2 = local_m2_threaded.getMine();
for (unsigned int i=0; i<ncomp(fc); i++)
local_m2[0] += this->data[rindex*ncomp(fc)+i]*this->data[rindex*ncomp(fc)+i];
@@ -1507,11 +1507,11 @@ int compute_gradient(
{
case ONE:
kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 < kk->kM2)
{
dst->cval(cindex, 0, 0) = -kk->kx[xindex]*src->cval(cindex, 1);
@@ -1524,11 +1524,11 @@ int compute_gradient(
break;
case THREE:
kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 < kk->kM2)
{
for (unsigned int field_component = 0;
@@ -1561,11 +1561,11 @@ int compute_divergence(
*dst = 0.0;
dst->real_space_representation = false;
kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
dst->cval(cindex, 0) = -(kk->kx[xindex]*src->cval(cindex, 0, 1) +
kk->ky[yindex]*src->cval(cindex, 1, 1) +
kk->kz[zindex]*src->cval(cindex, 2, 1));
@@ -1589,24 +1589,27 @@ int invert_curl(
std::fill_n(dst->get_rdata(), dst->rmemlayout->local_size, 0);
dst->real_space_representation = false;
kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
- if (k2 <= kk->kM2 && k2 > 0)
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ const double kk2 = (k2 > 0) ? k2 : 1.0;
+ if (k2 <= kk->kM2)
{
- dst->cval(cindex,0,0) = -(kk->ky[yindex]*src->cval(cindex,2,1) - kk->kz[zindex]*src->cval(cindex,1,1)) / k2;
- dst->cval(cindex,0,1) = (kk->ky[yindex]*src->cval(cindex,2,0) - kk->kz[zindex]*src->cval(cindex,1,0)) / k2;
- dst->cval(cindex,1,0) = -(kk->kz[zindex]*src->cval(cindex,0,1) - kk->kx[xindex]*src->cval(cindex,2,1)) / k2;
- dst->cval(cindex,1,1) = (kk->kz[zindex]*src->cval(cindex,0,0) - kk->kx[xindex]*src->cval(cindex,2,0)) / k2;
- dst->cval(cindex,2,0) = -(kk->kx[xindex]*src->cval(cindex,1,1) - kk->ky[yindex]*src->cval(cindex,0,1)) / k2;
- dst->cval(cindex,2,1) = (kk->kx[xindex]*src->cval(cindex,1,0) - kk->ky[yindex]*src->cval(cindex,0,0)) / k2;
+ dst->cval(cindex,0,0) = -(kk->ky[yindex]*src->cval(cindex,2,1) - kk->kz[zindex]*src->cval(cindex,1,1)) / kk2;
+ dst->cval(cindex,0,1) = (kk->ky[yindex]*src->cval(cindex,2,0) - kk->kz[zindex]*src->cval(cindex,1,0)) / kk2;
+ dst->cval(cindex,1,0) = -(kk->kz[zindex]*src->cval(cindex,0,1) - kk->kx[xindex]*src->cval(cindex,2,1)) / kk2;
+ dst->cval(cindex,1,1) = (kk->kz[zindex]*src->cval(cindex,0,0) - kk->kx[xindex]*src->cval(cindex,2,0)) / kk2;
+ dst->cval(cindex,2,0) = -(kk->kx[xindex]*src->cval(cindex,1,1) - kk->ky[yindex]*src->cval(cindex,0,1)) / kk2;
+ dst->cval(cindex,2,1) = (kk->kx[xindex]*src->cval(cindex,1,0) - kk->ky[yindex]*src->cval(cindex,0,0)) / kk2;
}
else
std::fill_n((rnumber*)(dst->get_cdata()+3*cindex), 6, 0.0);
}
);
+ if (kk->layout->myrank == 0)
+ std::fill_n((rnumber*)(dst->get_cdata()), 6, 0.0);
dst->symmetrize();
return EXIT_SUCCESS;
}
@@ -1720,10 +1723,10 @@ int joint_rspace_PDF(
{
TIMEZONE("field::RLOOP");
f1->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
ptrdiff_t *local_histm = local_histm_threaded.getMine();
int bin1 = 0;
int bin2 = 0;
@@ -1908,10 +1911,10 @@ int joint_rspace_3PDF(
{
TIMEZONE("field::RLOOP");
f1->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
ptrdiff_t *local_histm = local_histm_threaded.getMine();
int bin1 = 0;
int bin2 = 0;
@@ -2130,19 +2133,19 @@ int make_gaussian_random_field(
rseed*omp_get_max_threads()*output_field->clayout->nprocs +
output_field->clayout->myrank*omp_get_max_threads() +
thread_id);
- DEBUG_MSG("in make_gaussian_random_field, thread_id = %d, current_seed = %d\n", thread_id, current_seed);
+ //DEBUG_MSG("in make_gaussian_random_field, thread_id = %d, current_seed = %d\n", thread_id, current_seed);
rgen[thread_id].seed(current_seed);
}
output_field->real_space_representation = false;
*output_field = 0.0;
- DEBUG_MSG("slope: %g\n", slope);
+ //DEBUG_MSG("slope: %g\n", slope);
// inside loop use only thread-local random number generator
kk->CLOOP_K2([&](
- ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2)
+ const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2)
{
if (k2 > 0)
switch(fc)
diff --git a/cpp/field.hpp b/cpp/field.hpp
index b1dbc5e03e8d9c8b8a306518956993adec4abe7e..0364e8d527a9f4ae8b6c30b7026b10c5245b1038 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -242,9 +242,14 @@ class field
inline field<rnumber, be, fc>& operator=(const rnumber value)
{
- std::fill_n(this->data,
- this->rmemlayout->local_size,
- value);
+ #pragma omp parallel
+ {
+ const hsize_t start = OmpUtils::ForIntervalStart(this->rmemlayout->local_size);
+ const hsize_t end = OmpUtils::ForIntervalEnd(this->rmemlayout->local_size);
+ std::fill_n(this->data + start,
+ end - start,
+ value);
+ }
return *this;
}
@@ -269,6 +274,33 @@ class field
}
}
template <class func_type>
+ void RLOOP_simd(func_type expression)
+ {
+ 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]);
+
+ #pragma omp simd
+ for (hsize_t zindex = 0; zindex < this->rlayout->subsizes[0]; zindex++)
+ for (hsize_t yindex = start; yindex < end; yindex++)
+ {
+ const ptrdiff_t rindex = (
+ zindex * this->rlayout->subsizes[1] + yindex)*(
+ this->rmemlayout->subsizes[2]);
+ for (hsize_t xindex = 0; xindex < this->rlayout->subsizes[2]; xindex++)
+ {
+ expression(rindex + xindex, xindex, yindex, zindex);
+ }
+ }
+ }
+ break;
+ }
+ }
+ template <class func_type>
void RLOOP(func_type expression)
{
switch(be)
@@ -280,15 +312,15 @@ class field
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++)
{
- ptrdiff_t rindex = (
+ const ptrdiff_t rindex = (
zindex * this->rlayout->subsizes[1] + yindex)*(
this->rmemlayout->subsizes[2]);
for (hsize_t xindex = 0; xindex < this->rlayout->subsizes[2]; xindex++)
{
- expression(rindex, xindex, yindex, zindex);
- rindex++;
+ expression(rindex + xindex, xindex, yindex, zindex);
}
}
}
diff --git a/cpp/field_binary_IO.cpp b/cpp/field_binary_IO.cpp
index 52ca21f5947689408265423e0c4f075a670fa578..288ef26a9cdb490895e908c5c3f377ed13c756de 100644
--- a/cpp/field_binary_IO.cpp
+++ b/cpp/field_binary_IO.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/field_binary_IO.hpp b/cpp/field_binary_IO.hpp
index b06e901c27acd99855a314dc7b5afb8b43f190b7..66234ab39f7e0a1172c1c80ad464c33ab4a8d8e5 100644
--- a/cpp/field_binary_IO.hpp
+++ b/cpp/field_binary_IO.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/field_layout.cpp b/cpp/field_layout.cpp
index e2872a64fd157532bf60d036fc2dc71797f10e71..0f19a6cd98a5dbe6d2be54d80f887d77a33f67ea 100644
--- a/cpp/field_layout.cpp
+++ b/cpp/field_layout.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/field_layout.hpp b/cpp/field_layout.hpp
index 770119c2dcb05017d495b62559f050646872dc84..5029ae30af83dcf65796d0f717b209e7816ae5d6 100644
--- a/cpp/field_layout.hpp
+++ b/cpp/field_layout.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
index 2ac29889aaaad33895ecf145c62610feee952f2b..744832c4aff23f6f2a65797af58429ed3e579679 100644
--- a/cpp/full_code/Gauss_field_test.cpp
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -145,6 +145,14 @@ int Gauss_field_test<rnumber>::do_work(void)
this->vector_field,
this->scalar_field);
+ /// compute integral of premultiplied energy spectrum
+ this->kk->template cospectrum<rnumber, THREE>(
+ this->vector_field->get_cdata(),
+ stat_group,
+ "k*velocity_k*velocity",
+ 0,
+ 2.);
+
/// compute basic statistics of Gaussian field
this->vector_field->compute_stats(
this->kk,
diff --git a/cpp/full_code/NSVE_Stokes_particles.cpp b/cpp/full_code/NSVE_Stokes_particles.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..a4b3c41e8908a853333d6041b073d8742021b87a
--- /dev/null
+++ b/cpp/full_code/NSVE_Stokes_particles.cpp
@@ -0,0 +1,214 @@
+/**********************************************************************
+* *
+* 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 "NSVE_Stokes_particles.hpp"
+#include "scope_timer.hpp"
+#include "particles/particles_sampling.hpp"
+#include "particles/p2p_ghost_collisions.hpp"
+#include "particles/particles_inner_computer_2nd_order.hpp"
+
+template <typename rnumber>
+int NSVE_Stokes_particles<rnumber>::initialize(void)
+{
+ TIMEZONE("NSVE_Stokes_particles::intialize");
+ this->NSVE<rnumber>::initialize();
+ this->pressure = new field<rnumber, FFTW, ONE>(
+ this->fs->cvelocity->rlayout->sizes[2],
+ this->fs->cvelocity->rlayout->sizes[1],
+ this->fs->cvelocity->rlayout->sizes[0],
+ this->fs->cvelocity->rlayout->comm,
+ this->fs->cvelocity->fftw_plan_rigor);
+
+ particles_inner_computer_2nd_order_Stokes<double, long long int> current_particles_inner_computer;
+ current_particles_inner_computer.set_drag_coefficient(this->drag_coefficient);
+ this->ps = particles_system_builder_with_p2p(
+ this->fs->cvelocity, // (field object)
+ this->fs->kk, // (kspace object, contains dkx, dky, dkz)
+ tracers0_integration_steps, // to check coherency between parameters and hdf input file (nb rhs)
+ (long long int)nparticles, // to check coherency between parameters and hdf input file
+ this->fs->get_current_fname(), // particles input filename
+ std::string("/tracers0/state/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ std::string("/tracers0/rhs/") + std::to_string(this->fs->iteration), // dataset name for initial input
+ tracers0_neighbours, // parameter (interpolation no neighbours)
+ tracers0_smoothness, // parameter
+ this->comm,
+ this->fs->iteration+1,
+ std::move(p2p_ghost_collisions<double, long long int>()),
+ std::move(current_particles_inner_computer),
+ this->tracers0_cutoff);
+ //DEBUG_MSG_WAIT(MPI_COMM_WORLD, "after call to particles_system_builder\n");
+ this->particles_output_writer_mpi = new particles_output_hdf5<
+ long long int, double, 6>(
+ MPI_COMM_WORLD,
+ "tracers0",
+ nparticles,
+ tracers0_integration_steps);
+ this->particles_output_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
+ this->particles_sample_writer_mpi = new particles_output_sampling_hdf5<
+ long long int, double, 3>(
+ MPI_COMM_WORLD,
+ this->ps->getGlobalNbParticles(),
+ (this->simname + "_particles.h5"),
+ "tracers0",
+ "position/0");
+ this->particles_sample_writer_mpi->setParticleFileLayout(this->ps->getParticleFileLayout());
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_particles<rnumber>::step(void)
+{
+ TIMEZONE("NSVE_Stokes_particles::step");
+ this->fs->compute_velocity(this->fs->cvorticity);
+ this->fs->cvelocity->ift();
+ this->ps->complete2ndOrderLoop(this->dt);
+ this->NSVE<rnumber>::step();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_particles<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("NSVE_Stokes_particles::write_checkpoint");
+ this->NSVE<rnumber>::write_checkpoint();
+ this->particles_output_writer_mpi->open_file(this->fs->get_current_fname());
+ this->particles_output_writer_mpi->template save<6>(
+ this->ps->getParticlesState(),
+ this->ps->getParticlesRhs(),
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->fs->iteration);
+ this->particles_output_writer_mpi->close_file();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_particles<rnumber>::finalize(void)
+{
+ TIMEZONE("NSVE_Stokes_particles::finalize");
+ delete this->pressure;
+ delete this->ps.release();
+ delete this->particles_output_writer_mpi;
+ delete this->particles_sample_writer_mpi;
+ this->NSVE<rnumber>::finalize();
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute fluid stats and sample fields at particle locations.
+ */
+
+template <typename rnumber>
+int NSVE_Stokes_particles<rnumber>::do_stats()
+{
+ TIMEZONE("NSVE_Stokes_particles::do_stats");
+ /// fluid stats go here
+ this->NSVE<rnumber>::do_stats();
+
+
+ /// either one of two conditions suffices to compute statistics:
+ /// 1) current iteration is a multiple of niter_part
+ /// 2) we are within niter_part_fine_duration/2 of a multiple of niter_part_fine_period
+ if (!(this->iteration % this->niter_part == 0 ||
+ ((this->iteration + this->niter_part_fine_duration/2) % this->niter_part_fine_period <=
+ this->niter_part_fine_duration)))
+ return EXIT_SUCCESS;
+
+ /// allocate temporary data array
+ /// initialize pdata0 with the positions, and pdata1 with the orientations
+ std::unique_ptr<double[]> pdata0 = this->ps->extractParticlesState(0, 3);
+ std::unique_ptr<double[]> pdata1 = this->ps->extractParticlesState(3, 6);
+ std::unique_ptr<double[]> pdata2(new double[9*this->ps->getLocalNbParticles()]);
+
+ /// sample position
+ this->particles_sample_writer_mpi->template save_dataset<3>(
+ "tracers0",
+ "position",
+ pdata0.get(), // we need to use pdata0.get() here, because it's 3D, whereas getParticlesState may give something else
+ &pdata0,
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->ps->get_step_idx()-1);
+
+ /// sample particle momentum
+ this->particles_sample_writer_mpi->template save_dataset<3>(
+ "tracers0",
+ "momentum",
+ pdata0.get(),
+ &pdata1,
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->ps->get_step_idx()-1);
+
+ /// sample velocity
+ std::fill_n(pdata1.get(), 3*this->ps->getLocalNbParticles(), 0);
+ if (!(this->iteration % this->niter_stat == 0))
+ {
+ // we need to compute velocity field manually, because it didn't happen in NSVE::do_stats()
+ this->fs->compute_velocity(this->fs->cvorticity);
+ *this->tmp_vec_field = this->fs->cvelocity->get_cdata();
+ this->tmp_vec_field->ift();
+ }
+ this->ps->sample_compute_field(*this->tmp_vec_field, pdata1.get());
+ this->particles_sample_writer_mpi->template save_dataset<3>(
+ "tracers0",
+ "velocity",
+ pdata0.get(),
+ &pdata1,
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->ps->get_step_idx()-1);
+
+ // deallocate temporary data array
+ delete[] pdata0.release();
+ delete[] pdata1.release();
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int NSVE_Stokes_particles<rnumber>::read_parameters(void)
+{
+ TIMEZONE("NSVE_Stokes_particles::read_parameters");
+ this->NSVE<rnumber>::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ this->niter_part = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part");
+ this->niter_part_fine_period = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_period");
+ this->niter_part_fine_duration = hdf5_tools::read_value<int>(parameter_file, "parameters/niter_part_fine_duration");
+ this->nparticles = hdf5_tools::read_value<int>(parameter_file, "parameters/nparticles");
+ this->tracers0_integration_steps = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_integration_steps");
+ this->tracers0_neighbours = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_neighbours");
+ this->tracers0_smoothness = hdf5_tools::read_value<int>(parameter_file, "parameters/tracers0_smoothness");
+ this->tracers0_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/tracers0_cutoff");
+ this->drag_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/drag_coefficient");
+ H5Fclose(parameter_file);
+ return EXIT_SUCCESS;
+}
+
+template class NSVE_Stokes_particles<float>;
+template class NSVE_Stokes_particles<double>;
+
diff --git a/cpp/full_code/NSVE_Stokes_particles.hpp b/cpp/full_code/NSVE_Stokes_particles.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e818041b966fe26a6163609e5f2fb50ee8f2489
--- /dev/null
+++ b/cpp/full_code/NSVE_Stokes_particles.hpp
@@ -0,0 +1,89 @@
+/**********************************************************************
+* *
+* 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 NSVE_STOKES_PARTICLES_HPP
+#define NSVE_STOKES_PARTICLES_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+
+/** \brief Navier-Stokes solver that includes simple Lagrangian tracers.
+ *
+ * Child of Navier Stokes vorticity equation solver, this class calls all the
+ * methods from `NSVE`, and in addition integrates simple Lagrangian tracers
+ * in the resulting velocity field.
+ */
+
+template <typename rnumber>
+class NSVE_Stokes_particles: public NSVE<rnumber>
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ int niter_part;
+ int niter_part_fine_period;
+ int niter_part_fine_duration;
+ int nparticles;
+ int tracers0_integration_steps;
+ int tracers0_neighbours;
+ int tracers0_smoothness;
+ double tracers0_cutoff;
+ double drag_coefficient;
+
+ /* other stuff */
+ std::unique_ptr<abstract_particles_system<long long int, double>> ps;
+ field<rnumber, FFTW, ONE> *pressure;
+
+ particles_output_hdf5<long long int, double,6> *particles_output_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+
+
+ NSVE_Stokes_particles(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ ~NSVE_Stokes_particles(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+};
+
+#endif//NSVE_STOKES_PARTICLES_HPP
+
diff --git a/cpp/full_code/NSVE_field_stats.cpp b/cpp/full_code/NSVE_field_stats.cpp
index 0969175cc75530e2dad2c3c5dd9e6a0449416ed0..7a4f1312ff15e291321058dbcc0907b1d681f374 100644
--- a/cpp/full_code/NSVE_field_stats.cpp
+++ b/cpp/full_code/NSVE_field_stats.cpp
@@ -101,6 +101,33 @@ int NSVE_field_stats<rnumber>::read_current_cvorticity(void)
}
return EXIT_SUCCESS;
}
+template <typename rnumber>
+int NSVE_field_stats<rnumber>::read_arbitrary_cvorticity(int arbitrary_iteration)
+{
+ TIMEZONE("NSVE_field_stats::read_arbitrary_cvorticity");
+ this->vorticity->real_space_representation = false;
+ if (this->bin_IO != NULL)
+ {
+ char itername[16];
+ sprintf(itername, "i%.5x", arbitrary_iteration);
+ std::string native_binary_fname = (
+ this->simname +
+ std::string("_cvorticity_") +
+ std::string(itername));
+ this->bin_IO->read(
+ native_binary_fname,
+ this->vorticity->get_cdata());
+ }
+ else
+ {
+ this->vorticity->io(
+ this->simname + std::string("_fields.h5"),
+ "vorticity",
+ arbitrary_iteration,
+ true);
+ }
+ return EXIT_SUCCESS;
+}
template <typename rnumber>
int NSVE_field_stats<rnumber>::finalize(void)
diff --git a/cpp/full_code/NSVE_field_stats.hpp b/cpp/full_code/NSVE_field_stats.hpp
index 28a2376f17ac2ac837cbacac828cd91572bb3a17..0939ba03b60e03f66073d1713d5c3905d9110dd6 100644
--- a/cpp/full_code/NSVE_field_stats.hpp
+++ b/cpp/full_code/NSVE_field_stats.hpp
@@ -59,6 +59,7 @@ class NSVE_field_stats: public postprocess
virtual int finalize(void);
int read_current_cvorticity(void);
+ int read_arbitrary_cvorticity(int iter);
};
#endif//NSVE_FIELD_STATS_HPP
diff --git a/cpp/full_code/NSVEp_extra_sampling.cpp b/cpp/full_code/NSVEp_extra_sampling.cpp
index 7b3e5a76c6d47c990df9698ccb5f8ef22770a70d..5a8cef55edf62575f950208f22ca028de1649932 100644
--- a/cpp/full_code/NSVEp_extra_sampling.cpp
+++ b/cpp/full_code/NSVEp_extra_sampling.cpp
@@ -149,6 +149,14 @@ int NSVEp_extra_sampling<rnumber>::do_stats()
return EXIT_SUCCESS;
}
+template <typename rnumber>
+int NSVEp_extra_sampling<rnumber>::read_parameters(void)
+{
+ TIMEZONE("NSVEp_extra_sampling::read_parameters");
+ this->NSVEparticles<rnumber>::read_parameters();
+ return EXIT_SUCCESS;
+}
+
template class NSVEp_extra_sampling<float>;
template class NSVEp_extra_sampling<double>;
diff --git a/cpp/full_code/bandpass_stats.cpp b/cpp/full_code/bandpass_stats.cpp
index 29d261e7ac62d3d791a7b6ed85827305df40b32d..26381daa9ab0d6243cea46fb3833c18bdab2164d 100644
--- a/cpp/full_code/bandpass_stats.cpp
+++ b/cpp/full_code/bandpass_stats.cpp
@@ -146,7 +146,7 @@ int bandpass_stats<rnumber>::work_on_current_iteration(void)
stat_group = 0;
/// loop over all bands
- for (unsigned int band_counter = 0; band_counter < this->k0list.size(); band_counter++)
+ for (int band_counter = 0; band_counter < int(this->k0list.size()); band_counter++)
{
*filtered_vel0 = *vel;
*filtered_vel1 = *vel;
diff --git a/cpp/full_code/field_single_to_double.cpp b/cpp/full_code/field_single_to_double.cpp
index 93a03aed5a494138ba8279792788a7bb19105325..bb234e8b8e227c7c9e322cec29d457b218bc6655 100644
--- a/cpp/full_code/field_single_to_double.cpp
+++ b/cpp/full_code/field_single_to_double.cpp
@@ -74,7 +74,9 @@ template <typename rnumber>
int field_single_to_double<rnumber>::work_on_current_iteration(void)
{
TIMEZONE("field_single_to_double::work_on_current_iteration");
+ DEBUG_MSG("before read_vorticity at iteration %d\n", this->iteration);
this->read_current_cvorticity();
+ DEBUG_MSG("after read_vorticity at iteration %d\n", this->iteration);
// using CLOOP as opposed to a global std::copy because CLOOP
// is openmp parallelized.
diff --git a/cpp/full_code/get_rfields.cpp b/cpp/full_code/get_rfields.cpp
index 45f6b5dce95b5d4fbb9edc2ce353fdde51f0fba8..c3cf6c9b1d7a802f831c4faac939ed07d5552f3c 100644
--- a/cpp/full_code/get_rfields.cpp
+++ b/cpp/full_code/get_rfields.cpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -35,8 +35,27 @@ int get_rfields<rnumber>::initialize(void)
TIMEZONE("get_rfields::initialize");
this->NSVE_field_stats<rnumber>::initialize();
DEBUG_MSG("after NSVE_field_stats::initialize\n");
+
+ // allocate kspace
this->kk = new kspace<FFTW, SMOOTH>(
this->vorticity->clayout, this->dkx, this->dky, this->dkz);
+ // allocate field for TrS2
+ this->traceS2 = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ this->vorticity->fftw_plan_rigor);
+ this->traceS2->real_space_representation = true;
+ // allocate field for velocity
+ this->vel = new field<rnumber, FFTW, THREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ this->vorticity->fftw_plan_rigor);
+ this->vel->real_space_representation = false;
+ // allocate field for velocity gradient
+ this->grad_vel = new field<rnumber, FFTW, THREExTHREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ this->vorticity->fftw_plan_rigor);
hid_t parameter_file = H5Fopen(
(this->simname + std::string(".h5")).c_str(),
H5F_ACC_RDONLY,
@@ -53,6 +72,9 @@ int get_rfields<rnumber>::initialize(void)
this->iteration_list = hdf5_tools::read_vector<int>(
parameter_file,
"/get_rfields/parameters/iteration_list");
+ this->TrS2_on = hdf5_tools::read_value<int>(
+ parameter_file,
+ "/get_rfields/parameters/TrS2_on");
H5Fclose(parameter_file);
return EXIT_SUCCESS;
}
@@ -62,53 +84,66 @@ int get_rfields<rnumber>::work_on_current_iteration(void)
{
TIMEZONE("get_rfields::work_on_current_iteration");
this->read_current_cvorticity();
- field<rnumber, FFTW, THREE> *vel = new field<rnumber, FFTW, THREE>(
- this->nx, this->ny, this->nz,
- this->comm,
- this->vorticity->fftw_plan_rigor);
- vel->real_space_representation = false;
- this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
- if (k2 <= this->kk->kM2 && k2 > 0)
- {
- vel->cval(cindex,0,0) = -(this->kk->ky[yindex]*this->vorticity->cval(cindex,2,1) - this->kk->kz[zindex]*this->vorticity->cval(cindex,1,1)) / k2;
- vel->cval(cindex,0,1) = (this->kk->ky[yindex]*this->vorticity->cval(cindex,2,0) - this->kk->kz[zindex]*this->vorticity->cval(cindex,1,0)) / k2;
- vel->cval(cindex,1,0) = -(this->kk->kz[zindex]*this->vorticity->cval(cindex,0,1) - this->kk->kx[xindex]*this->vorticity->cval(cindex,2,1)) / k2;
- vel->cval(cindex,1,1) = (this->kk->kz[zindex]*this->vorticity->cval(cindex,0,0) - this->kk->kx[xindex]*this->vorticity->cval(cindex,2,0)) / k2;
- vel->cval(cindex,2,0) = -(this->kk->kx[xindex]*this->vorticity->cval(cindex,1,1) - this->kk->ky[yindex]*this->vorticity->cval(cindex,0,1)) / k2;
- vel->cval(cindex,2,1) = (this->kk->kx[xindex]*this->vorticity->cval(cindex,1,0) - this->kk->ky[yindex]*this->vorticity->cval(cindex,0,0)) / k2;
- }
- else
- std::fill_n((rnumber*)(vel->get_cdata()+3*cindex), 6, 0.0);
- }
- );
- vel->symmetrize();
- vel->ift();
+ invert_curl(
+ this->kk,
+ this->vorticity,
+ this->vel);
+
+ /// compute velocity gradient
+ compute_gradient<rnumber, FFTW, THREE, THREExTHREE>(
+ this->kk,
+ this->vel,
+ this->grad_vel);
+
+ // compute TrS2
+ this->grad_vel->ift();
+ this->traceS2->RLOOP(
+ [&](ptrdiff_t rindex,
+ ptrdiff_t xindex,
+ ptrdiff_t yindex,
+ ptrdiff_t zindex){
+ rnumber AxxAxx = this->grad_vel->rval(rindex, 0, 0)*this->grad_vel->rval(rindex, 0, 0);
+ rnumber AyyAyy = this->grad_vel->rval(rindex, 1, 1)*this->grad_vel->rval(rindex, 1, 1);
+ rnumber AzzAzz = this->grad_vel->rval(rindex, 2, 2)*this->grad_vel->rval(rindex, 2, 2);
+ rnumber Sxy = this->grad_vel->rval(rindex, 0, 1) + this->grad_vel->rval(rindex, 1, 0);
+ rnumber Syz = this->grad_vel->rval(rindex, 1, 2) + this->grad_vel->rval(rindex, 2, 1);
+ rnumber Szx = this->grad_vel->rval(rindex, 2, 0) + this->grad_vel->rval(rindex, 0, 2);
+ this->traceS2->rval(rindex) = (
+ AxxAxx + AyyAyy + AzzAzz +
+ (Sxy*Sxy + Syz*Syz + Szx*Szx)/2);
+ });
std::string fname = (
this->simname +
std::string("_checkpoint_") +
std::to_string(this->iteration / (this->niter_out*this->checkpoints_per_file)) +
std::string(".h5"));
- vel->io(
+
+ /// output velocity field
+ this->vel->ift();
+ this->vel->io(
fname,
"velocity",
this->iteration,
false);
- delete vel;
-
+ /// output vorticity field
this->vorticity->ift();
this->vorticity->io(
fname,
"vorticity",
this->iteration,
false);
+
+ if (this->TrS2_on)
+ {
+ this->traceS2->io(
+ fname,
+ "TrS2",
+ this->iteration,
+ false);
+ }
return EXIT_SUCCESS;
}
@@ -117,6 +152,9 @@ int get_rfields<rnumber>::finalize(void)
{
TIMEZONE("get_rfields::finalize");
delete this->kk;
+ delete this->traceS2;
+ delete this->vel;
+ delete this->grad_vel;
this->NSVE_field_stats<rnumber>::finalize();
return EXIT_SUCCESS;
}
diff --git a/cpp/full_code/get_rfields.hpp b/cpp/full_code/get_rfields.hpp
index ae669aa3012ce492835a9737ca443cdc846e00ba..2d266d390322fbf3c2ecb631c95e074aef298491 100644
--- a/cpp/full_code/get_rfields.hpp
+++ b/cpp/full_code/get_rfields.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -42,7 +42,11 @@ class get_rfields: public NSVE_field_stats<rnumber>
public:
int checkpoints_per_file;
int niter_out;
+ int TrS2_on;
kspace<FFTW, SMOOTH> *kk;
+ field<rnumber, FFTW, ONE> *traceS2;
+ field<rnumber, FFTW, THREE> *vel;
+ field<rnumber, FFTW, THREExTHREE> *grad_vel;
get_rfields(
const MPI_Comm COMMUNICATOR,
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
index 212ba9b8184019b0ab24b45336ac15f7e5cf66d5..3ec76efc03ce3c39ab2151cfaac27346e71ed950 100644
--- a/cpp/full_code/kraichnan_field.cpp
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -104,26 +104,17 @@ int kraichnan_field<rnumber>::initialize(void)
DEBUG_MSG("Coefficient: %g\n",
this->spectrum_coefficient);
+ this->generate_random_velocity();
// is this the first iteration?
if ((this->iteration == 0) and (this->output_velocity == 1))
{
// if yes, generate random field and save it
- this->generate_random_velocity();
this->velocity->io(
this->get_current_fname(),
"velocity",
this->iteration,
false);
}
- else
- {
- // if not, read the random field that exists in the checkpoint file
- this->velocity->io(
- this->get_current_fname(),
- "velocity",
- this->iteration,
- true);
- }
return EXIT_SUCCESS;
}
@@ -136,12 +127,19 @@ int kraichnan_field<rnumber>::generate_random_velocity(void)
make_gaussian_random_field(
this->kk,
this->velocity,
- this->iteration, // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration
+ this->iteration, // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration. See note below
this->spectrum_slope,
this->spectrum_k_cutoff,
this->spectrum_coefficient * 3./2.); // incompressibility projection factor
// this->velocity is now in Fourier space
// project divfree, requires field in Fourier space
+ // Note on the choice of random seed:
+ // If in the future the simulation will continue with a smaller number of total threads (number of processes times number of threads per process),
+ // then during that run some of the threads will be seeded with a seed that has already been used for a previous iteration.
+ // So some sequences of Fourier modes will be identical to sequences of Fourier modes that occured in the past.
+ // Also see implementation of "make_gaussian_random_field".
+ // One work-around would be to multiply "this->iteration" with 10 or so ---
+ // it is unlikely the simulation will be continued with less than 0.1 of the initial total number of threads.
DEBUG_MSG("L2Norm before: %g\n",
this->velocity->L2norm(this->kk));
this->kk->template project_divfree<rnumber>(this->velocity->get_cdata());
@@ -329,14 +327,23 @@ void kraichnan_field<rnumber>::update_checkpoint()
hsize_t fields_stored;
hid_t fid, group_id;
fid = H5Fopen(fname.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
- group_id = H5Gopen(fid, "velocity/real", H5P_DEFAULT);
- H5Gget_num_objs(
+ group_id = H5Gopen(fid, "tracers0/state", H5P_DEFAULT);
+ bool dset_exists;
+ if (group_id > 0)
+ {
+ H5Gget_num_objs(
group_id,
&fields_stored);
- bool dset_exists = H5Lexists(
+ dset_exists = H5Lexists(
group_id,
std::to_string(this->iteration).c_str(),
H5P_DEFAULT);
+ }
+ else
+ {
+ fields_stored = 0;
+ dset_exists = false;
+ }
H5Gclose(group_id);
H5Fclose(fid);
if ((int(fields_stored) >= this->checkpoints_per_file) &&
@@ -353,13 +360,13 @@ void kraichnan_field<rnumber>::update_checkpoint()
H5P_DEFAULT);
hid_t gg = H5Gcreate(
fid,
- "velocity",
+ "tracers0",
H5P_DEFAULT,
H5P_DEFAULT,
H5P_DEFAULT);
hid_t ggg = H5Gcreate(
gg,
- "real",
+ "state",
H5P_DEFAULT,
H5P_DEFAULT,
H5P_DEFAULT);
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.cpp b/cpp/full_code/ornstein_uhlenbeck_process.cpp
index b6af42d324a27cddedd56aefbf7c57632ce186c6..653c2540f01ee18661346b04945c6cade09450d2 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.cpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.cpp
@@ -4,7 +4,7 @@
#include <cassert>
#include "scope_timer.hpp"
#include <algorithm>
-#define NDEBUG
+#include <chrono>
template <class rnumber,field_backend be>
@@ -16,6 +16,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
double ou_kmin,
double ou_kmax,
double ou_energy_amplitude,
+ double ou_gamma_factor,
double DKX,
double DKY,
double DKZ,
@@ -28,7 +29,7 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
this->ou_field = new field<rnumber,be,THREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
*this->ou_field = 0.0;
- this->ou_field->dft();
+ this->ou_field->dft();
this->ou_field_vort = new field<rnumber,be,THREE>(
nx,ny,nz, MPI_COMM_WORLD, FFTW_PLAN_RIGOR);
@@ -44,21 +45,29 @@ ornstein_uhlenbeck_process<rnumber,be>::ornstein_uhlenbeck_process(
this->ou_kmin_squ = pow(ou_kmin,2);
this->ou_kmax_squ = pow(ou_kmax,2);
this->ou_energy_amplitude = ou_energy_amplitude;
+ this->ou_gamma_factor = ou_gamma_factor;
this->epsilon = pow((this->ou_energy_amplitude/this->kolmogorov_constant), 3./2.);
- assert(this->kk->kM2 >= this->ou_kmax_squ);
+ assert(this->kk->kM2 >= this->ou_kmax_squ);
gen.resize(omp_get_max_threads());
+ long now;
+
+ if (this->ou_field->clayout->myrank == 0){
+ now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
+ }
+ MPI_Bcast(&now,1,MPI_LONG,0,this->ou_field->comm);
+
for(int thread=0;thread<omp_get_max_threads();thread++)
{
- int current_seed =
+ long current_seed =
this->ou_field->clayout->myrank*omp_get_max_threads() +
- thread;
+ thread+now;
gen[thread].seed(current_seed);
}
-
+
this->initialize_B();
}
@@ -179,10 +188,13 @@ void ornstein_uhlenbeck_process<rnumber,be>::initialize_B()
template <class rnumber, field_backend be>
void ornstein_uhlenbeck_process<rnumber, be>::let_converge(void)
{
- //add some logic here TODO
- for(int i=0; i<1000; i++)
+ double ou_kmin = sqrt(this->ou_kmin_squ);
+ double tau = 1.0/this->gamma(ou_kmin);
+ double dt = tau/1000.;
+
+ for(int i=0; i<2000; i++)
{
- this->step(2e-3);
+ this->step(dt);
}
}
diff --git a/cpp/full_code/ornstein_uhlenbeck_process.hpp b/cpp/full_code/ornstein_uhlenbeck_process.hpp
index a1aeaa98dd948ca39b2c96177341890b345c9953..5a1288c929afac4cad4a3a6f36ff5b2a0436b1ca 100644
--- a/cpp/full_code/ornstein_uhlenbeck_process.hpp
+++ b/cpp/full_code/ornstein_uhlenbeck_process.hpp
@@ -19,6 +19,7 @@ class ornstein_uhlenbeck_process{
double ou_kmin_squ;
double ou_kmax_squ;
double ou_energy_amplitude;
+ double ou_gamma_factor;
double kolmogorov_constant = 2;
double epsilon;
@@ -38,6 +39,7 @@ class ornstein_uhlenbeck_process{
double ou_kmin,
double ou_kmax,
double ou_energy_amplitude,
+ double ou_gamma_factor,
double DKX = 1.0,
double DKY = 1.0,
double DKZ = 1.0,
@@ -52,8 +54,7 @@ class ornstein_uhlenbeck_process{
inline double gamma(double kabs)
{
- return this->kolmogorov_constant*sqrt(this->kolmogorov_constant/
- this->ou_energy_amplitude);
+ return this->ou_gamma_factor*pow(kabs,2./3.)*sqrt(this->ou_energy_amplitude);
}
diff --git a/cpp/full_code/ou_vorticity_equation.cpp b/cpp/full_code/ou_vorticity_equation.cpp
index 37d3570d6f0c19849f5b584d4eaac5c7ed60c43f..561b70201e8167c08b3b4b14d457a20e5b583539 100644
--- a/cpp/full_code/ou_vorticity_equation.cpp
+++ b/cpp/full_code/ou_vorticity_equation.cpp
@@ -3,7 +3,6 @@
#include <cstring>
#include <cassert>
#include "scope_timer.hpp"
-#define NDEBUG
template <class rnumber, field_backend be>
diff --git a/cpp/full_code/ou_vorticity_equation.hpp b/cpp/full_code/ou_vorticity_equation.hpp
index e36292c2271616794815d43d084810b9daeb3fc9..8f8e110c64ca29bd1a5039a96f7cc8e632b5fc7b 100644
--- a/cpp/full_code/ou_vorticity_equation.hpp
+++ b/cpp/full_code/ou_vorticity_equation.hpp
@@ -24,6 +24,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
double ou_kmin,
double ou_kmax,
double ou_energy_amplitude,
+ double ou_gamma_factor,
double DKX = 1.0,
double DKY = 1.0,
double DKZ = 1.0,
@@ -34,7 +35,7 @@ class ou_vorticity_equation : public vorticity_equation<rnumber, be>
this->ou = new ornstein_uhlenbeck_process<rnumber,be>(
NAME,
nx, ny, nz,
- ou_kmin, ou_kmax, ou_energy_amplitude,
+ ou_kmin, ou_kmax, ou_energy_amplitude,ou_gamma_factor,
DKX, DKY, DKZ, FFTW_PLAN_RIGOR);
this->ou_forcing_type = "replace";
}
diff --git a/cpp/full_code/write_rpressure.cpp b/cpp/full_code/write_rpressure.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5835ee33625b6ac0240bb3c46e5aa9932b3a514a
--- /dev/null
+++ b/cpp/full_code/write_rpressure.cpp
@@ -0,0 +1,138 @@
+#include <string>
+#include <cmath>
+#include "write_rpressure.hpp"
+#include "scope_timer.hpp"
+
+
+template <typename rnumber>
+int write_rpressure<rnumber>::initialize(void)
+{
+
+ // initialize
+ this->NSVE_field_stats<rnumber>::initialize();
+
+ //iteration list needed by postprocess.cpp for the main loop
+ hid_t parameter_file = H5Fopen(
+ (this->simname + std::string("_post.h5")).c_str(),
+ H5F_ACC_RDONLY,
+ H5P_DEFAULT);
+ this->iteration_list = hdf5_tools::read_vector<int>(
+ parameter_file,
+ "/write_rpressure/parameters/iteration_list");
+ H5Fclose(parameter_file);
+ if (this->myrank==0) DEBUG_MSG("Iteration list[0] = %d \n", this->iteration_list[0]);
+
+ //get the pressure from here
+ this->ve = new vorticity_equation<rnumber, FFTW>(
+ this->simname.c_str(),
+ this->nx,
+ this->ny,
+ this->nz,
+ this->dkx,
+ this->dky,
+ this->dkz,
+ this->vorticity->fftw_plan_rigor);
+
+ // output field
+ this->pressure = new field<rnumber, FFTW, ONE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ this->vorticity->fftw_plan_rigor);
+
+ // write to binary
+ this->bin_IO_scalar = new field_binary_IO<rnumber, REAL, ONE>(
+ this->pressure->rlayout->sizes,
+ this->pressure->rlayout->subsizes,
+ this->pressure->rlayout->starts,
+ this->pressure->rlayout->comm);
+
+ if (this->myrank==0) DEBUG_MSG("Initialize end\n");
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int write_rpressure<rnumber>::clip_zero_padding(
+ field <rnumber, FFTW, ONE>* f)
+{
+ rnumber *a = f->get_rdata();
+ rnumber *b = f->get_rdata();
+ ptrdiff_t copy_size = f->rlayout->sizes[2] * ncomp(ONE);
+ ptrdiff_t skip_size = copy_size + 2*ncomp(ONE);
+ for (int i0 = 0; i0 < int(f->rlayout->subsizes[0]); i0++)
+ for (int i1 = 0; i1 < int(f->rlayout->sizes[1]); i1++)
+ {
+ std::copy(a, a + copy_size, b);
+ a += skip_size;
+ b += copy_size;
+ }
+ return EXIT_SUCCESS;
+}
+
+
+template <typename rnumber>
+int write_rpressure<rnumber>::work_on_current_iteration(void)
+{
+ if (this->myrank==0) DEBUG_MSG("\nComputing the pressure\n");
+ //compute_pressure
+ this->read_current_cvorticity();
+ *this->ve->cvorticity = this->vorticity->get_cdata();
+
+ this->ve->compute_velocity(this->ve->cvorticity);
+ this->ve->cvelocity->ift();
+ this->ve->compute_pressure(this->pressure);
+
+
+ if (this->myrank==0) DEBUG_MSG("Computed the pressure\n");
+
+ if (this->myrank==0) DEBUG_MSG("vorticity in ve fourier = [%lG, %lG; ...; %lG, %lG; %lG, %lG;.....; %lG, %lG;] \n", this->vorticity->get_cdata()[0][0], this->vorticity->get_cdata()[0][1], this->vorticity->get_cdata()[1][0], this->vorticity->get_cdata()[1][1], this->vorticity->get_cdata()[2][0], this->vorticity->get_cdata()[2][1],this->vorticity->get_cdata()[65][0], this->vorticity->get_cdata()[65][1]);
+
+ if (this->myrank==0) DEBUG_MSG("vorticity in ve fourier = [%lG, %lG; ...; %lG, %lG; %lG, %lG;.....; %lG, %lG;] \n", this->ve->cvorticity->get_cdata()[0][0], this->ve->cvorticity->get_cdata()[0][1], this->ve->cvorticity->get_cdata()[1][0], this->ve->cvorticity->get_cdata()[1][1], this->ve->cvorticity->get_cdata()[2][0], this->ve->cvorticity->get_cdata()[2][1],this->ve->cvorticity->get_cdata()[65][0], this->ve->cvorticity->get_cdata()[65][1]);
+ this->ve->compute_velocity(this->vorticity);
+
+ if (this->myrank==0) DEBUG_MSG("velocity in ve fourier = [%lG, %lG; ...; %lG, %lG; %lG, %lG;.....; %lG, %lG;] \n", this->ve->u->get_cdata()[0][0], this->ve->u->get_cdata()[0][1], this->ve->u->get_cdata()[1][0], this->ve->u->get_cdata()[1][1], this->ve->u->get_cdata()[2][0], this->ve->u->get_cdata()[2][1],this->ve->u->get_cdata()[65][0], this->ve->u->get_cdata()[65][1]);
+
+ if (this->myrank==0) DEBUG_MSG("pressure fourier = [%lG, %lG; ...; %lG, %lG; %lG, %lG;.....; %lG, %lG;] \n", this->pressure->get_cdata()[0][0], this->pressure->get_cdata()[0][1], this->pressure->get_cdata()[1][0], this->pressure->get_cdata()[1][1], this->pressure->get_cdata()[2][0], this->pressure->get_cdata()[2][1],this->pressure->get_cdata()[65][0], this->pressure->get_cdata()[65][1]);
+
+ //transform to real space
+ if (this->pressure->real_space_representation == false)
+ this->pressure->ift();
+
+ if (this->myrank==0) DEBUG_MSG("Real space\n");
+
+ if (this->myrank==0) DEBUG_MSG("pressure before clipping = [%lG; ...; %lG;.....; %lG;] \n", this->pressure->get_rdata()[0], this->pressure->get_rdata()[2], this->pressure->get_rdata()[13549]);
+
+ //write out to file
+ this->clip_zero_padding(pressure);
+ if (this->myrank==0) DEBUG_MSG("pressure after clipping = [%lG; ...; %lG;.....; %lG;] \n", this->pressure->get_rdata()[0], this->pressure->get_rdata()[2], this->pressure->get_rdata()[13549]);
+
+ //write to file -- THIS IS OK
+ char itername[16];
+ sprintf(itername, "i%.5x", this->iteration);
+ std::string native_binary_fname;
+
+ native_binary_fname = (
+ this->simname +
+ std::string("_rpressure_") +
+ std::string(itername));
+ if (this->myrank ==0) DEBUG_MSG("%s\n", native_binary_fname.c_str());
+
+ this->bin_IO_scalar->write(
+ native_binary_fname,
+ this->pressure->get_rdata());
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int write_rpressure<rnumber>::finalize(void)
+{
+
+ delete this->bin_IO_scalar;
+ delete this->pressure;
+ delete this->ve;
+ this->NSVE_field_stats<rnumber>::finalize();
+ return EXIT_SUCCESS;
+}
+
+template class write_rpressure<float>;
+template class write_rpressure<double>;
+
diff --git a/cpp/full_code/write_rpressure.hpp b/cpp/full_code/write_rpressure.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7e2d6d9c86c16bf3297fb7c5856ecfd86c23824
--- /dev/null
+++ b/cpp/full_code/write_rpressure.hpp
@@ -0,0 +1,68 @@
+/**********************************************************************
+* *
+* 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 WRITE_RPRESSURE
+#define WRITE_RPRESSURE
+
+#include <cstdlib>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <vector>
+#include "base.hpp"
+#include "field.hpp"
+#include "field_binary_IO.hpp"
+#include "vorticity_equation.hpp"
+#include "full_code/NSVE_field_stats.hpp"
+
+/** generate and write out pressure field in real space **/
+
+template <typename rnumber>
+class write_rpressure: public NSVE_field_stats<rnumber>
+{
+ public:
+ int niter_out;
+
+ vorticity_equation<rnumber, FFTW> *ve;
+ field<rnumber, FFTW, ONE> *pressure;
+ field_binary_IO <rnumber, REAL, ONE> *bin_IO_scalar;
+ int clip_zero_padding(field <rnumber, FFTW, ONE>* );
+
+ write_rpressure(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ NSVE_field_stats<rnumber>(
+ COMMUNICATOR,
+ simulation_name){}
+ virtual ~write_rpressure(){}
+
+ int initialize(void);
+ int work_on_current_iteration(void);
+ int finalize(void);
+};
+
+
+#endif//WRITE_PRESSURE
+
diff --git a/cpp/hdf5_tools.hpp b/cpp/hdf5_tools.hpp
index b217147ba3590d73ac7135458c5ce156ca8fc386..57a257625c5e50eef39438f85487962745893185 100644
--- a/cpp/hdf5_tools.hpp
+++ b/cpp/hdf5_tools.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/kspace.cpp b/cpp/kspace.cpp
index d601345742f6c0845e15b1b0cea2b2f1770c7a71..5f1d431550c3e9d28469a55e43dbb07bdf13dd69 100644
--- a/cpp/kspace.cpp
+++ b/cpp/kspace.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -118,15 +118,15 @@ kspace<be, dt>::kspace(
});
this->CLOOP_K2_NXMODES(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2,
- int nxmodes){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2,
+ const int nxmodes){
if (k2 < this->kM2)
{
- double knorm = sqrt(k2);
+ const double knorm = sqrt(k2);
kshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes*knorm;
nshell_local_thread.getMine()[int(knorm/this->dk)] += nxmodes;
}
@@ -230,11 +230,11 @@ void kspace<be, dt>::low_pass(
{
const double km2 = kmax*kmax;
this->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 >= km2)
std::fill_n((rnumber*)(a + ncomp(fc)*cindex), 2*ncomp(fc), 0);
});
@@ -264,15 +264,15 @@ void kspace<be, dt>::ball_filter(
{
const double prefactor0 = double(3) / pow(ell/2, 3);
this->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 > 0)
{
- double argument = sqrt(k2)*ell / 2;
- double prefactor = prefactor0 / pow(k2, 1.5);
+ const double argument = sqrt(k2)*ell / 2;
+ const double prefactor = prefactor0 / pow(k2, 1.5);
for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
prefactor *
@@ -300,18 +300,19 @@ void kspace<be, dt>::general_M_filter(
{
const double prefactor0 = 1.0;
this->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 > 0)
{
- double argument = sqrt(k2)*ell;
- double prefactor = prefactor0;
+ const double argument = sqrt(k2)*ell;
+ const double prefactor = prefactor0;
for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= (
- prefactor * (exp(-0.5*pow((2.9*argument),(68.0*(pow(ell,0.74))))) * sqrt(1.0 + (pow((argument/0.06),3.8))/(1.0 + (pow((argument/0.057),3.8))))));
+// prefactor * (exp(-0.5*pow((2.9*argument),(68.0*(pow(ell,0.74))))) * sqrt(1.0 + (pow((argument/0.06),3.8))/(1.0 + (pow((argument/0.057),3.8))))));
+ prefactor * (exp(-1.7*argument)) * sqrt(1.0 + 0.68*(pow((argument/0.065),4.6))/(1.0 + (pow((argument/0.065),4.6)))));
}
});
}
@@ -334,11 +335,11 @@ void kspace<be, dt>::Gauss_filter(
{
const double prefactor = - sigma*sigma/2;
this->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
{
for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= exp(prefactor*k2);
@@ -529,15 +530,15 @@ void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restri
this->low_pass<rnumber, fc>(a, this->kM);
break;
case SMOOTH:
- this->CLOOP(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
- double kk2 = (pow(this->kx[xindex]/this->kMx, 2) +
- pow(this->ky[yindex]/this->kMy, 2) +
- pow(this->kz[zindex]/this->kMz, 2));
- double tval = exp(-36.0 * (pow(kk2, 18)));
+ this->CLOOP_simd(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ const double kk2 = (pow(this->kx[xindex]/this->kMx, 2) +
+ pow(this->ky[yindex]/this->kMy, 2) +
+ pow(this->kz[zindex]/this->kMz, 2));
+ const double tval = exp(-36.0 * (pow(kk2, 18)));
for (unsigned int tcounter=0; tcounter<2*ncomp(fc); tcounter++)
((rnumber*)a)[2*ncomp(fc)*cindex + tcounter] *= tval;
});
@@ -554,11 +555,11 @@ void kspace<be, dt>::project_divfree(
{
TIMEZONE("kspace::project_divfree");
this->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 > 0)
{
typename fftw_interface<rnumber>::complex tval;
@@ -596,6 +597,7 @@ void kspace<be, dt>::project_divfree(
((*(a + cindex*3+1))[1])*((*(a + cindex*3+1))[1]) +
((*(a + cindex*3+2))[1])*((*(a + cindex*3+2))[1]));
if (projected_size > 0)
+ #pragma omp simd
for (int component=0; component<3; component++)
for (int imag_part=0; imag_part<2; imag_part++)
{
@@ -620,12 +622,12 @@ template <typename rnumber>
void kspace<be, dt>::rotate_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a)
{
TIMEZONE("kspace::rotate_divfree");
- this->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ this->CLOOP_K2_simd(
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 > 0)
{
double cosTheta;
@@ -684,7 +686,8 @@ void kspace<be, dt>::cospectrum(
const rnumber(* __restrict b)[2],
const hid_t group,
const std::string dset_name,
- const hsize_t toffset)
+ const hsize_t toffset,
+ const double wavenumber_exp)
{
TIMEZONE("field::cospectrum2");
shared_array<double> spec_local_thread(this->nshells*ncomp(fc)*ncomp(fc),[&](double* spec_local){
@@ -692,19 +695,20 @@ void kspace<be, dt>::cospectrum(
});
this->CLOOP_K2_NXMODES(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2,
- int nxmodes){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2,
+ const int nxmodes){
if (k2 <= this->kM2)
{
double* spec_local = spec_local_thread.getMine();
- int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
+ const int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
for (hsize_t i=0; i<ncomp(fc); i++)
for (hsize_t j=0; j<ncomp(fc); j++){
- spec_local[tmp_int + i*ncomp(fc)+j] += nxmodes * (
+ spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2,wavenumber_exp/2.)*
+ nxmodes * (
(a[ncomp(fc)*cindex + i][0] * b[ncomp(fc)*cindex + j][0]) +
(a[ncomp(fc)*cindex + i][1] * b[ncomp(fc)*cindex + j][1]));
}
@@ -762,7 +766,8 @@ void kspace<be, dt>::cospectrum(
const rnumber(* __restrict a)[2],
const hid_t group,
const std::string dset_name,
- const hsize_t toffset)
+ const hsize_t toffset,
+ const double wavenumber_exp)
{
TIMEZONE("field::cospectrum1");
shared_array<double> spec_local_thread(this->nshells*ncomp(fc)*ncomp(fc),[&](double* spec_local){
@@ -770,19 +775,21 @@ void kspace<be, dt>::cospectrum(
});
this->CLOOP_K2_NXMODES(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2,
- int nxmodes){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2,
+ const int nxmodes){
if (k2 <= this->kM2)
{
double* spec_local = spec_local_thread.getMine();
- int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
+ const int tmp_int = int(sqrt(k2) / this->dk)*ncomp(fc)*ncomp(fc);
+
for (hsize_t i=0; i<ncomp(fc); i++)
for (hsize_t j=0; j<ncomp(fc); j++){
- spec_local[tmp_int + i*ncomp(fc)+j] += nxmodes * (
+ spec_local[tmp_int + i*ncomp(fc)+j] += pow(k2, wavenumber_exp/2.)*
+ nxmodes * (
(a[ncomp(fc)*cindex + i][0] * a[ncomp(fc)*cindex + j][0]) +
(a[ncomp(fc)*cindex + i][1] * a[ncomp(fc)*cindex + j][1]));
}
@@ -845,12 +852,12 @@ double kspace<be, dt>::L2norm(
});
this->CLOOP_K2_NXMODES(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2,
- int nxmodes){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2,
+ const int nxmodes){
{
double* L2_local = L2_local_thread.getMine();
for (hsize_t i=0; i<ncomp(fc); i++){
@@ -1007,136 +1014,160 @@ template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const typename fftw_interface<float>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const typename fftw_interface<double>::complex *__restrict__ b,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, TWO_THIRDS>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<float, THREExTHREE>(
const typename fftw_interface<float>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, ONE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template void kspace<FFTW, SMOOTH>::cospectrum<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a,
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp);
template double kspace<FFTW, TWO_THIRDS>::L2norm<float, ONE>(
const typename fftw_interface<float>::complex *__restrict__ a);
diff --git a/cpp/kspace.hpp b/cpp/kspace.hpp
index 80d89f2593d7462b0dd8e28c229e3c7919a9e1ae..70342187672b8853dc5e7d54944d09497a1c51b9 100644
--- a/cpp/kspace.hpp
+++ b/cpp/kspace.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -25,7 +25,6 @@
#include <hdf5.h>
-#include <unordered_map>
#include <vector>
#include <string>
#include "omputils.hpp"
@@ -129,7 +128,8 @@ class kspace
const rnumber(* __restrict__ b)[2],
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp = 0);
template <typename rnumber,
field_components fc>
@@ -137,7 +137,8 @@ class kspace
const rnumber(* __restrict__ a)[2],
const hid_t group,
const std::string dset_name,
- const hsize_t toffset);
+ const hsize_t toffset,
+ const double wavenumber_exp = 0);
template <typename rnumber,
field_components fc>
@@ -154,12 +155,34 @@ class kspace
for (hsize_t yindex = 0; yindex < this->layout->subsizes[0]; yindex++){
for (hsize_t zindex = start; zindex < end; zindex++){
- ptrdiff_t cindex = yindex*this->layout->subsizes[1]*this->layout->subsizes[2]
- + zindex*this->layout->subsizes[2];
+ const ptrdiff_t cindex = (
+ yindex*this->layout->subsizes[1]*this->layout->subsizes[2] +
+ zindex*this->layout->subsizes[2]);
+ for (hsize_t xindex = 0; xindex < this->layout->subsizes[2]; xindex++)
+ {
+ expression(cindex + xindex, xindex, yindex, zindex);
+ }
+ }
+ }
+ }
+ }
+ template <class func_type>
+ void CLOOP_simd(func_type expression)
+ {
+ #pragma omp parallel
+ {
+ const hsize_t start = OmpUtils::ForIntervalStart(this->layout->subsizes[1]);
+ const hsize_t end = OmpUtils::ForIntervalEnd(this->layout->subsizes[1]);
+
+ for (hsize_t yindex = 0; yindex < this->layout->subsizes[0]; yindex++){
+ #pragma omp simd
+ for (hsize_t zindex = start; zindex < end; zindex++){
+ const ptrdiff_t cindex = (
+ yindex*this->layout->subsizes[1]*this->layout->subsizes[2] +
+ zindex*this->layout->subsizes[2]);
for (hsize_t xindex = 0; xindex < this->layout->subsizes[2]; xindex++)
{
- expression(cindex, xindex, yindex, zindex);
- cindex++;
+ expression(cindex + xindex, xindex, yindex, zindex);
}
}
}
@@ -175,15 +198,46 @@ class kspace
for (hsize_t yindex = 0; yindex < this->layout->subsizes[0]; yindex++){
for (hsize_t zindex = start; zindex < end; zindex++){
- ptrdiff_t cindex = yindex*this->layout->subsizes[1]*this->layout->subsizes[2]
+ const ptrdiff_t cindex = yindex*this->layout->subsizes[1]*this->layout->subsizes[2]
+ + zindex*this->layout->subsizes[2];
+ for (hsize_t xindex = 0; xindex < this->layout->subsizes[2]; xindex++)
+ {
+ expression(
+ cindex+xindex,
+ xindex,
+ yindex,
+ zindex,
+ (this->kx[xindex]*this->kx[xindex] +
+ this->ky[yindex]*this->ky[yindex] +
+ this->kz[zindex]*this->kz[zindex]));
+ }
+ }
+ }
+ }
+ }
+ template <class func_type>
+ void CLOOP_K2_simd(func_type expression)
+ {
+ #pragma omp parallel
+ {
+ const hsize_t start = OmpUtils::ForIntervalStart(this->layout->subsizes[1]);
+ const hsize_t end = OmpUtils::ForIntervalEnd(this->layout->subsizes[1]);
+
+ for (hsize_t yindex = 0; yindex < this->layout->subsizes[0]; yindex++){
+ #pragma omp simd
+ for (hsize_t zindex = start; zindex < end; zindex++){
+ const ptrdiff_t cindex = yindex*this->layout->subsizes[1]*this->layout->subsizes[2]
+ zindex*this->layout->subsizes[2];
for (hsize_t xindex = 0; xindex < this->layout->subsizes[2]; xindex++)
{
- double k2 = (this->kx[xindex]*this->kx[xindex] +
- this->ky[yindex]*this->ky[yindex] +
- this->kz[zindex]*this->kz[zindex]);
- expression(cindex, xindex, yindex, zindex, k2);
- cindex++;
+ expression(
+ cindex+xindex,
+ xindex,
+ yindex,
+ zindex,
+ (this->kx[xindex]*this->kx[xindex] +
+ this->ky[yindex]*this->ky[yindex] +
+ this->kz[zindex]*this->kz[zindex]));
}
}
}
@@ -199,21 +253,15 @@ class kspace
for (hsize_t yindex = 0; yindex < this->layout->subsizes[0]; yindex++){
for (hsize_t zindex = start; zindex < end; zindex++){
- ptrdiff_t cindex = yindex*this->layout->subsizes[1]*this->layout->subsizes[2]
+ const ptrdiff_t cindex = yindex*this->layout->subsizes[1]*this->layout->subsizes[2]
+ zindex*this->layout->subsizes[2];
- hsize_t xindex = 0;
- double k2 = (
+ const double k2 = (
this->ky[yindex]*this->ky[yindex] +
this->kz[zindex]*this->kz[zindex]);
- expression(cindex, xindex, yindex, zindex, k2, 1);
- cindex++;
- for (xindex = 1; xindex < this->layout->subsizes[2]; xindex++)
+ expression(cindex, 0, yindex, zindex, k2, 1);
+ for (hsize_t xindex = 1; xindex < this->layout->subsizes[2]; xindex++)
{
- k2 = (this->kx[xindex]*this->kx[xindex] +
- this->ky[yindex]*this->ky[yindex] +
- this->kz[zindex]*this->kz[zindex]);
- expression(cindex, xindex, yindex, zindex, k2, 2);
- cindex++;
+ expression(cindex+xindex, xindex, yindex, zindex, k2 + this->kx[xindex]*this->kx[xindex], 2);
}
}
}
diff --git a/cpp/particles/abstract_particles_input.hpp b/cpp/particles/abstract_particles_input.hpp
index 48c38bc592ddc442489d437327b421312bfd3f55..baac781c96413135875da8845e618eac7a3f7197 100644
--- a/cpp/particles/abstract_particles_input.hpp
+++ b/cpp/particles/abstract_particles_input.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/abstract_particles_output.hpp b/cpp/particles/abstract_particles_output.hpp
index 7a6a364c9ac21adc94d916f51fa74559df604d24..10ebbdd5c5dfae6185d318cbb699143e92bc6b7e 100644
--- a/cpp/particles/abstract_particles_output.hpp
+++ b/cpp/particles/abstract_particles_output.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -58,7 +58,7 @@ class abstract_particles_output {
std::unique_ptr<real_number[]> buffer_particles_positions_recv;
std::vector<std::unique_ptr<real_number[]>> buffer_particles_rhs_recv;
std::unique_ptr<partsize_t[]> buffer_indexes_recv;
- partsize_t size_buffers_recv;
+ partsize_t size_buffers_recv;
int buffers_size_particle_rhs_recv;
int nb_processes_involved;
diff --git a/cpp/particles/abstract_particles_system.hpp b/cpp/particles/abstract_particles_system.hpp
index 8a1be218346783f1ef9d20f8011610fd3527d28d..abcb8e684c5f269f10f47f337e36260486fba931 100644
--- a/cpp/particles/abstract_particles_system.hpp
+++ b/cpp/particles/abstract_particles_system.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -64,6 +64,7 @@ public:
virtual void shift_rhs_vectors() = 0;
virtual void completeLoop(const real_number dt) = 0;
+ virtual void complete2ndOrderLoop(const real_number dt) = 0;
virtual void completeLoopWithVorticity(
const real_number dt,
diff --git a/cpp/particles/abstract_particles_system_with_p2p.hpp b/cpp/particles/abstract_particles_system_with_p2p.hpp
index 3c31abf6f145ea1f130f64cae4eb281a1a83bdd2..8ce95b0c35cd36db153176b6f2ec6de6c0099839 100644
--- a/cpp/particles/abstract_particles_system_with_p2p.hpp
+++ b/cpp/particles/abstract_particles_system_with_p2p.hpp
@@ -1,3 +1,28 @@
+/******************************************************************************
+* *
+* Copyright 2019 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@ds.mpg.de *
+* *
+******************************************************************************/
+
+
+
#ifndef ABSTRACT_PARTICLES_SYSTEM_WITH_P2P_HPP
#define ABSTRACT_PARTICLES_SYSTEM_WITH_P2P_HPP
diff --git a/cpp/particles/alltoall_exchanger.hpp b/cpp/particles/alltoall_exchanger.hpp
index d3423523d9b9d02347514972c3bcb3f92129df56..288ff47e3b68e9fc5c759ce3547ec530b76e74d8 100644
--- a/cpp/particles/alltoall_exchanger.hpp
+++ b/cpp/particles/alltoall_exchanger.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/lock_free_bool_array.hpp b/cpp/particles/lock_free_bool_array.hpp
index 5e32a7d41bec3ddc7d56962d14d338a78f2b084a..4e5ad31b6274ee82b1dd99b0a8b87eba97fb64b7 100644
--- a/cpp/particles/lock_free_bool_array.hpp
+++ b/cpp/particles/lock_free_bool_array.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/p2p_computer.hpp b/cpp/particles/p2p_computer.hpp
index e5ddda47505d13b776f96191e925f580be3933af..47d193e095403d329cd704e8623567c0f080a39b 100644
--- a/cpp/particles/p2p_computer.hpp
+++ b/cpp/particles/p2p_computer.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/p2p_computer_empty.hpp b/cpp/particles/p2p_computer_empty.hpp
index cb9710ccf7e6255adeeee47fa9ea336e2caf8a75..ce4e282c3bdb7ff4bfc479b5a0b16d106f0492c1 100644
--- a/cpp/particles/p2p_computer_empty.hpp
+++ b/cpp/particles/p2p_computer_empty.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/p2p_cylinder_collisions.hpp b/cpp/particles/p2p_cylinder_collisions.hpp
index 0888349227d4ad7c60d8b3067dcfa1845d2dd3fe..44736ffc88643dd5664cdcd8440f19f0ea5778d2 100644
--- a/cpp/particles/p2p_cylinder_collisions.hpp
+++ b/cpp/particles/p2p_cylinder_collisions.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/p2p_distr_mpi.hpp b/cpp/particles/p2p_distr_mpi.hpp
index 8d2bb40919022c1a113bd519a6d3f8507d156df0..41f68f3ec8a9cab8ca7b1fabd6cf7c87752b0175 100644
--- a/cpp/particles/p2p_distr_mpi.hpp
+++ b/cpp/particles/p2p_distr_mpi.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/p2p_ghost_collisions.hpp b/cpp/particles/p2p_ghost_collisions.hpp
index fd3511c1eca2fa625dc7b43ec916cd495315dc11..82662695af316fd7bbce7e311feec2c9c1c48ee0 100644
--- a/cpp/particles/p2p_ghost_collisions.hpp
+++ b/cpp/particles/p2p_ghost_collisions.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/p2p_tree.hpp b/cpp/particles/p2p_tree.hpp
index cdb3089174ff888cbfc13810d18c617b4a8358e7..928957c02eb51fb02b9ad526ae05d61cd3052389 100644
--- a/cpp/particles/p2p_tree.hpp
+++ b/cpp/particles/p2p_tree.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_adams_bashforth.hpp b/cpp/particles/particles_adams_bashforth.hpp
index 21412e3530408a5980c376453cd6f5199466d830..f760d473890d8210a5b9927e524bad0a6d1df761 100644
--- a/cpp/particles/particles_adams_bashforth.hpp
+++ b/cpp/particles/particles_adams_bashforth.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -46,7 +46,7 @@ public:
TIMEZONE("particles_adams_bashforth::move_particles");
if(Max_steps < nb_rhs){
- throw std::runtime_error("Error, in bfps particles_adams_bashforth.\n"
+ throw std::runtime_error("Error, in TurTLE particles_adams_bashforth.\n"
"Step in particles_adams_bashforth is too large,"
"you must add formulation up this number or limit the number of steps.");
}
diff --git a/cpp/particles/particles_distr_mpi.hpp b/cpp/particles/particles_distr_mpi.hpp
index 57d8067059a1408afb90c53e32c83c0731030bd2..2cc73cd5d732ec3fd5877ae5d5a5fbc694f2f913 100644
--- a/cpp/particles/particles_distr_mpi.hpp
+++ b/cpp/particles/particles_distr_mpi.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -45,7 +45,7 @@ protected:
static const int MaxNbRhs = 10;
enum MpiTag{
- TAG_LOW_UP_NB_PARTICLES,
+ TAG_LOW_UP_NB_PARTICLES = 999,
TAG_UP_LOW_NB_PARTICLES,
TAG_LOW_UP_PARTICLES,
TAG_UP_LOW_PARTICLES,
@@ -65,6 +65,8 @@ protected:
TAG_UP_LOW_MOVED_PARTICLES_RHS = TAG_LOW_UP_MOVED_PARTICLES_RHS_MAX,
TAG_UP_LOW_MOVED_PARTICLES_RHS_MAX = TAG_UP_LOW_MOVED_PARTICLES_RHS+MaxNbRhs,
+
+ TAG_SHIFT_OFFSET
};
struct NeighborDescriptor{
@@ -114,6 +116,8 @@ protected:
std::vector<MPI_Request> mpiRequests;
std::vector<NeighborDescriptor> neigDescriptors;
+ int counter_shift_tags;
+
public:
////////////////////////////////////////////////////////////////////////////
@@ -124,7 +128,8 @@ public:
my_rank(-1), nb_processes(-1),nb_processes_involved(-1),
current_partition_interval(in_current_partitions),
current_partition_size(current_partition_interval.second-current_partition_interval.first),
- field_grid_dim(in_field_grid_dim){
+ field_grid_dim(in_field_grid_dim),
+ counter_shift_tags(0){
AssertMpi(MPI_Comm_rank(current_com, &my_rank));
AssertMpi(MPI_Comm_size(current_com, &nb_processes));
@@ -271,14 +276,14 @@ public:
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Isend(const_cast<partsize_t*>(&descriptor.nbParticlesToSend), 1, particles_utils::GetMpiType(partsize_t()),
- descriptor.destProc, TAG_LOW_UP_NB_PARTICLES,
+ descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_NB_PARTICLES,
current_com, &mpiRequests.back()));
if(descriptor.nbParticlesToSend){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(descriptor.nbParticlesToSend*size_particle_positions < std::numeric_limits<int>::max());
- AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_positions[0]), int(descriptor.nbParticlesToSend*size_particle_positions), particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_LOW_UP_PARTICLES,
+ AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_positions[0]), int(descriptor.nbParticlesToSend*size_particle_positions), particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_PARTICLES,
current_com, &mpiRequests.back()));
assert(descriptor.toRecvAndMerge == nullptr);
@@ -286,7 +291,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
mpiRequests.emplace_back();
assert(descriptor.nbParticlesToSend*size_particle_rhs < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(descriptor.toRecvAndMerge.get(), int(descriptor.nbParticlesToSend*size_particle_rhs), particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_UP_LOW_RESULTS,
+ AssertMpi(MPI_Irecv(descriptor.toRecvAndMerge.get(), int(descriptor.nbParticlesToSend*size_particle_rhs), particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_RESULTS,
current_com, &mpiRequests.back()));
}
}
@@ -295,7 +300,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{RECV_PARTICLES, idxDescr});
mpiRequests.emplace_back();
AssertMpi(MPI_Irecv(&descriptor.nbParticlesToRecv,
- 1, particles_utils::GetMpiType(partsize_t()), descriptor.destProc, TAG_UP_LOW_NB_PARTICLES,
+ 1, particles_utils::GetMpiType(partsize_t()), descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_NB_PARTICLES,
current_com, &mpiRequests.back()));
}
else{
@@ -303,7 +308,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Isend(const_cast<partsize_t*>(&descriptor.nbParticlesToSend), 1, particles_utils::GetMpiType(partsize_t()),
- descriptor.destProc, TAG_UP_LOW_NB_PARTICLES,
+ descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_NB_PARTICLES,
current_com, &mpiRequests.back()));
if(descriptor.nbParticlesToSend){
@@ -312,7 +317,7 @@ public:
assert(descriptor.nbParticlesToSend*size_particle_positions < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(const_cast<real_number*>(&particles_positions[(current_offset_particles_for_partition[current_partition_size-descriptor.nbPartitionsToSend])*size_particle_positions]),
int(descriptor.nbParticlesToSend*size_particle_positions), particles_utils::GetMpiType(real_number()),
- descriptor.destProc, TAG_UP_LOW_PARTICLES,
+ descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_PARTICLES,
current_com, &mpiRequests.back()));
assert(descriptor.toRecvAndMerge == nullptr);
@@ -320,7 +325,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{MERGE_PARTICLES, idxDescr});
mpiRequests.emplace_back();
assert(descriptor.nbParticlesToSend*size_particle_rhs < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(descriptor.toRecvAndMerge.get(), int(descriptor.nbParticlesToSend*size_particle_rhs), particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_LOW_UP_RESULTS,
+ AssertMpi(MPI_Irecv(descriptor.toRecvAndMerge.get(), int(descriptor.nbParticlesToSend*size_particle_rhs), particles_utils::GetMpiType(real_number()), descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_RESULTS,
current_com, &mpiRequests.back()));
}
@@ -328,15 +333,15 @@ public:
whatNext.emplace_back(std::pair<Action,int>{RECV_PARTICLES, idxDescr});
mpiRequests.emplace_back();
AssertMpi(MPI_Irecv(&descriptor.nbParticlesToRecv,
- 1, particles_utils::GetMpiType(partsize_t()), descriptor.destProc, TAG_LOW_UP_NB_PARTICLES,
+ 1, particles_utils::GetMpiType(partsize_t()), descriptor.destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_NB_PARTICLES,
current_com, &mpiRequests.back()));
}
}
}
const bool more_than_one_thread = (omp_get_max_threads() > 1);
- MPI_Barrier(MPI_COMM_WORLD);
- //DEBUG_MSG_WAIT(MPI_COMM_WORLD, "line 338 of particles_distr_mpi.hpp\n");
+ /// MPI_Barrier(MPI_COMM_WORLD);
+ /// DEBUG_MSG_WAIT(MPI_COMM_WORLD, "line 338 of particles_distr_mpi.hpp\n");
TIMEZONE_OMP_INIT_PREPARALLEL(omp_get_max_threads())
#pragma omp parallel default(shared)
@@ -376,7 +381,7 @@ public:
mpiRequests.emplace_back();
assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions),
- particles_utils::GetMpiType(real_number()), destProc, TAG_UP_LOW_PARTICLES,
+ particles_utils::GetMpiType(real_number()), destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_PARTICLES,
current_com, &mpiRequests.back()));
}
}
@@ -393,7 +398,7 @@ public:
mpiRequests.emplace_back();
assert(NbParticlesToReceive*size_particle_positions < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(descriptor.toCompute.get(), int(NbParticlesToReceive*size_particle_positions),
- particles_utils::GetMpiType(real_number()), destProc, TAG_LOW_UP_PARTICLES,
+ particles_utils::GetMpiType(real_number()), destProc, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_PARTICLES,
current_com, &mpiRequests.back()));
}
}
@@ -434,7 +439,7 @@ public:
const int destProc = descriptor.destProc;
whatNext.emplace_back(std::pair<Action,int>{RELEASE_BUFFER_PARTICLES, releasedAction.second});
mpiRequests.emplace_back();
- const int tag = descriptor.isLower? TAG_LOW_UP_RESULTS : TAG_UP_LOW_RESULTS;
+ const int tag = descriptor.isLower? TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_RESULTS : TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_RESULTS;
assert(NbParticlesToReceive*size_particle_rhs < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(descriptor.results.get(), int(NbParticlesToReceive*size_particle_rhs), particles_utils::GetMpiType(real_number()), destProc, tag,
current_com, &mpiRequests.back()));
@@ -526,6 +531,8 @@ public:
assert(whatNext.size() == 0);
assert(mpiRequests.size() == 0);
+
+ counter_shift_tags += 1;
}
@@ -635,35 +642,35 @@ public:
whatNext.emplace_back(std::pair<Action,int>{RECV_MOVE_NB_LOW, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Irecv(&nbNewFromLow, 1, particles_utils::GetMpiType(partsize_t()),
- (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_NB_PARTICLES,
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_NB_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
eventsBeforeWaitall += 1;
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Isend(const_cast<partsize_t*>(&nbOutLower), 1, particles_utils::GetMpiType(partsize_t()),
- (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_NB_PARTICLES,
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_NB_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
if(nbOutLower){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(nbOutLower*size_particle_positions < std::numeric_limits<int>::max());
- AssertMpi(MPI_Isend(&(*inout_positions_particles)[0], int(nbOutLower*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES,
+ AssertMpi(MPI_Isend(&(*inout_positions_particles)[0], int(nbOutLower*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(nbOutLower*size_particle_index < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(&(*inout_index_particles)[0], int(nbOutLower*size_particle_index), particles_utils::GetMpiType(partsize_t()),
- (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(nbOutLower*size_particle_rhs < std::numeric_limits<int>::max());
- AssertMpi(MPI_Isend(&inout_rhs_particles[idx_rhs][0], int(nbOutLower*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_RHS+idx_rhs,
+ AssertMpi(MPI_Isend(&inout_rhs_particles[idx_rhs][0], int(nbOutLower*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_RHS+idx_rhs,
MPI_COMM_WORLD, &mpiRequests.back()));
}
}
@@ -671,14 +678,14 @@ public:
whatNext.emplace_back(std::pair<Action,int>{RECV_MOVE_NB_UP, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Irecv(&nbNewFromUp, 1, particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved,
- TAG_LOW_UP_MOVED_NB_PARTICLES,
+ TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_NB_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
eventsBeforeWaitall += 1;
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
AssertMpi(MPI_Isend(const_cast<partsize_t*>(&nbOutUpper), 1, particles_utils::GetMpiType(partsize_t()),
- (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_NB_PARTICLES,
+ (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_NB_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
if(nbOutUpper){
@@ -686,14 +693,14 @@ public:
mpiRequests.emplace_back();
assert(nbOutUpper*size_particle_positions < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(&(*inout_positions_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_positions],
- int(nbOutUpper*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES,
+ int(nbOutUpper*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(nbOutUpper*size_particle_index < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(&(*inout_index_particles)[(myTotalNbParticles-nbOutUpper)*size_particle_index], int(nbOutUpper*size_particle_index),
- particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
+ particles_utils::GetMpiType(partsize_t()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
@@ -701,7 +708,7 @@ public:
mpiRequests.emplace_back();
assert(nbOutUpper*size_particle_rhs < std::numeric_limits<int>::max());
AssertMpi(MPI_Isend(&inout_rhs_particles[idx_rhs][(myTotalNbParticles-nbOutUpper)*size_particle_rhs],
- int(nbOutUpper*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_RHS+idx_rhs,
+ int(nbOutUpper*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_RHS+idx_rhs,
MPI_COMM_WORLD, &mpiRequests.back()));
}
}
@@ -726,7 +733,7 @@ public:
mpiRequests.emplace_back();
assert(nbNewFromLow*size_particle_positions < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(&newParticlesLow[0], int(nbNewFromLow*size_particle_positions), particles_utils::GetMpiType(real_number()),
- (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES,
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
newParticlesLowIndexes.reset(new partsize_t[nbNewFromLow*size_particle_index]);
@@ -735,7 +742,7 @@ public:
assert(nbNewFromLow*size_particle_index < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(&newParticlesLowIndexes[0], int(nbNewFromLow*size_particle_index),
particles_utils::GetMpiType(partsize_t()),
- (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
+ (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
@@ -743,7 +750,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(nbNewFromLow*size_particle_rhs < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(&newParticlesLowRhs[idx_rhs][0], int(nbNewFromLow*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_UP_LOW_MOVED_PARTICLES_RHS+idx_rhs,
+ AssertMpi(MPI_Irecv(&newParticlesLowRhs[idx_rhs][0], int(nbNewFromLow*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank-1+nb_processes_involved)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_UP_LOW_MOVED_PARTICLES_RHS+idx_rhs,
MPI_COMM_WORLD, &mpiRequests.back()));
}
}
@@ -756,7 +763,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{RECV_MOVE_UP, -1});
mpiRequests.emplace_back();
assert(nbNewFromUp*size_particle_positions < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(&newParticlesUp[0], int(nbNewFromUp*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES,
+ AssertMpi(MPI_Irecv(&newParticlesUp[0], int(nbNewFromUp*size_particle_positions), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES,
MPI_COMM_WORLD, &mpiRequests.back()));
newParticlesUpIndexes.reset(new partsize_t[nbNewFromUp*size_particle_index]);
@@ -765,7 +772,7 @@ public:
assert(nbNewFromUp*size_particle_index < std::numeric_limits<int>::max());
AssertMpi(MPI_Irecv(&newParticlesUpIndexes[0], int(nbNewFromUp*size_particle_index),
particles_utils::GetMpiType(partsize_t()),
- (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
+ (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_INDEXES,
MPI_COMM_WORLD, &mpiRequests.back()));
for(int idx_rhs = 0 ; idx_rhs < in_nb_rhs ; ++idx_rhs){
@@ -773,7 +780,7 @@ public:
whatNext.emplace_back(std::pair<Action,int>{NOTHING_TODO, -1});
mpiRequests.emplace_back();
assert(nbNewFromUp*size_particle_rhs < std::numeric_limits<int>::max());
- AssertMpi(MPI_Irecv(&newParticlesUpRhs[idx_rhs][0], int(nbNewFromUp*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_LOW_UP_MOVED_PARTICLES_RHS+idx_rhs,
+ AssertMpi(MPI_Irecv(&newParticlesUpRhs[idx_rhs][0], int(nbNewFromUp*size_particle_rhs), particles_utils::GetMpiType(real_number()), (my_rank+1)%nb_processes_involved, TAG_SHIFT_OFFSET*counter_shift_tags + TAG_LOW_UP_MOVED_PARTICLES_RHS+idx_rhs,
MPI_COMM_WORLD, &mpiRequests.back()));
}
}
@@ -880,6 +887,8 @@ public:
(*nb_particles) = myTotalNbParticles;
assert(mpiRequests.size() == 0);
+
+ counter_shift_tags += 1;
}
};
diff --git a/cpp/particles/particles_field_computer.hpp b/cpp/particles/particles_field_computer.hpp
index 054adab854f09f4a5e3c138dd810def68f7ca3aa..146e18073e660ae330404ffb4305df023523f620 100644
--- a/cpp/particles/particles_field_computer.hpp
+++ b/cpp/particles/particles_field_computer.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_generic_interp.hpp b/cpp/particles/particles_generic_interp.hpp
index da48641ca543dd853c24d675c1fea8b96f9da449..e75ae9b3d92b34b170f13bde065eb24afa9dafe8 100644
--- a/cpp/particles/particles_generic_interp.hpp
+++ b/cpp/particles/particles_generic_interp.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -26,12 +26,15 @@
#ifndef PARTICLES_GENERIC_INTERP_HPP
#define PARTICLES_GENERIC_INTERP_HPP
-template <class real_number, int interp_neighbours, int mode>
+template <class real_number, int interp_neighbours, int smoothness>
class particles_generic_interp;
#include "Lagrange_polys.hpp"
#include "spline.hpp"
+/*****************************************************************************/
+/* Lagrange interpolation */
+
template <>
class particles_generic_interp<double, 1,0>{
public:
@@ -43,122 +46,137 @@ public:
};
template <>
-class particles_generic_interp<double, 1,1>{
+class particles_generic_interp<double, 2,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n1_m1(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n2(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 1,2>{
+class particles_generic_interp<double, 3,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n1_m2(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n3(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 2,0>{
+class particles_generic_interp<double, 4,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n2(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n4(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 2,1>{
+class particles_generic_interp<double, 5,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n2_m1(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n5(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 2,2>{
+class particles_generic_interp<double, 6,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n2_m2(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n6(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 3,0>{
+class particles_generic_interp<double, 7,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n3(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n7(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 3,1>{
+class particles_generic_interp<double, 8,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n3_m1(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n8(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 3,2>{
+class particles_generic_interp<double, 9,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n3_m2(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n9(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 4,0>{
+class particles_generic_interp<double, 10,0>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n4(in_derivative, in_part_val, poly_val);
+ beta_Lagrange_n10(in_derivative, in_part_val, poly_val);
}
};
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* spline C1 */
+
template <>
-class particles_generic_interp<double, 4,1>{
+class particles_generic_interp<double, 1,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n4_m1(in_derivative, in_part_val, poly_val);
+ beta_n1_m1(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 4,2>{
+class particles_generic_interp<double, 2,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n4_m2(in_derivative, in_part_val, poly_val);
+ beta_n2_m1(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 5,0>{
+class particles_generic_interp<double, 3,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n5(in_derivative, in_part_val, poly_val);
+ beta_n3_m1(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 4,1>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n4_m1(in_derivative, in_part_val, poly_val);
}
};
@@ -173,95 +191,127 @@ public:
};
template <>
-class particles_generic_interp<double, 5,2>{
+class particles_generic_interp<double, 6,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n5_m2(in_derivative, in_part_val, poly_val);
+ beta_n6_m1(in_derivative, in_part_val, poly_val);
}
};
+template <>
+class particles_generic_interp<double, 7,1>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n7_m1(in_derivative, in_part_val, poly_val);
+ }
+};
template <>
-class particles_generic_interp<double, 6,0>{
+class particles_generic_interp<double, 8,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n6(in_derivative, in_part_val, poly_val);
+ beta_n8_m1(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 6,1>{
+class particles_generic_interp<double, 9,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n6_m1(in_derivative, in_part_val, poly_val);
+ beta_n9_m1(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 6,2>{
+class particles_generic_interp<double, 10,1>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n6_m2(in_derivative, in_part_val, poly_val);
+ beta_n10_m1(in_derivative, in_part_val, poly_val);
}
};
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* spline C2 */
template <>
-class particles_generic_interp<double, 7,0>{
+class particles_generic_interp<double, 1,2>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n7(in_derivative, in_part_val, poly_val);
+ beta_n1_m2(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 7,1>{
+class particles_generic_interp<double, 2,2>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n7_m1(in_derivative, in_part_val, poly_val);
+ beta_n2_m2(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 7,2>{
+class particles_generic_interp<double, 3,2>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n7_m2(in_derivative, in_part_val, poly_val);
+ beta_n3_m2(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 4,2>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n4_m2(in_derivative, in_part_val, poly_val);
}
};
+template <>
+class particles_generic_interp<double, 5,2>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n5_m2(in_derivative, in_part_val, poly_val);
+ }
+};
template <>
-class particles_generic_interp<double, 8,0>{
+class particles_generic_interp<double, 6,2>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n8(in_derivative, in_part_val, poly_val);
+ beta_n6_m2(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 8,1>{
+class particles_generic_interp<double, 7,2>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n8_m1(in_derivative, in_part_val, poly_val);
+ beta_n7_m2(in_derivative, in_part_val, poly_val);
}
};
@@ -275,67 +325,241 @@ public:
}
};
+template <>
+class particles_generic_interp<double, 9,2>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n9_m2(in_derivative, in_part_val, poly_val);
+ }
+};
template <>
-class particles_generic_interp<double, 9, 0>{
+class particles_generic_interp<double, 10,2>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n9(in_derivative, in_part_val, poly_val);
+ beta_n10_m2(in_derivative, in_part_val, poly_val);
}
};
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* spline C3 */
+
template <>
-class particles_generic_interp<double, 9,1>{
+class particles_generic_interp<double, 1,3>{
+ // no such thing as C3 interpolant with just a 4-point kernel, so we just use C2
+ // this specialisation must exist, otherwise the template_double_for_if will fail
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n9_m1(in_derivative, in_part_val, poly_val);
+ beta_n1_m2(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 9,2>{
+class particles_generic_interp<double, 2,3>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n9_m2(in_derivative, in_part_val, poly_val);
+ beta_n2_m3(in_derivative, in_part_val, poly_val);
}
};
+template <>
+class particles_generic_interp<double, 3,3>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n3_m3(in_derivative, in_part_val, poly_val);
+ }
+};
template <>
-class particles_generic_interp<double, 10,0>{
+class particles_generic_interp<double, 4,3>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_Lagrange_n10(in_derivative, in_part_val, poly_val);
+ beta_n4_m3(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 10,1>{
+class particles_generic_interp<double, 5,3>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n10_m1(in_derivative, in_part_val, poly_val);
+ beta_n5_m3(in_derivative, in_part_val, poly_val);
}
};
template <>
-class particles_generic_interp<double, 10,2>{
+class particles_generic_interp<double, 6,3>{
public:
using real_number = double;
void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
- beta_n10_m2(in_derivative, in_part_val, poly_val);
+ beta_n6_m3(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 7,3>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n7_m3(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 8,3>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n8_m3(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 9,3>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n9_m3(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 10,3>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n10_m3(in_derivative, in_part_val, poly_val);
+ }
+};
+
+/*****************************************************************************/
+
+/*****************************************************************************/
+/* spline C4 */
+
+template <>
+class particles_generic_interp<double, 1,4>{
+ // no such thing as C4 interpolant with just a 4-point kernel, so we just use C2
+ // this specialisation must exist, otherwise the template_double_for_if will fail
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n1_m2(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 2,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n2_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 3,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n3_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 4,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n4_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 5,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n5_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 6,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n6_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 7,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n7_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 8,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n8_m4(in_derivative, in_part_val, poly_val);
}
};
+template <>
+class particles_generic_interp<double, 9,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n9_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+template <>
+class particles_generic_interp<double, 10,4>{
+public:
+ using real_number = double;
+
+ void compute_beta(const int in_derivative, const double in_part_val, double poly_val[]) const {
+ beta_n10_m4(in_derivative, in_part_val, poly_val);
+ }
+};
+
+/*****************************************************************************/
+
#endif//PARTICLES_INTERP_SPLINE_HPP
diff --git a/cpp/particles/particles_inner_computer.cpp b/cpp/particles/particles_inner_computer.cpp
index 3a841bee50f2849ef981cb5c585ee448570ae2ca..d1e346c8376105aefc307870224d418c304d6a0d 100644
--- a/cpp/particles/particles_inner_computer.cpp
+++ b/cpp/particles/particles_inner_computer.cpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_inner_computer.hpp b/cpp/particles/particles_inner_computer.hpp
index 7f30ad6829e5cfa0ac40bd59db7a9a09cbe8ac6f..fd9dd678c27c2f1f878a40cecd7b2a1df4a0c727 100644
--- a/cpp/particles/particles_inner_computer.hpp
+++ b/cpp/particles/particles_inner_computer.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_inner_computer_2nd_order.hpp b/cpp/particles/particles_inner_computer_2nd_order.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..20ffda07fe3cbd7d47786a290718774f938ff7fe
--- /dev/null
+++ b/cpp/particles/particles_inner_computer_2nd_order.hpp
@@ -0,0 +1,84 @@
+/******************************************************************************
+* *
+* Copyright 2019 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@ds.mpg.de *
+* *
+******************************************************************************/
+
+
+
+#ifndef PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
+#define PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
+
+#include <cstring>
+#include <cassert>
+
+template <class real_number, class partsize_t>
+class particles_inner_computer_2nd_order_Stokes{
+ double drag_coefficient;
+public:
+ template <int size_particle_state, int size_particle_rhs>
+ void compute_interaction(const partsize_t number_of_particles, real_number particle_state[], real_number particle_rhs[]) const{
+ }
+
+ template <int size_particle_state>
+ void enforce_unit_orientation(const partsize_t /*nb_particles*/, real_number /*pos_part*/[]) const{
+ }
+
+ template <int size_particle_state, int size_particle_rhs>
+ void add_Lagrange_multipliers(const partsize_t /*nb_particles*/, real_number /*pos_part*/[], real_number /*rhs_part*/[]) const{
+ }
+
+ template <int size_particle_state, int size_particle_rhs, int size_particle_rhs_extra>
+ void compute_interaction_with_extra(
+ const partsize_t number_of_particles,
+ real_number particle_state[],
+ real_number particle_rhs[],
+ const real_number sampled_velocity[]) const{
+ assert(size_particle_state == 6);
+ assert(size_particle_rhs == 6);
+ assert(size_particle_rhs_extra == 3);
+ #pragma omp parallel for
+ for(partsize_t idx_part = 0 ; idx_part < number_of_particles ; ++idx_part){
+ particle_rhs[idx_part*size_particle_rhs + IDXC_X] = particle_state[idx_part*size_particle_state + 3 + IDXC_X];
+ particle_rhs[idx_part*size_particle_rhs + IDXC_Y] = particle_state[idx_part*size_particle_state + 3 + IDXC_Y];
+ particle_rhs[idx_part*size_particle_rhs + IDXC_Z] = particle_state[idx_part*size_particle_state + 3 + IDXC_Z];
+ particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_X] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_X] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_X]);
+ particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Y] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Y] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Y]);
+ particle_rhs[idx_part*size_particle_rhs + 3 + IDXC_Z] = - this->drag_coefficient * (particle_state[idx_part*size_particle_state + 3 + IDXC_Z] - sampled_velocity[idx_part*size_particle_rhs_extra + IDXC_Z]);
+ }
+ }
+
+ constexpr static bool isEnable() {
+ return true;
+ }
+
+ void set_drag_coefficient(double mu)
+ {
+ this->drag_coefficient = mu;
+ }
+
+ double get_drag_coefficient()
+ {
+ return this->drag_coefficient;
+ }
+};
+
+#endif//PARTICLES_INNER_COMPUTER_2ND_ORDER_HPP
+
diff --git a/cpp/particles/particles_inner_computer_empty.hpp b/cpp/particles/particles_inner_computer_empty.hpp
index a90d3aa1b9f5ca5e9e2085173c6c55a25809b469..aeff11c4e1c87256f216b3c70eda6fe78ea39f65 100644
--- a/cpp/particles/particles_inner_computer_empty.hpp
+++ b/cpp/particles/particles_inner_computer_empty.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_input_hdf5.hpp b/cpp/particles/particles_input_hdf5.hpp
index 3f895be3613030fca0a0fce1a786bb6fc541fe9c..1bd715c129942a3cbf3a78465eeefba02b9f021e 100644
--- a/cpp/particles/particles_input_hdf5.hpp
+++ b/cpp/particles/particles_input_hdf5.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_output_hdf5.hpp b/cpp/particles/particles_output_hdf5.hpp
index 6978a45c351bdb0b80e844dbb4e4fd7832390cc4..367e11c9ea4f7d457067d58459b3514188044971 100644
--- a/cpp/particles/particles_output_hdf5.hpp
+++ b/cpp/particles/particles_output_hdf5.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -29,6 +29,7 @@
#include <memory>
#include <vector>
#include <hdf5.h>
+#include <sys/stat.h>
#include "abstract_particles_output.hpp"
#include "scope_timer.hpp"
@@ -140,10 +141,23 @@ public:
if(Parent::isInvolved()){
if (Parent::getMyRank() == 0)
{
- hid_t file_id = H5Fopen(
+ bool file_exists = false;
+ {
+ struct stat file_buffer;
+ file_exists = (stat(filename.c_str(), &file_buffer) == 0);
+ }
+ hid_t file_id;
+ if (file_exists)
+ file_id = H5Fopen(
filename.c_str(),
H5F_ACC_RDWR | H5F_ACC_DEBUG,
H5P_DEFAULT);
+ else
+ file_id = H5Fcreate(
+ filename.c_str(),
+ H5F_ACC_EXCL | H5F_ACC_DEBUG,
+ H5P_DEFAULT,
+ H5P_DEFAULT);
assert(file_id >= 0);
bool group_exists = H5Lexists(
file_id,
diff --git a/cpp/particles/particles_output_mpiio.hpp b/cpp/particles/particles_output_mpiio.hpp
index b1c17898c3c2941e0ed161e40113a0d13c99b524..b05579e2e3cf013689ef0b2e30b763d285952e86 100644
--- a/cpp/particles/particles_output_mpiio.hpp
+++ b/cpp/particles/particles_output_mpiio.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_output_sampling_hdf5.hpp b/cpp/particles/particles_output_sampling_hdf5.hpp
index b3a4c51f8f03772ed5017aabdaa351e5453662d5..04174b15f46730de1f42d6cee15e0d36b255b58e 100644
--- a/cpp/particles/particles_output_sampling_hdf5.hpp
+++ b/cpp/particles/particles_output_sampling_hdf5.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_sampling.hpp b/cpp/particles/particles_sampling.hpp
index 672c080aea1b59e250109f94fd5fee388e199755..1606b0228b334e15355644610e02410ec4253653 100644
--- a/cpp/particles/particles_sampling.hpp
+++ b/cpp/particles/particles_sampling.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/particles/particles_system.hpp b/cpp/particles/particles_system.hpp
index 1b14faa5f7dd79c533054453d96e575d45e7d3c3..c5230f41d7267c1a8aaefc5e629d9535f5eb4268 100644
--- a/cpp/particles/particles_system.hpp
+++ b/cpp/particles/particles_system.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -368,7 +368,7 @@ public:
void move(const real_number dt) final {
TIMEZONE("particles_system::move");
positions_updater.move_particles(my_particles_positions.get(), my_nb_particles,
- my_particles_rhs.data(), std::min(step_idx,int(my_particles_rhs.size())),
+ my_particles_rhs.data(), std::min(step_idx, int(my_particles_rhs.size())),
dt);
}
@@ -415,6 +415,21 @@ public:
shift_rhs_vectors();
}
+ void complete2ndOrderLoop(const real_number dt) final {
+ assert(size_particle_positions == 6);
+ assert(size_particle_rhs == 6);
+ std::unique_ptr<real_number[]> sampled_velocity(new real_number[getLocalNbParticles()*3]());
+ std::fill_n(sampled_velocity.get(), 3*getLocalNbParticles(), 0);
+ sample_compute_field(default_field, sampled_velocity.get());
+ this->computer_particules_inner.template compute_interaction_with_extra<size_particle_positions, size_particle_rhs, 3>(
+ my_nb_particles, my_particles_positions.get(), my_particles_rhs.front().get(),
+ sampled_velocity.get());
+ move(dt);
+ redistribute();
+ inc_step_idx();
+ shift_rhs_vectors();
+ }
+
void completeLoopWithVorticity(
const real_number dt,
const real_number particle_extra_field[]) final {
diff --git a/cpp/particles/particles_system_builder.hpp b/cpp/particles/particles_system_builder.hpp
index dfc20dc6cac10ffcf98faf682026ea2e4cd745ee..e7db94985e225c475830ae4a7f03347df42f2462 100644
--- a/cpp/particles/particles_system_builder.hpp
+++ b/cpp/particles/particles_system_builder.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -45,6 +45,13 @@
///
/// Double template "for"
///
+/// "Template_double_for_if" is used to systematically generate specialized
+/// templates for two ranges of template parameters.
+/// For the interpolation we have `n` and `m` that designate "number of
+/// neighbours" and "smoothness of interpolant".
+/// Calling Template_double_for_if is essentially equivalent to having a couple
+/// of nested "switch" statements, but without all the repetitive code lines.
+///
//////////////////////////////////////////////////////////////////////////////
namespace Template_double_for_if{
@@ -291,7 +298,7 @@ inline std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>
const int in_current_iteration){
return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system<partsize_t, particles_rnumber>>,
int, 1, 11, 1, // interpolation_size
- int, 0, 3, 1, // spline_mode
+ int, 0, 5, 1, // spline_mode
particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
p2p_computer_empty<particles_rnumber,partsize_t>,
particles_inner_computer_empty<particles_rnumber,partsize_t>,
@@ -321,7 +328,7 @@ inline std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_
const particles_rnumber cutoff){
return Template_double_for_if::evaluate<std::unique_ptr<abstract_particles_system_with_p2p<partsize_t, particles_rnumber, p2p_computer_class>>,
int, 1, 11, 1, // interpolation_size
- int, 0, 3, 1, // spline_mode
+ int, 0, 5, 1, // spline_mode
particles_system_build_container<partsize_t, field_rnumber,be,fc,particles_rnumber,
p2p_computer_class,
particles_inner_computer_class,
diff --git a/cpp/particles/particles_utils.hpp b/cpp/particles/particles_utils.hpp
index 3b215035916b66dd5f53e7eb01b51554c4d6f35d..2f7e149568ff304818c8ab70e4023c6258c26f16 100644
--- a/cpp/particles/particles_utils.hpp
+++ b/cpp/particles/particles_utils.hpp
@@ -2,20 +2,20 @@
* *
* Copyright 2019 Max Planck Institute for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/scope_timer.hpp b/cpp/scope_timer.hpp
index 890f522c415d7a102a0fff25c5292502cbcb459c..98c23c8019d7003f2c123cab3cdeaefad793f25c 100644
--- a/cpp/scope_timer.hpp
+++ b/cpp/scope_timer.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n1.cpp b/cpp/spline_n1.cpp
index 62b2928001438876cce92ff9523ee0a9d57ad469..f311f9e6b14761f20ec859b14e3bfe6de1b95ea6 100644
--- a/cpp/spline_n1.cpp
+++ b/cpp/spline_n1.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n1.hpp b/cpp/spline_n1.hpp
index 707e9fbcd3a1772c07ba59e6494260f0d8740a7d..f090e07a8a0ee074ea394171de3de9142aca87d8 100644
--- a/cpp/spline_n1.hpp
+++ b/cpp/spline_n1.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n10.cpp b/cpp/spline_n10.cpp
index 1712663933addeecba17955f8632d1a20e8c32c6..b40be0f557ceb972faaf57cb656e596bcc5973de 100644
--- a/cpp/spline_n10.cpp
+++ b/cpp/spline_n10.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n10.hpp b/cpp/spline_n10.hpp
index 4853f1975e8759fbcaf4c499df02a5b3bc934a91..c0fa6d0a3f8d64caa4ed734726811f3dd8213baa 100644
--- a/cpp/spline_n10.hpp
+++ b/cpp/spline_n10.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n2.cpp b/cpp/spline_n2.cpp
index b2fc17f4c6be63883bbc4346d43bb59f238d1bad..be015b15b45e9e1944d52fd8405ecd21cae93025 100644
--- a/cpp/spline_n2.cpp
+++ b/cpp/spline_n2.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n2.hpp b/cpp/spline_n2.hpp
index b29e7be13f663d9bdbca310b91b7e3a1db92a205..033d4a42429ddb947054512326c88bab97a9dfa4 100644
--- a/cpp/spline_n2.hpp
+++ b/cpp/spline_n2.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n3.cpp b/cpp/spline_n3.cpp
index ef83a28e31da17d204693b217affc493dd5e632a..b1ac72d768b537365392ee5a00e2f07606b6d48e 100644
--- a/cpp/spline_n3.cpp
+++ b/cpp/spline_n3.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n3.hpp b/cpp/spline_n3.hpp
index 159f8b970947b507205f40fcfb5140ca683770ce..57206746240b948eeaf7d32db8c04bbbbf6f4937 100644
--- a/cpp/spline_n3.hpp
+++ b/cpp/spline_n3.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n4.cpp b/cpp/spline_n4.cpp
index 1e42c7838d5bbc0c3fe9a78ebc45993112d72568..ccb73c4534333c76c9e6e29ed55a380d0a99737a 100644
--- a/cpp/spline_n4.cpp
+++ b/cpp/spline_n4.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n4.hpp b/cpp/spline_n4.hpp
index 4b80d16ea110fe1738f1eecca5f2f9cbe6e43e84..13417a3d4ae49d99608c47a45c9ae49bc4c53dac 100644
--- a/cpp/spline_n4.hpp
+++ b/cpp/spline_n4.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n5.cpp b/cpp/spline_n5.cpp
index 90e3a6e46424eca376ed85abf4f67a559e3cb910..dee1a574427a606777ba6d55ceb9c238cfcac9a1 100644
--- a/cpp/spline_n5.cpp
+++ b/cpp/spline_n5.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n5.hpp b/cpp/spline_n5.hpp
index fccd512076dc571fdded3fe3ccf44f7e4132ec48..fa4371acc2c9f41da1200f7d5096e4d84413f918 100644
--- a/cpp/spline_n5.hpp
+++ b/cpp/spline_n5.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n6.cpp b/cpp/spline_n6.cpp
index 4fe0101d171a53519ec4ae05cda2cce8f7da9b0b..3e9285470c60970f6404db0fb3fd40e517f2b374 100644
--- a/cpp/spline_n6.cpp
+++ b/cpp/spline_n6.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n6.hpp b/cpp/spline_n6.hpp
index 8eb99092eafd0f1f7c783a023eae49d714479ce1..2a0a0194892df1982089d0ab94fa67d37d1a647e 100644
--- a/cpp/spline_n6.hpp
+++ b/cpp/spline_n6.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n7.cpp b/cpp/spline_n7.cpp
index 84a3fdbc99a96cd3f5a13501414cbaa0a908096d..01407629782e07ba9b473da9226f51c594c1c9be 100644
--- a/cpp/spline_n7.cpp
+++ b/cpp/spline_n7.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n7.hpp b/cpp/spline_n7.hpp
index 2c0b86f6ff8bd722fd6bd038d8dfee39727d5fc1..cceab8ffe259534ae39799082e21b0a6bf56c249 100644
--- a/cpp/spline_n7.hpp
+++ b/cpp/spline_n7.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n8.cpp b/cpp/spline_n8.cpp
index d6cea769acbc3a74eee1dfd53933a45441318016..53841835a32c1f1ebb14a4a3a1cc11f6c70c8719 100644
--- a/cpp/spline_n8.cpp
+++ b/cpp/spline_n8.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n8.hpp b/cpp/spline_n8.hpp
index 39b2e03ad321b495a1da04f1fae8526daa4d3536..5377f5402a086778f8ad04571b19507da6d86e01 100644
--- a/cpp/spline_n8.hpp
+++ b/cpp/spline_n8.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n9.cpp b/cpp/spline_n9.cpp
index 4587a3f91c27937f553430463d19ea2221a173dc..ed92aeb14a0120001c39bbafebef2bb91437a3fc 100644
--- a/cpp/spline_n9.cpp
+++ b/cpp/spline_n9.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/spline_n9.hpp b/cpp/spline_n9.hpp
index b3770a824b2faec1f9196a2b366ed729e6ad21c7..35882791ceb3a486a5eff3d3f5e62ebf83a33287 100644
--- a/cpp/spline_n9.hpp
+++ b/cpp/spline_n9.hpp
@@ -3,20 +3,20 @@
* Copyright 2017 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
diff --git a/cpp/vorticity_equation.cpp b/cpp/vorticity_equation.cpp
index 58de8a5eca1d331b9fe4b4915603504c7aaab289..8518fd1e9c0d36310cb063b68e92a201903135f8 100644
--- a/cpp/vorticity_equation.cpp
+++ b/cpp/vorticity_equation.cpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -176,11 +176,11 @@ void vorticity_equation<rnumber, be>::compute_vorticity()
TIMEZONE("vorticity_equation::compute_vorticity");
this->cvorticity->real_space_representation = false;
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
this->cvorticity->cval(cindex,0,0) = -(this->kk->ky[yindex]*this->u->cval(cindex,2,1) - this->kk->kz[zindex]*this->u->cval(cindex,1,1));
@@ -204,24 +204,27 @@ void vorticity_equation<rnumber, be>::compute_velocity(field<rnumber, be, THREE>
TIMEZONE("vorticity_equation::compute_velocity");
this->u->real_space_representation = false;
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
- if (k2 <= this->kk->kM2 && k2 > 0)
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ const double kk2 = (k2 > 0) ? k2 : 1.;
+ if (k2 <= this->kk->kM2)
{
- this->u->cval(cindex,0,0) = -(this->kk->ky[yindex]*vorticity->cval(cindex,2,1) - this->kk->kz[zindex]*vorticity->cval(cindex,1,1)) / k2;
- this->u->cval(cindex,0,1) = (this->kk->ky[yindex]*vorticity->cval(cindex,2,0) - this->kk->kz[zindex]*vorticity->cval(cindex,1,0)) / k2;
- this->u->cval(cindex,1,0) = -(this->kk->kz[zindex]*vorticity->cval(cindex,0,1) - this->kk->kx[xindex]*vorticity->cval(cindex,2,1)) / k2;
- this->u->cval(cindex,1,1) = (this->kk->kz[zindex]*vorticity->cval(cindex,0,0) - this->kk->kx[xindex]*vorticity->cval(cindex,2,0)) / k2;
- this->u->cval(cindex,2,0) = -(this->kk->kx[xindex]*vorticity->cval(cindex,1,1) - this->kk->ky[yindex]*vorticity->cval(cindex,0,1)) / k2;
- this->u->cval(cindex,2,1) = (this->kk->kx[xindex]*vorticity->cval(cindex,1,0) - this->kk->ky[yindex]*vorticity->cval(cindex,0,0)) / k2;
+ this->u->cval(cindex,0,0) = -(this->kk->ky[yindex]*vorticity->cval(cindex,2,1) - this->kk->kz[zindex]*vorticity->cval(cindex,1,1)) / kk2;
+ this->u->cval(cindex,0,1) = (this->kk->ky[yindex]*vorticity->cval(cindex,2,0) - this->kk->kz[zindex]*vorticity->cval(cindex,1,0)) / kk2;
+ this->u->cval(cindex,1,0) = -(this->kk->kz[zindex]*vorticity->cval(cindex,0,1) - this->kk->kx[xindex]*vorticity->cval(cindex,2,1)) / kk2;
+ this->u->cval(cindex,1,1) = (this->kk->kz[zindex]*vorticity->cval(cindex,0,0) - this->kk->kx[xindex]*vorticity->cval(cindex,2,0)) / kk2;
+ this->u->cval(cindex,2,0) = -(this->kk->kx[xindex]*vorticity->cval(cindex,1,1) - this->kk->ky[yindex]*vorticity->cval(cindex,0,1)) / kk2;
+ this->u->cval(cindex,2,1) = (this->kk->kx[xindex]*vorticity->cval(cindex,1,0) - this->kk->ky[yindex]*vorticity->cval(cindex,0,0)) / kk2;
}
else
std::fill_n((rnumber*)(this->u->get_cdata()+3*cindex), 6, 0.0);
}
);
+ if (this->kk->layout->myrank == 0)
+ std::fill_n((rnumber*)this->u->get_cdata(), 6, 0.0);
this->u->symmetrize();
}
@@ -229,8 +232,8 @@ template <class rnumber,
field_backend be>
void vorticity_equation<rnumber, be>::add_Kolmogorov_forcing(
field<rnumber, be, THREE> *dst,
- int fmode,
- double famplitude)
+ const int fmode,
+ const double famplitude)
{
TIMEZONE("vorticity_equation::add_Kolmogorov_forcing");
ptrdiff_t cindex;
@@ -251,18 +254,18 @@ template <class rnumber,
void vorticity_equation<rnumber, be>::add_field_band(
field<rnumber, be, THREE> *dst,
field<rnumber, be, THREE> *src,
- double k0, double k1,
- double prefactor)
+ const double k0, const double k1,
+ const double prefactor)
{
TIMEZONE("vorticity_equation::add_field_band");
this->kk->CLOOP(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
- 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]);
+ [&](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++)
@@ -332,13 +335,13 @@ void vorticity_equation<rnumber, be>::add_forcing(
shared_array<double> local_energy_in_shell(1);
double energy_in_shell = 0;
this->kk->CLOOP_K2_NXMODES(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2,
- int nxmodes){
- double knorm = sqrt(k2);
+ [&](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))
@@ -396,15 +399,15 @@ void vorticity_equation<rnumber, be>::impose_forcing(
shared_array<double> local_total_energy(1);
double energy_in_shell, total_energy;
this->kk->CLOOP_K2_NXMODES(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2,
- int nxmodes){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2,
+ const int nxmodes){
if (k2 > 0)
{
- double mode_energy = nxmodes*(
+ const double mode_energy = nxmodes*(
onew->cval(cindex, 0, 0)*onew->cval(cindex, 0, 0) + onew->cval(cindex, 0, 1)*onew->cval(cindex, 0, 1) +
onew->cval(cindex, 1, 0)*onew->cval(cindex, 1, 0) + onew->cval(cindex, 1, 1)*onew->cval(cindex, 1, 1) +
onew->cval(cindex, 2, 0)*onew->cval(cindex, 2, 0) + onew->cval(cindex, 2, 1)*onew->cval(cindex, 2, 1)
@@ -439,12 +442,12 @@ void vorticity_equation<rnumber, be>::impose_forcing(
// see Michael's thesis, page 38
double temp_famplitude = sqrt((this->energy - total_energy + energy_in_shell) / energy_in_shell);
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
- double knorm = sqrt(k2);
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
+ const double knorm = sqrt(k2);
if ((this->fk0 <= knorm) &&
(this->fk1 >= knorm))
for (int c=0; c<3; c++)
@@ -461,7 +464,7 @@ template <class rnumber,
void vorticity_equation<rnumber, be>::omega_nonlin(
int src)
{
- DEBUG_MSG("vorticity_equation::omega_nonlin(%d)\n", src);
+ //DEBUG_MSG("vorticity_equation::omega_nonlin(%d)\n", src);
assert(src >= 0 && src < 3);
this->compute_velocity(this->v[src]);
/* get fields from Fourier space to real space */
@@ -471,10 +474,10 @@ void vorticity_equation<rnumber, be>::omega_nonlin(
this->rvorticity->ift();
/* compute cross product $u \times \omega$, and normalize */
this->u->RLOOP(
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
rnumber tmp[3];
for (int cc=0; cc<3; cc++)
tmp[cc] = (this->u->rval(rindex,(cc+1)%3)*this->rvorticity->rval(rindex,(cc+2)%3) -
@@ -489,19 +492,29 @@ void vorticity_equation<rnumber, be>::omega_nonlin(
this->kk->template dealias<rnumber, THREE>(this->u->get_cdata());
/* $\imath k \times Fourier(u \times \omega)$ */
this->kk->CLOOP(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
- rnumber tmp[3][2];
- {
- tmp[0][0] = -(this->kk->ky[yindex]*this->u->cval(cindex,2,1) - this->kk->kz[zindex]*this->u->cval(cindex,1,1));
- tmp[1][0] = -(this->kk->kz[zindex]*this->u->cval(cindex,0,1) - this->kk->kx[xindex]*this->u->cval(cindex,2,1));
- tmp[2][0] = -(this->kk->kx[xindex]*this->u->cval(cindex,1,1) - this->kk->ky[yindex]*this->u->cval(cindex,0,1));
- tmp[0][1] = (this->kk->ky[yindex]*this->u->cval(cindex,2,0) - this->kk->kz[zindex]*this->u->cval(cindex,1,0));
- tmp[1][1] = (this->kk->kz[zindex]*this->u->cval(cindex,0,0) - this->kk->kx[xindex]*this->u->cval(cindex,2,0));
- tmp[2][1] = (this->kk->kx[xindex]*this->u->cval(cindex,1,0) - this->kk->ky[yindex]*this->u->cval(cindex,0,0));
- }
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
+ // this tmp needs to be defined as const
+ const rnumber tmp[3][2] = {
+ {(rnumber)(-(this->kk->ky[yindex]*this->u->cval(cindex,2,1) - this->kk->kz[zindex]*this->u->cval(cindex,1,1))),
+ (rnumber)( (this->kk->ky[yindex]*this->u->cval(cindex,2,0) - this->kk->kz[zindex]*this->u->cval(cindex,1,0)))},
+ {(rnumber)(-(this->kk->kz[zindex]*this->u->cval(cindex,0,1) - this->kk->kx[xindex]*this->u->cval(cindex,2,1))),
+ (rnumber)( (this->kk->kz[zindex]*this->u->cval(cindex,0,0) - this->kk->kx[xindex]*this->u->cval(cindex,2,0)))},
+ {(rnumber)(-(this->kk->kx[xindex]*this->u->cval(cindex,1,1) - this->kk->ky[yindex]*this->u->cval(cindex,0,1))),
+ (rnumber)( (this->kk->kx[xindex]*this->u->cval(cindex,1,0) - this->kk->ky[yindex]*this->u->cval(cindex,0,0)))}
+ };
+ //rnumber tmp[3][2];
+ //{
+ // tmp[0][0] = -(this->kk->ky[yindex]*this->u->cval(cindex,2,1) - this->kk->kz[zindex]*this->u->cval(cindex,1,1));
+ // tmp[1][0] = -(this->kk->kz[zindex]*this->u->cval(cindex,0,1) - this->kk->kx[xindex]*this->u->cval(cindex,2,1));
+ // tmp[2][0] = -(this->kk->kx[xindex]*this->u->cval(cindex,1,1) - this->kk->ky[yindex]*this->u->cval(cindex,0,1));
+ // tmp[0][1] = (this->kk->ky[yindex]*this->u->cval(cindex,2,0) - this->kk->kz[zindex]*this->u->cval(cindex,1,0));
+ // tmp[1][1] = (this->kk->kz[zindex]*this->u->cval(cindex,0,0) - this->kk->kx[xindex]*this->u->cval(cindex,2,0));
+ // tmp[2][1] = (this->kk->kx[xindex]*this->u->cval(cindex,1,0) - this->kk->ky[yindex]*this->u->cval(cindex,0,0));
+ //}
+ //#pragma omp simd
for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
this->u->cval(cindex, cc, i) = tmp[cc][i];
}
@@ -515,20 +528,19 @@ template <class rnumber,
field_backend be>
void vorticity_equation<rnumber, be>::step(double dt)
{
- DEBUG_MSG("vorticity_equation::step\n");
+ //DEBUG_MSG("vorticity_equation::step\n");
TIMEZONE("vorticity_equation::step");
*this->v[1] = 0.0;
this->omega_nonlin(0);
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
- double factor0;
- factor0 = exp(-this->nu * k2 * dt);
+ const double factor0 = exp(-this->nu * k2 * dt);
for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
this->v[1]->cval(cindex,cc,i) = (
this->v[0]->cval(cindex,cc,i) +
@@ -540,16 +552,15 @@ void vorticity_equation<rnumber, be>::step(double dt)
this->omega_nonlin(1);
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
- double factor0, factor1;
- factor0 = exp(-this->nu * k2 * dt/2);
- factor1 = exp( this->nu * k2 * dt/2);
+ const double factor0 = exp(-this->nu * k2 * dt/2);
+ const double factor1 = exp( this->nu * k2 * dt/2);
for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
this->v[2]->cval(cindex, cc, i) = (
3*this->v[0]->cval(cindex,cc,i)*factor0 +
@@ -564,15 +575,14 @@ void vorticity_equation<rnumber, be>::step(double dt)
// store old vorticity
*this->v[1] = *this->v[0];
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
- double factor0;
- factor0 = exp(-this->nu * k2 * dt * 0.5);
+ const double factor0 = exp(-this->nu * k2 * dt * 0.5);
for (int cc=0; cc<3; cc++) for (int i=0; i<2; i++)
this->v[3]->cval(cindex,cc,i) = (
this->v[0]->cval(cindex,cc,i)*factor0 +
@@ -595,14 +605,16 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
TIMEZONE("vorticity_equation::compute_pressure");
pressure->real_space_representation = false;
/* assume velocity is already in real space representation */
-
+ /* set the pressure representation to Fourier space */
+ pressure->real_space_representation = false;
+
this->v[1]->real_space_representation = true;
/* diagonal terms 11 22 33 */
this->v[1]->RLOOP (
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
//ptrdiff_t tindex = 3*rindex;
for (int cc=0; cc<3; cc++)
this->v[1]->rval(rindex,cc) = this->u->rval(rindex,cc)*this->u->rval(rindex,cc);
@@ -612,11 +624,11 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
this->v[1]->dft();
this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2 && k2 > 0)
{
ptrdiff_t tindex = 3*cindex;
@@ -635,10 +647,10 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
/* off-diagonal terms 12 23 31 */
this->v[1]->real_space_representation = true;
this->v[1]->RLOOP (
- [&](ptrdiff_t rindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex){
+ [&](const ptrdiff_t rindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex){
//ptrdiff_t tindex = 3*rindex;
for (int cc=0; cc<3; cc++)
this->v[1]->rval(rindex,cc) = this->u->rval(rindex,cc)*this->u->rval(rindex,(cc+1)%3);
@@ -648,11 +660,11 @@ void vorticity_equation<rnumber, be>::compute_pressure(field<rnumber, be, ONE> *
this->v[1]->dft();
this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2 && k2 > 0)
{
ptrdiff_t tindex = 3*cindex;
@@ -698,11 +710,11 @@ void vorticity_equation<rnumber, be>::compute_Lagrangian_acceleration(
acceleration->real_space_representation = false;
*acceleration = 0.0;
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
ptrdiff_t tindex = 3*cindex;
@@ -741,11 +753,11 @@ void vorticity_equation<rnumber, be>::compute_Eulerian_acceleration(
acceleration->real_space_representation = false;
/* put in linear terms */
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
ptrdiff_t tindex = 3*cindex;
@@ -786,11 +798,11 @@ void vorticity_equation<rnumber, be>::compute_Eulerian_acceleration(
this->v[1]->dft();
this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
ptrdiff_t tindex = 3*cindex;
@@ -825,11 +837,11 @@ void vorticity_equation<rnumber, be>::compute_Eulerian_acceleration(
this->v[1]->dft();
this->kk->template dealias<rnumber, THREE>(this->v[1]->get_cdata());
this->kk->CLOOP_K2(
- [&](ptrdiff_t cindex,
- ptrdiff_t xindex,
- ptrdiff_t yindex,
- ptrdiff_t zindex,
- double k2){
+ [&](const ptrdiff_t cindex,
+ const ptrdiff_t xindex,
+ const ptrdiff_t yindex,
+ const ptrdiff_t zindex,
+ const double k2){
if (k2 <= this->kk->kM2)
{
ptrdiff_t tindex = 3*cindex;
diff --git a/cpp/vorticity_equation.hpp b/cpp/vorticity_equation.hpp
index 09119cc48db3c37c4faf3c54162dbe5ebaea99e9..a7e3531555354d86df5e6ea78eff4c53fe672cff 100644
--- a/cpp/vorticity_equation.hpp
+++ b/cpp/vorticity_equation.hpp
@@ -3,20 +3,20 @@
* Copyright 2015 Max Planck Institute *
* for Dynamics and Self-Organization *
* *
-* This file is part of bfps. *
+* This file is part of TurTLE. *
* *
-* bfps is free software: you can redistribute it and/or modify *
+* TurTLE is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published *
* by the Free Software Foundation, either version 3 of the License, *
* or (at your option) any later version. *
* *
-* bfps is distributed in the hope that it will be useful, *
+* TurTLE is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
-* along with bfps. If not, see <http://www.gnu.org/licenses/> *
+* along with TurTLE. If not, see <http://www.gnu.org/licenses/> *
* *
* Contact: Cristian.Lalescu@ds.mpg.de *
* *
@@ -84,7 +84,7 @@ class vorticity_equation
double DKY = 1.0,
double DKZ = 1.0,
unsigned FFTW_PLAN_RIGOR = FFTW_MEASURE);
- ~vorticity_equation(void);
+ virtual ~vorticity_equation(void);
/* solver essential methods */
virtual void omega_nonlin(int src);
@@ -102,13 +102,13 @@ class vorticity_equation
field<rnumber, be, THREE> *src_vorticity);
void add_Kolmogorov_forcing(field<rnumber, be, THREE> *dst,
- int fmode,
- double famplitude);
+ const int fmode,
+ const double famplitude);
void add_field_band(
field<rnumber, be, THREE> *dst,
field<rnumber, be, THREE> *src,
- double k0, double k1,
- double prefactor);
+ const double k0, const double k1,
+ const double prefactor);
/** \brief Method that imposes action of forcing on new vorticity field.
*
diff --git a/documentation/chapters/cpp_doxygen.rst b/documentation/chapters/cpp_doxygen.rst
index 1e4246a162083d4508f590d7370d2d29028cd598..51d1c9cca6c637fbce1335a3b8d6c358ddb6a2fd 100644
--- a/documentation/chapters/cpp_doxygen.rst
+++ b/documentation/chapters/cpp_doxygen.rst
@@ -21,3 +21,13 @@ CPP
.. doxygenclass:: NSVE
:project: TurTLE
:members:
+
+.. doxygenclass:: particles_distr_mpi
+ :project: TurTLE
+ :path: ...
+ :members: [...]
+ :protected-members:
+ :private-members:
+ :undoc-members:
+ :outline:
+ :no-link:
diff --git a/documentation/chapters/overview.rst b/documentation/chapters/overview.rst
index 2ed997037f07026be50ee2976781c5eb7f861f92..c710c59b5d7117ad84a514cb1df14e88dbcc1a22 100644
--- a/documentation/chapters/overview.rst
+++ b/documentation/chapters/overview.rst
@@ -8,44 +8,37 @@ General comments
The purpose of this code is to run pseudo-spectral DNS of turbulence,
and integrate particle trajectories in the resulting fields.
-In brief, the main aim of the code is to simplify the launching of
+An important aim of the code is to simplify the launching of
compute jobs and postprocessing, up to and including the generation of
publication-ready figures.
For research, people routinely write code from scratch because research
goals change to a point where modifying the previous code is too
expensive.
-With bfps, the desire is to identify core functionality that should be
+With TurTLE, the desire is to identify core functionality that should be
implemented in a library.
The core library can then be used by many problem-specific codes.
In this sense, the structuring of the code-base is non-standard.
The core functionality is implemented in C++ (classes useful for
-describing working with fields or sets of particles), while a python
+describing working with fields or sets of particles), while a Python3
wrapper is used for generating "main" programmes to be linked against
the core library.
-The core library uses MPI for parallelization, and the python wrapper
+The core library uses a hybrid MPI/OpenMP approach for parallelization, and the Python3 wrapper
compiles this core library when being installed.
-The compilation environment can be configured for different
-machines as required.
Python3 "wrapper"
-----------------
-In principle, users of the code should only need to use python3 for
+In principle, users of the code should only need to use Python3 for
launching jobs and postprocessing data.
-While python2 compatibility should not be too hard to maintain, the
-usage of strings makes it a bit cumbersome ---
-the code makes extensive usage of strings for `HDF5` I/O.
-Classes defined in the python package can be used to generate executable
+Classes defined in the Python3 package can be used to generate executable
codes, compile/launch them, and then for accessing and postprocessing
-data.
-Obviously, postprocessing methods can be optimized with C extensions or
-otherwise, as needed.
+data with a full Python3 environment.
Code generation is quite straightforward, with C++ code snippets handled
-as strings in the python code, such that they can be combined in
+as strings in the Python3 code, such that they can be combined in
different ways.
Once a "main" file has been written, it is compiled and linked against
diff --git a/pc_host_info.py b/pc_host_info.py
index 383a3bc3abf0660f300bbcd0bc6e006afa4b1724..8e66eb3535c0d13f12d1434f5407fb068f5683fc 100644
--- a/pc_host_info.py
+++ b/pc_host_info.py
@@ -32,7 +32,9 @@ host_info = {'type' : 'pc'}
# 'deltanprocs' : info_template_deltanprocs,
# 'mail_address': info_template_mail_address,
# 'account' : info_template_account,
-# 'executable_launcher': info_template_executable_launcher}
+# 'executable_launcher' : info_template_executable_launcher,
+# 'extra_slurm_lines' : info_template_extra_slurm_lines,
+# 'explicit_slurm_environment' : info_template_explicit_slurm_environment}
# info_template_type can be one of:
# 'pc' --- jobs run interactively
@@ -59,3 +61,13 @@ host_info = {'type' : 'pc'}
# info_template_executable_launcher, relevant for 'SLURM' clusters,
# is the binary that should be called to execute the hybrid MPI/OpenMP code.
# by default it is 'srun', but it may be 'mpiexec' for some machines.
+
+# info_template_extra_slurm_lines, relevant for 'SLURM' clusters,
+# is a list of strings, written as extra options within the SLURM file.
+# for example something like `#SBATCH --get-user-env`
+
+# info_template_explicit_slurm_environment, relevant for `SLURM` clusters,
+# is a bool specifying whether the environment must be mentioned explicitly
+# in the submission script, or whether the cluster chooses it based on the
+# resources requested
+
diff --git a/tests/run_all_tests.sh b/tests/run_all_tests.sh
index 0d4ada45e44dcbff7da7e74a934dd36a3405a76e..227df5ee533d37d7039b70cfa41ed892f1bd3e89 100644
--- a/tests/run_all_tests.sh
+++ b/tests/run_all_tests.sh
@@ -10,6 +10,6 @@ turtle DNS static_field_with_ghost_collisions --simname dumb_test_of_ghost_colli
# test postprocessing
turtle PP field_single_to_double --simname dns_nsveparticles --iter0 32 --iter1 32
-turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64
+turtle PP get_rfields --simname dns_nsveparticles --iter0 0 --iter1 64 --TrS2_on 1
turtle PP joint_acc_vel_stats --simname dns_nsveparticles --iter0 0 --iter1 64
turtle PP resize --simname dns_nsveparticles --new_nx 96 --new_ny 96 --new_nz 96 --new_simname dns_nsveparticles_resized