diff --git a/CMakeLists.txt b/CMakeLists.txt
index 0e0bb22ca538e6a1fb1c835403a78d42b78c0e80..2abff7d03cac86ee627842edbed532f0adf64e77 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -262,11 +262,13 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/static_field.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/kraichnan_field.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_test.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/write_filtered_test.cpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/Gauss_field_test.cpp
${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
@@ -308,9 +310,11 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/direct_numerical_simulation.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVE.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/static_field.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/kraichnan_field.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/joint_acc_vel_stats.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/filter_test.hpp
+ ${PROJECT_SOURCE_DIR}/cpp/full_code/Gauss_field_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/field_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/write_filtered_test.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/symmetrize_test.hpp
diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index 1d5bea0ca4cb3e01b03b8562bd006d3af2c2ced3..cd1bb79b968bb3855af32ccc05b90dddd7f5fb37 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -439,7 +439,7 @@ 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']:
+ if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles', 'static_field', '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)
@@ -627,16 +627,10 @@ class DNS(_code):
parser_NSVE = subparsers.add_parser(
'NSVE',
help = 'plain Navier-Stokes vorticity formulation')
- self.simulation_parser_arguments(parser_NSVE)
- self.job_parser_arguments(parser_NSVE)
- self.parameters_to_parser_arguments(parser_NSVE)
parser_NSVE_no_output = subparsers.add_parser(
'NSVE_no_output',
help = 'plain Navier-Stokes vorticity formulation, checkpoints are NOT SAVED')
- self.simulation_parser_arguments(parser_NSVE_no_output)
- self.job_parser_arguments(parser_NSVE_no_output)
- self.parameters_to_parser_arguments(parser_NSVE_no_output)
parser_NSVEparticles_no_output = subparsers.add_parser(
'NSVEparticles_no_output',
@@ -646,6 +640,10 @@ class DNS(_code):
'static_field',
help = 'static field with basic fluid tracers')
+ parser_kraichnan_field = subparsers.add_parser(
+ 'kraichnan_field',
+ help = 'Kraichnan field with basic fluid tracers')
+
parser_NSVEp2 = subparsers.add_parser(
'NSVEparticles',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers')
@@ -653,22 +651,47 @@ class DNS(_code):
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.')
- parser_NSVE_ou = subparsers.add_parser(
- 'NSVE_ou_forcing',
- help = 'plain Navier-Stokes vorticity formulation, with ornstein-uhlenbeck forcing')
- for parser in ['NSVEparticles_no_output', 'NSVEp2', 'NSVEp2p', 'NSVEp_extra', 'static_field']:
+ for parser in [
+ 'NSVE',
+ 'NSVE_no_output',
+ 'NSVEparticles_no_output',
+ 'NSVEp2',
+ 'NSVEp2p',
+ 'NSVEp_extra',
+ 'static_field',
+ 'kraichnan_field']:
eval('self.simulation_parser_arguments({0})'.format('parser_' + parser))
eval('self.job_parser_arguments({0})'.format('parser_' + parser))
- eval('self.particle_parser_arguments({0})'.format('parser_' + parser))
eval('self.parameters_to_parser_arguments({0})'.format('parser_' + parser))
+ eval('self.parameters_to_parser_arguments('
+ 'parser_{0},'
+ 'self.generate_extra_parameters(\'{0}\'))'.format(parser))
+ for parser in [
+ 'NSVEparticles_no_output',
+ 'NSVEp2',
+ 'NSVEp2p',
+ 'NSVEp_extra',
+ 'static_field',
+ 'kraichnan_field']:
+ eval('self.particle_parser_arguments({0})'.format('parser_' + parser))
eval('self.parameters_to_parser_arguments('
'parser_{0},'
'self.NSVEp_extra_parameters)'.format(parser))
return None
+ def generate_extra_parameters(
+ self,
+ dns_type):
+ pars = {}
+ if dns_type == 'kraichnan_field':
+ pars['output_velocity'] = int(1)
+ pars['field_random_seed'] = int(1)
+ pars['spectrum_slope'] = float(-5./3)
+ pars['spectrum_k_cutoff'] = float(16)
+ pars['spectrum_coefficient'] = float(1)
+ return pars
def prepare_launch(
self,
args = [],
@@ -703,13 +726,16 @@ class DNS(_code):
self.dns_type = opt.DNS_class
self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
# merge parameters if needed
- if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling', 'static_field']:
+ if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling', 'static_field', 'kraichnan_field']:
for k in self.NSVEp_extra_parameters.keys():
self.parameters[k] = self.NSVEp_extra_parameters[k]
if type(extra_parameters) != type(None):
if self.dns_type in extra_parameters.keys():
for k in extra_parameters[self.dns_type].keys():
self.parameters[k] = extra_parameters[self.dns_type][k]
+ additional_parameters = self.generate_extra_parameters(self.dns_type)
+ for k in additional_parameters.keys():
+ self.parameters[k] = additional_parameters[k]
if ((self.parameters['niter_todo'] % self.parameters['niter_out']) != 0):
self.parameters['niter_out'] = self.parameters['niter_todo']
if len(opt.src_work_dir) == 0:
@@ -739,7 +765,10 @@ class DNS(_code):
opt.nz > opt.n):
opt.n = min(opt.nx, opt.ny, opt.nz)
print("Warning: '-n' parameter changed to minimum of nx, ny, nz. This affects the computation of nu.")
- self.parameters['dt'] = (opt.dtfactor / opt.n)
+ if self.dns_type in ['kraichnan_field']:
+ self.parameters['dt'] = opt.dtfactor * 0.5 / opt.n**2
+ else:
+ self.parameters['dt'] = (opt.dtfactor / opt.n)
self.parameters['nu'] = (opt.kMeta * 2 / opt.n)**(4./3)
# check value of kMax
kM = opt.n * 0.5
@@ -768,7 +797,7 @@ 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']:
+ if self.dns_type in ['kraichnan_field', 'static_field', 'NSVEparticles', 'NSVEcomplex_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
@@ -1004,19 +1033,27 @@ class DNS(_code):
# take care of fields' initial condition
# first, check if initial field exists
need_field = False
- if not os.path.exists(self.get_checkpoint_0_fname()):
- need_field = True
- else:
- f = h5py.File(self.get_checkpoint_0_fname(), 'r')
- try:
- dset = f['vorticity/complex/0']
- need_field = (dset.shape == (self.parameters['ny'],
- self.parameters['nz'],
- self.parameters['nx']//2+1,
- 3))
- except:
+ if self.dns_type in [
+ 'NSVE',
+ 'NSVE_no_output',
+ 'static_field',
+ 'NSVEparticles',
+ 'NSVEcomplex_particles',
+ 'NSVEparticles_no_output',
+ 'NSVEp_extra_sampling']:
+ if not os.path.exists(self.get_checkpoint_0_fname()):
need_field = True
- f.close()
+ else:
+ f = h5py.File(self.get_checkpoint_0_fname(), 'r')
+ try:
+ dset = f['vorticity/complex/0']
+ need_field = (dset.shape == (self.parameters['ny'],
+ self.parameters['nz'],
+ self.parameters['nx']//2+1,
+ 3))
+ except:
+ need_field = True
+ f.close()
if need_field:
f = h5py.File(self.get_checkpoint_0_fname(), 'a')
if len(opt.src_simname) > 0:
@@ -1043,8 +1080,19 @@ class DNS(_code):
amplitude = 0.05)
f['vorticity/complex/{0}'.format(0)] = data
f.close()
+ if self.dns_type == 'kraichnan_field':
+ if not os.path.exists(self.get_checkpoint_0_fname()):
+ f = h5py.File(self.get_checkpoint_0_fname(), 'a')
+ f.create_group('velocity/real')
+ f.close()
# now take care of particles' initial condition
- if self.dns_type in ['static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']:
+ if self.dns_type in [
+ 'kraichnan_field',
+ 'static_field',
+ 'NSVEparticles',
+ 'NSVEcomplex_particles',
+ 'NSVEparticles_no_output',
+ 'NSVEp_extra_sampling']:
self.generate_particle_data(opt = opt)
return None
def launch_jobs(
@@ -1053,6 +1101,9 @@ class DNS(_code):
if not os.path.exists(self.get_data_file_name()):
self.generate_initial_condition(opt = opt)
self.write_par()
+ if (('test' in self.dns_type) or
+ (self.dns_type in ['kraichnan_field'])):
+ self.check_current_vorticity_exists = False
self.run(
nb_processes = opt.nb_processes,
nb_threads_per_process = opt.nb_threads_per_process,
diff --git a/TurTLE/TEST.py b/TurTLE/TEST.py
index 14f1e138ab3c2ce411d81990faa351557d820e94..4bf762e7ab9ff1159929fcee46a9247a3c2dff44 100644
--- a/TurTLE/TEST.py
+++ b/TurTLE/TEST.py
@@ -51,6 +51,7 @@ class TEST(_code):
work_dir = work_dir,
simname = simname)
self.generate_default_parameters()
+ self.check_current_vorticity_exists = False
return None
def set_precision(
self,
@@ -118,7 +119,7 @@ class TEST(_code):
self.parameters['dky'] = float(1.0)
self.parameters['dkz'] = float(1.0)
self.parameters['filter_length'] = float(1.0)
- self.parameters['random_seed'] = int(1)
+ self.parameters['field_random_seed'] = int(1)
return None
def generate_extra_parameters(
self,
@@ -129,6 +130,12 @@ class TEST(_code):
pars['tracers0_integration_steps'] = int(4)
pars['tracers0_neighbours'] = int(1)
pars['tracers0_smoothness'] = int(1)
+ if dns_type == 'Gauss_field_test':
+ pars['histogram_bins'] = int(129)
+ pars['max_velocity_estimate'] = float(8)
+ pars['spectrum_slope'] = float(-5./3)
+ pars['spectrum_k_cutoff'] = float(16)
+ pars['spectrum_coefficient'] = float(1)
return pars
def get_kspace(self):
kspace = {}
@@ -175,6 +182,83 @@ class TEST(_code):
kspace = self.get_kspace()
nshells = kspace['nshell'].shape[0]
ofile['checkpoint'] = int(0)
+ vec_spectra_stats = []
+ tens_rspace_stats = []
+ vec4_rspace_stats = []
+ scal_rspace_stats = []
+ if self.dns_type in ['Gauss_field_test']:
+ vec_spectra_stats.append('velocity')
+ vec4_rspace_stats.append('velocity')
+ tens_rspace_stats.append('velocity_gradient')
+ scal_rspace_stats.append('velocity_divergence')
+ for k in vec_spectra_stats:
+ time_chunk = 2**20//(8*3*3*nshells)
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/spectra/' + k + '_' + k,
+ (1, nshells, 3, 3),
+ chunks = (time_chunk, nshells, 3, 3),
+ maxshape = (None, nshells, 3, 3),
+ dtype = np.float64)
+ for k in scal_rspace_stats:
+ time_chunk = 2**20//(8*10)
+ time_chunk = max(time_chunk, 1)
+ a = ofile.create_dataset('statistics/moments/' + k,
+ (1, 10),
+ chunks = (time_chunk, 10),
+ maxshape = (None, 10),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*self.parameters['histogram_bins'])
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/histograms/' + k,
+ (1,
+ self.parameters['histogram_bins']),
+ chunks = (time_chunk,
+ self.parameters['histogram_bins']),
+ maxshape = (None,
+ self.parameters['histogram_bins']),
+ dtype = np.int64)
+ for k in vec4_rspace_stats:
+ time_chunk = 2**20//(8*4*10)
+ time_chunk = max(time_chunk, 1)
+ a = ofile.create_dataset('statistics/moments/' + k,
+ (1, 10, 4),
+ chunks = (time_chunk, 10, 4),
+ maxshape = (None, 10, 4),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*4*self.parameters['histogram_bins'])
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/histograms/' + k,
+ (1,
+ self.parameters['histogram_bins'],
+ 4),
+ chunks = (time_chunk,
+ self.parameters['histogram_bins'],
+ 4),
+ maxshape = (None,
+ self.parameters['histogram_bins'],
+ 4),
+ dtype = np.int64)
+ for k in tens_rspace_stats:
+ time_chunk = 2**20//(8*9*10)
+ time_chunk = max(time_chunk, 1)
+ a = ofile.create_dataset('statistics/moments/' + k,
+ (1, 10, 3, 3),
+ chunks = (time_chunk, 10, 3, 3),
+ maxshape = (None, 10, 3, 3),
+ dtype = np.float64)
+ time_chunk = 2**20//(8*9*self.parameters['histogram_bins'])
+ time_chunk = max(time_chunk, 1)
+ ofile.create_dataset('statistics/histograms/' + k,
+ (1,
+ self.parameters['histogram_bins'],
+ 3, 3),
+ chunks = (time_chunk,
+ self.parameters['histogram_bins'],
+ 3, 3),
+ maxshape = (None,
+ self.parameters['histogram_bins'],
+ 3, 3),
+ dtype = np.int64)
return None
def job_parser_arguments(
self,
@@ -257,6 +341,9 @@ class TEST(_code):
dest = 'TEST_class',
help = 'type of simulation to run')
subparsers.required = True
+ parser_Gauss_field_test = subparsers.add_parser(
+ 'Gauss_field_test',
+ help = 'test incompressible Gaussian random fields')
parser_filter_test = subparsers.add_parser(
'filter_test',
help = 'plain filter test')
@@ -275,23 +362,21 @@ class TEST(_code):
parser_test_interpolation = subparsers.add_parser(
'test_interpolation',
help = 'test velocity gradient interpolation')
-
parser_ornstein_uhlenbeck_test = subparsers.add_parser(
'ornstein_uhlenbeck_test',
help = 'test ornstein uhlenbeck process')
- for parser in ['parser_filter_test',
- 'parser_write_filtered_test',
- 'parser_field_test',
- 'parser_symmetrize_test',
- 'parser_field_output_test',
- 'parser_test_interpolation',
- 'parser_ornstein_uhlenbeck_test']:
- eval('self.simulation_parser_arguments(' + parser + ')')
- eval('self.job_parser_arguments(' + parser + ')')
- eval('self.parameters_to_parser_arguments(' + parser + ')')
- eval('self.parameters_to_parser_arguments(' + parser + ',' +
- 'parameters = self.generate_extra_parameters(dns_type = \'' + parser + '\'))')
+ for parser in ['Gauss_field_test',
+ 'filter_test',
+ 'write_filtered_test',
+ 'field_test',
+ 'symmetrize_test',
+ 'field_output_test',
+ 'test_interpolation',
+ 'ornstein_uhlenbeck_test']:
+ eval('self.simulation_parser_arguments(parser_' + parser + ')')
+ eval('self.job_parser_arguments(parser_' + parser + ')')
+ eval('self.parameters_to_parser_arguments(parser_' + parser + ')')
return None
def prepare_launch(
self,
@@ -300,6 +385,8 @@ class TEST(_code):
self.set_precision(opt.precision)
self.dns_type = opt.TEST_class
self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
+ self.parameters.update(
+ self.generate_extra_parameters(dns_type = self.dns_type))
# merge parameters if needed
self.pars_from_namespace(opt)
return opt
@@ -308,8 +395,6 @@ class TEST(_code):
args = [],
**kwargs):
opt = self.prepare_launch(args = args)
- self.parameters.update(
- self.generate_extra_parameters(dns_type = self.dns_type))
self.launch_jobs(opt = opt, **kwargs)
return None
def launch_jobs(
diff --git a/TurTLE/_code.py b/TurTLE/_code.py
index 05b36518ae81205a8e6af088402d73b41add25b3..be898972439e485f3ad3f6b3e4a947c28bc14951 100644
--- a/TurTLE/_code.py
+++ b/TurTLE/_code.py
@@ -176,6 +176,7 @@ class _code(_base):
"""
self.host_info = host_info
self.main = ''
+ self.check_current_vorticity_exists = True
return None
def write_src(self):
with open(self.name + '.cpp', 'w') as outfile:
@@ -278,7 +279,7 @@ class _code(_base):
if os.path.exists(os.path.join(self.work_dir, 'stop_' + self.simname)):
warnings.warn("Found stop_<simname> file, I will remove it.")
os.remove(os.path.join(self.work_dir, 'stop_' + self.simname))
- if not 'test' in self.name:
+ if self.check_current_vorticity_exists:
with h5py.File(os.path.join(self.work_dir, self.simname + '.h5'), 'r') as data_file:
iter0 = data_file['iteration'][...]
checkpoint0 = data_file['checkpoint'][()]
diff --git a/TurTLE/test/test_Gaussian_field.py b/TurTLE/test/test_Gaussian_field.py
new file mode 100644
index 0000000000000000000000000000000000000000..149073e7c1abe931368df1a026ef8f02ac4bc070
--- /dev/null
+++ b/TurTLE/test/test_Gaussian_field.py
@@ -0,0 +1,109 @@
+#! /usr/bin/env python
+
+import numpy as np
+from scipy import trapz
+from scipy.stats import norm
+from scipy.integrate import quad
+import h5py
+import sys
+import time
+
+import TurTLE
+from TurTLE import TEST
+
+try:
+ import matplotlib.pyplot as plt
+except:
+ plt = None
+
+def main():
+ c = TEST()
+ # size of grid
+ n = 256
+ slope = -5./3.
+ k_cutoff = 30.
+ 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]
+ coeff = 1./total_energy
+ bin_no = 100
+ rseed = int(time.time())
+
+ c.launch(
+ ['Gauss_field_test',
+ '--nx', str(n),
+ '--ny', str(n),
+ '--nz', str(n),
+ '--simname', 'Gaussianity_test',
+ '--np', '4',
+ '--ntpp', '1',
+ '--wd', './',
+ '--histogram_bins', str(bin_no),
+ '--max_velocity_estimate', '8.',
+ '--spectrum_slope', str(slope),
+ '--spectrum_k_cutoff', str(k_cutoff),
+ '--spectrum_coefficient', str(coeff),
+ '--field_random_seed', str(rseed)] +
+ sys.argv[1:])
+ plot_stuff(c.simname, total_energy = total_energy)
+ return None
+
+def plot_stuff(simname, total_energy = 1.):
+ df = h5py.File(simname + '.h5', 'r')
+ for kk in ['spectrum_slope',
+ 'spectrum_k_cutoff',
+ 'spectrum_coefficient',
+ 'field_random_seed',
+ 'histogram_bins']:
+ print(kk, df['parameters/' + kk][...])
+ coeff = df['parameters/spectrum_coefficient'][()]
+ k_cutoff = df['parameters/spectrum_k_cutoff'][()]
+ slope = df['parameters/spectrum_slope'][()]
+ bin_no = df['parameters/histogram_bins'][()]
+
+ f = plt.figure()
+ # test spectrum
+ a = f.add_subplot(121)
+ kk = df['kspace/kshell'][...]
+ print('dk: {}'.format(df['kspace/dk'][()]))
+ phi_ij = df['statistics/spectra/velocity_velocity'][0] / df['kspace/dk'][()]
+ energy = (phi_ij[..., 0, 0] + phi_ij[..., 1, 1] + phi_ij[..., 2, 2])/2
+
+ a.scatter(kk[1:-2], energy[1:-2])
+ a.plot(kk[1:-2], coeff*kk[1:-2]**slope*np.exp(-kk[1:-2]/k_cutoff), ls='--', c='C0')
+ a.set_xscale('log')
+ a.set_yscale('log')
+ # test isotropy
+ a = f.add_subplot(122)
+
+ max_vel_estimate = df['parameters/max_velocity_estimate'][()]
+ velbinsize = 2*max_vel_estimate/bin_no
+ vel = np.linspace(-max_vel_estimate+velbinsize/2, max_vel_estimate-velbinsize/2, bin_no)
+ hist_vel = df['statistics/histograms/velocity'][0, :, :3]
+ f_vel = hist_vel / np.sum(hist_vel, axis=0, keepdims=True).astype(float) / velbinsize
+
+ print('Energy analytically: {}'.format(total_energy))
+ 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)))
+
+ vel_variance = df['statistics/moments/velocity'][0,2,3]/3.
+ a.plot(vel[:,None]/np.sqrt(vel_variance), f_vel*np.sqrt(vel_variance))
+ a.plot(vel/np.sqrt(vel_variance), norm.pdf(vel/np.sqrt(vel_variance)), ls='--')
+ a.set_yscale('log')
+ a.set_xlim(-5,5)
+ a.set_ylim(1e-5,1)
+ f.tight_layout()
+ f.savefig('spectrum_isotropy_test.pdf')
+ plt.close(f)
+
+ ### check divergence
+ print('Divergence second moment is: {0}'.format(
+ df['statistics/moments/velocity_divergence'][0, 2]))
+ print('Gradient second moment is: {0}'.format(
+ df['statistics/moments/velocity_gradient'][0, 2].mean()))
+ df.close()
+ return None
+
+if __name__ == '__main__':
+ main()
+
diff --git a/cpp/field.cpp b/cpp/field.cpp
index 04de2b0cfeef1992117b8493953dda6cff9520a8..81f71e10d330ea0886091d13348584378e01011c 100644
--- a/cpp/field.cpp
+++ b/cpp/field.cpp
@@ -29,6 +29,7 @@
#include <cstdlib>
#include <algorithm>
#include <cassert>
+#include <random>
#include "field.hpp"
#include "scope_timer.hpp"
#include "shared_array.hpp"
@@ -1500,7 +1501,7 @@ int compute_gradient(
assert(!src->real_space_representation);
assert((fc1 == ONE && fc2 == THREE) ||
(fc1 == THREE && fc2 == THREExTHREE));
- std::fill_n(dst->get_rdata(), dst->rmemlayout->local_size, 0);
+ *dst = 0.0;
dst->real_space_representation = false;
switch(fc1)
{
@@ -1547,6 +1548,34 @@ int compute_gradient(
return EXIT_SUCCESS;
}
+template <typename rnumber,
+ field_backend be,
+ kspace_dealias_type dt>
+int compute_divergence(
+ kspace<be, dt> *kk,
+ field<rnumber, be, THREE> *src,
+ field<rnumber, be, ONE> *dst)
+{
+ TIMEZONE("compute_divergence");
+ assert(!src->real_space_representation);
+ *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){
+ 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));
+ dst->cval(cindex, 1) = (kk->kx[xindex]*src->cval(cindex, 0, 0) +
+ kk->ky[yindex]*src->cval(cindex, 1, 0) +
+ kk->kz[zindex]*src->cval(cindex, 2, 0));
+ });
+ return EXIT_SUCCESS;
+}
+
template <typename rnumber,
field_backend be,
kspace_dealias_type dt>
@@ -2075,6 +2104,77 @@ field<rnumber, be, fc> &field<rnumber, be, fc>::operator=(
return *this;
}
+template <typename rnumber,
+ field_backend be,
+ field_components fc,
+ kspace_dealias_type dt>
+int make_gaussian_random_field(
+ kspace<be, dt> *kk,
+ field<rnumber, be, fc> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient)
+{
+ TIMEZONE("make_gaussian_random_field");
+ // initialize a separate random number generator for each thread
+ std::vector<std::mt19937_64> rgen;
+ std::normal_distribution<rnumber> rdist;
+ rgen.resize(omp_get_max_threads());
+ // seed random number generators such that no seed is ever repeated if we change the value of rseed.
+ // basically use a multi-dimensional array indexing technique to come up with actual seed.
+ // Note: this method IS NOT MPI/OpenMP-invariant!
+ for (int thread_id=0; thread_id < omp_get_max_threads(); thread_id++)
+ {
+ int current_seed = (
+ 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);
+ rgen[thread_id].seed(current_seed);
+ }
+ output_field->real_space_representation = false;
+ *output_field = 0.0;
+ 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)
+ {
+ if (k2 > 0)
+ switch(fc)
+ {
+ case ONE:
+ {
+ output_field->cval(cindex,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ output_field->cval(cindex,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ break;
+ }
+ case THREE:
+ for (int cc = 0; cc<3; cc++)
+ {
+ output_field->cval(cindex,cc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ output_field->cval(cindex,cc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ }
+ break;
+ case THREExTHREE:
+ for (int cc = 0; cc<3; cc++)
+ for (int ccc = 0; ccc<3; ccc++)
+ {
+ output_field->cval(cindex,cc,ccc,0) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ output_field->cval(cindex,cc,ccc,1) = rdist(rgen[omp_get_thread_num()]) * sqrt(coefficient * pow(k2, (slope - 2.)/2.) * exp(-sqrt(k2)/k_cutoff) / 12 / M_PI);
+ }
+ break;
+ }
+ });
+ output_field->symmetrize();
+
+ return EXIT_SUCCESS;
+}
+
template class field<float, FFTW, ONE>;
template class field<float, FFTW, THREE>;
template class field<float, FFTW, THREExTHREE>;
@@ -2168,6 +2268,15 @@ template int compute_gradient<double, FFTW, ONE, THREE, SMOOTH>(
field<double, FFTW, ONE> *,
field<double, FFTW, THREE> *);
+template int compute_divergence(
+ kspace<FFTW, SMOOTH> *,
+ field<float, FFTW, THREE> *,
+ field<float, FFTW, ONE> *);
+template int compute_divergence(
+ kspace<FFTW, SMOOTH> *,
+ field<double, FFTW, THREE> *,
+ field<double, FFTW, ONE> *);
+
template int invert_curl<float, FFTW, SMOOTH>(
kspace<FFTW, SMOOTH> *,
field<float, FFTW, THREE> *,
@@ -2232,3 +2341,49 @@ template int joint_rspace_3PDF<double, FFTW>(
const std::vector<double>,
const std::vector<double>);
+
+template int make_gaussian_random_field<
+ float,
+ FFTW,
+ ONE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<float, FFTW, ONE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+template int make_gaussian_random_field<
+ double,
+ FFTW,
+ ONE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<double, FFTW, ONE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+template int make_gaussian_random_field<
+ float,
+ FFTW,
+ THREE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<float, FFTW, THREE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+template int make_gaussian_random_field<
+ double,
+ FFTW,
+ THREE,
+ SMOOTH>(
+ kspace<FFTW, SMOOTH> *kk,
+ field<double, FFTW, THREE> *output_field,
+ const int rseed,
+ const double slope,
+ const double k_cutoff,
+ const double coefficient);
+
diff --git a/cpp/field.hpp b/cpp/field.hpp
index af2136b193c4b22d65972fac287c31ac7d20ec38..b1dbc5e03e8d9c8b8a306518956993adec4abe7e 100644
--- a/cpp/field.hpp
+++ b/cpp/field.hpp
@@ -325,7 +325,8 @@ class field
}
};
-/*
+/** \brief Compute gradient of a field
+ *
* if $A_{ij} = \partial_j u_i$
* then `compute_gradient` for a vector field stores the data in the resulting
* THREExTHREE field with $i$ the fast counter and $j$ the slow counter.
@@ -341,6 +342,17 @@ int compute_gradient(
field<rnumber, be, fc1> *source,
field<rnumber, be, fc2> *destination);
+/** \brief Compute divergence of a vector field
+ *
+ * */
+template <typename rnumber,
+ field_backend be,
+ kspace_dealias_type dt>
+int compute_divergence(
+ kspace<be, dt> *kk,
+ field<rnumber, be, THREE> *source,
+ field<rnumber, be, ONE> *destination);
+
template <typename rnumber,
field_backend be,
kspace_dealias_type dt>
@@ -374,5 +386,21 @@ int joint_rspace_3PDF(
const std::vector<double> max_f2_estimate,
const std::vector<double> max_f3_estimate);
+/** \brief Generate a Gaussian random field
+ *
+ * \note This method IS NOT MPI/OpenMP-invariant!
+ * */
+template <typename rnumber,
+ field_backend be,
+ field_components fc,
+ kspace_dealias_type dt>
+int make_gaussian_random_field(
+ kspace<be, dt> *kk,
+ field<rnumber, be, fc> *output_field,
+ const int rseed = 0,
+ const double slope = -5./3.,
+ const double k2_cutoff = 10.,
+ const double coefficient = 1.);
+
#endif//FIELD_HPP
diff --git a/cpp/full_code/Gauss_field_test.cpp b/cpp/full_code/Gauss_field_test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6822b6d8c02f87a532ffb0d0f63ed9bc767a8564
--- /dev/null
+++ b/cpp/full_code/Gauss_field_test.cpp
@@ -0,0 +1,182 @@
+/**********************************************************************
+* *
+* 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 *
+* *
+**********************************************************************/
+
+
+
+#include <string>
+#include <cmath>
+#include <random>
+#include "Gauss_field_test.hpp"
+#include "scope_timer.hpp"
+
+
+
+
+
+template <typename rnumber>
+int Gauss_field_test<rnumber>::initialize(void)
+{
+ TIMEZONE("Gauss_field_test::initialize");
+ this->read_parameters();
+ this->scalar_field = new field<rnumber, FFTW, ONE>(
+ nx, ny, nz,
+ this->comm,
+ FFTW_ESTIMATE);
+ this->vector_field = new field<rnumber, FFTW, THREE>(
+ nx, ny, nz,
+ this->comm,
+ FFTW_ESTIMATE);
+ this->kk = new kspace<FFTW, SMOOTH>(
+ this->scalar_field->clayout, this->dkx, this->dky, this->dkz);
+
+ if (this->myrank == 0)
+ {
+ hid_t stat_file = H5Fopen(
+ (this->simname + std::string(".h5")).c_str(),
+ H5F_ACC_RDWR,
+ H5P_DEFAULT);
+ this->kk->store(stat_file);
+ H5Fclose(stat_file);
+ }
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Gauss_field_test<rnumber>::finalize(void)
+{
+ TIMEZONE("Gauss_field_test::finalize");
+ delete this->vector_field;
+ delete this->scalar_field;
+ delete this->kk;
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Gauss_field_test<rnumber>::read_parameters()
+{
+ TIMEZONE("Gauss_field_test::read_parameters");
+ this->test::read_parameters();
+ hid_t parameter_file = H5Fopen(
+ (this->simname + std::string(".h5")).c_str(),
+ H5F_ACC_RDONLY,
+ H5P_DEFAULT);
+ this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "/parameters/max_velocity_estimate");
+ this->spectrum_slope = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_slope");
+ this->spectrum_k_cutoff = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_k_cutoff");
+ this->spectrum_coefficient = hdf5_tools::read_value<double>(parameter_file, "/parameters/spectrum_coefficient");
+ this->random_seed = hdf5_tools::read_value<int>(parameter_file, "/parameters/field_random_seed");
+ H5Fclose(parameter_file);
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int Gauss_field_test<rnumber>::do_work(void)
+{
+ TIMEZONE("Gauss_field_test::do_work");
+
+ /// initialize vector Gaussian random field
+ make_gaussian_random_field(
+ this->kk,
+ this->vector_field,
+ this->random_seed,
+ this->spectrum_slope,
+ this->spectrum_k_cutoff,
+ this->spectrum_coefficient);
+
+ /// impose divergence free condition while maintaining the energy of the field
+ this->kk->template project_divfree<rnumber>(this->vector_field->get_cdata(), true);
+
+ /// initialize statistics file.
+ hid_t stat_group, stat_file;
+ if (this->myrank == 0)
+ {
+ // set caching parameters
+ hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
+ herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
+ DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
+ stat_file = H5Fopen(
+ (this->simname + ".h5").c_str(),
+ H5F_ACC_RDWR,
+ fapl);
+ stat_group = H5Gopen(
+ stat_file,
+ "statistics",
+ H5P_DEFAULT);
+ }
+ else
+ {
+ stat_file = 0;
+ stat_group = 0;
+ }
+
+ /// compute gradient and divergence of field
+ field<rnumber, FFTW, THREExTHREE> *tensor_field = new field<rnumber, FFTW, THREExTHREE>(
+ nx, ny, nz,
+ this->comm,
+ FFTW_ESTIMATE);
+ compute_gradient(
+ this->kk,
+ this->vector_field,
+ tensor_field);
+ compute_divergence(
+ this->kk,
+ this->vector_field,
+ this->scalar_field);
+
+ /// compute basic statistics of Gaussian field
+ this->vector_field->compute_stats(
+ this->kk,
+ stat_group,
+ "velocity",
+ 0,
+ this->max_velocity_estimate);
+
+ /// verify divergence free method
+ tensor_field->ift();
+ /// for the stats max value doesn't really matter, I just want the moments
+ tensor_field->compute_rspace_stats(
+ stat_group,
+ "velocity_gradient",
+ 0,
+ std::vector<double>({1, 1, 1, 1, 1, 1, 1, 1, 1}));
+ this->scalar_field->ift();
+ this->scalar_field->compute_rspace_stats(
+ stat_group,
+ "velocity_divergence",
+ 0,
+ std::vector<double>({1}));
+ delete tensor_field;
+
+ /// close stat file
+ if (this->myrank == 0)
+ {
+ H5Gclose(stat_group);
+ H5Fclose(stat_file);
+ }
+ return EXIT_SUCCESS;
+}
+
+template class Gauss_field_test<float>;
+template class Gauss_field_test<double>;
+
diff --git a/cpp/full_code/Gauss_field_test.hpp b/cpp/full_code/Gauss_field_test.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..855b19279c27da7a058b5c010ef2197308d317ef
--- /dev/null
+++ b/cpp/full_code/Gauss_field_test.hpp
@@ -0,0 +1,74 @@
+/**********************************************************************
+* *
+* 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 GAUSS_FIELD_TEST_HPP
+#define GAUSS_FIELD_TEST_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "field.hpp"
+#include "full_code/test.hpp"
+
+/** \brief A test of incompressible Gaussian random fields.
+ *
+ */
+
+template <typename rnumber>
+class Gauss_field_test: public test
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ double max_velocity_estimate;
+ double spectrum_slope;
+ double spectrum_k_cutoff;
+ double spectrum_coefficient;
+ int random_seed;
+
+ /* other stuff */
+ kspace<FFTW, SMOOTH> *kk;
+ field<rnumber, FFTW, ONE> *scalar_field;
+ field<rnumber, FFTW, THREE> *vector_field;
+
+ Gauss_field_test(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ test(
+ COMMUNICATOR,
+ simulation_name){}
+ ~Gauss_field_test(){}
+
+ int initialize(void);
+ int do_work(void);
+ int finalize(void);
+ int read_parameters(void);
+};
+
+#endif//GAUSS_FIELD_TEST_HPP
+
diff --git a/cpp/full_code/kraichnan_field.cpp b/cpp/full_code/kraichnan_field.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..212ba9b8184019b0ab24b45336ac15f7e5cf66d5
--- /dev/null
+++ b/cpp/full_code/kraichnan_field.cpp
@@ -0,0 +1,377 @@
+/******************************************************************************
+* *
+* 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 *
+* *
+******************************************************************************/
+
+
+
+#include <string>
+#include <cmath>
+#include "kraichnan_field.hpp"
+#include "scope_timer.hpp"
+#include "fftw_tools.hpp"
+
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::initialize(void)
+{
+ TIMEZONE("kraichnan_file::initialize");
+ this->read_iteration();
+ this->read_parameters();
+ if (this->myrank == 0)
+ {
+ // set caching parameters
+ hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
+ herr_t cache_err = H5Pset_cache(fapl, 0, 521, 134217728, 1.0);
+ variable_used_only_in_assert(cache_err);
+ DEBUG_MSG("when setting stat_file cache I got %d\n", cache_err);
+ this->stat_file = H5Fopen(
+ (this->simname + ".h5").c_str(),
+ H5F_ACC_RDWR,
+ fapl);
+ }
+ int data_file_problem;
+ if (this->myrank == 0)
+ data_file_problem = this->grow_file_datasets();
+ MPI_Bcast(&data_file_problem, 1, MPI_INT, 0, this->comm);
+ if (data_file_problem > 0)
+ {
+ std::cerr <<
+ data_file_problem <<
+ " problems growing file datasets.\ntrying to exit now." <<
+ std::endl;
+ return EXIT_FAILURE;
+ }
+
+ this->velocity = new field<rnumber, FFTW, THREE>(
+ nx, ny, nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+ this->velocity->real_space_representation = true;
+
+ // the velocity layout should be the same as the vorticity one
+ this->kk = new kspace<FFTW, SMOOTH>(
+ this->velocity->clayout, this->dkx, this->dky, this->dkz);
+
+ // initialize particles
+ this->ps = particles_system_builder(
+ this->velocity, // (field object)
+ this->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->get_current_fname(), // particles input filename
+ std::string("/tracers0/state/") + std::to_string(this->iteration), // dataset name for initial input
+ std::string("/tracers0/rhs/") + std::to_string(this->iteration), // dataset name for initial input
+ tracers0_neighbours, // parameter (interpolation no neighbours)
+ tracers0_smoothness, // parameter
+ this->comm,
+ this->iteration+1);
+ // initialize output objects
+ this->particles_output_writer_mpi = new particles_output_hdf5<
+ long long int, double, 3>(
+ 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());
+
+ DEBUG_MSG("Coefficient: %g\n",
+ this->spectrum_coefficient);
+
+ // 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;
+}
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::generate_random_velocity(void)
+{
+ TIMEZONE("kraichnan_file::make_random_field");
+ // generate Gaussian random field
+ make_gaussian_random_field(
+ this->kk,
+ this->velocity,
+ this->iteration, // not an ideal choice because resulting field sequence depends on MPI/OpenMP configuration
+ 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
+ DEBUG_MSG("L2Norm before: %g\n",
+ this->velocity->L2norm(this->kk));
+ this->kk->template project_divfree<rnumber>(this->velocity->get_cdata());
+ DEBUG_MSG("L2Norm after: %g\n",
+ this->velocity->L2norm(this->kk));
+ // transform velocity to real space representation
+ this->velocity->ift();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::step(void)
+{
+ TIMEZONE("kraichnan_file::step");
+ this->ps->completeLoop(sqrt(this->dt));
+ this->generate_random_velocity();
+ this->iteration++;
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::write_checkpoint(void)
+{
+ TIMEZONE("kraichnan_file::write_checkpoint");
+ this->update_checkpoint();
+ // output field
+ if (this->output_velocity == 1)
+ {
+ assert(this->velocity->real_space_representation);
+ std::string fname = this->get_current_fname();
+ this->velocity->io(
+ fname,
+ "velocity",
+ this->iteration,
+ false);
+ }
+ // output particles
+ this->particles_output_writer_mpi->open_file(this->get_current_fname());
+ this->particles_output_writer_mpi->template save<3>(
+ this->ps->getParticlesState(),
+ this->ps->getParticlesRhs(),
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->iteration);
+ this->particles_output_writer_mpi->close_file();
+ this->write_iteration();
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::finalize(void)
+{
+ TIMEZONE("kraichnan_field::finalize");
+ if (this->myrank == 0)
+ H5Fclose(this->stat_file);
+ // good practice rule number n+1: always delete in inverse order of allocation
+ delete this->ps.release();
+ delete this->particles_output_writer_mpi;
+ delete this->particles_sample_writer_mpi;
+ delete this->kk;
+ delete this->velocity;
+ return EXIT_SUCCESS;
+}
+
+/** \brief Compute statistics.
+ *
+ */
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::do_stats()
+{
+ TIMEZONE("kraichnan_field::do_stats");
+ /// compute field statistics
+ if (this->iteration % this->niter_stat == 0)
+ {
+ hid_t stat_group;
+ if (this->myrank == 0)
+ stat_group = H5Gopen(
+ this->stat_file,
+ "statistics",
+ H5P_DEFAULT);
+ else
+ stat_group = 0;
+ field<rnumber, FFTW, THREE> *tvelocity = new field<rnumber, FFTW, THREE>(
+ this->nx, this->ny, this->nz,
+ this->comm,
+ fftw_planner_string_to_flag[this->fftw_plan_rigor]);
+ *tvelocity = *this->velocity;
+ tvelocity->compute_stats(
+ this->kk,
+ stat_group,
+ "velocity",
+ this->iteration / this->niter_stat,
+ this->max_velocity_estimate/sqrt(3));
+ if (this->myrank == 0)
+ H5Gclose(stat_group);
+ delete tvelocity;
+ }
+ /// 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
+ std::unique_ptr<double[]> pdata(new double[3*this->ps->getLocalNbParticles()]);
+
+ /// copy position data
+
+ /// sample position
+ std::copy(this->ps->getParticlesState(),
+ this->ps->getParticlesState()+3*this->ps->getLocalNbParticles(),
+ pdata.get());
+ this->particles_sample_writer_mpi->template save_dataset<3>(
+ "tracers0",
+ "position",
+ this->ps->getParticlesState(),
+ &pdata,
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->ps->get_step_idx()-1);
+
+ /// sample velocity
+ std::fill_n(pdata.get(), 3*this->ps->getLocalNbParticles(), 0);
+ this->ps->sample_compute_field(*this->velocity, pdata.get());
+ this->particles_sample_writer_mpi->template save_dataset<3>(
+ "tracers0",
+ "velocity",
+ this->ps->getParticlesState(),
+ &pdata,
+ this->ps->getParticlesIndexes(),
+ this->ps->getLocalNbParticles(),
+ this->ps->get_step_idx()-1);
+
+ // deallocate temporary data array
+ delete[] pdata.release();
+
+ return EXIT_SUCCESS;
+}
+
+template <typename rnumber>
+int kraichnan_field<rnumber>::read_parameters(void)
+{
+ TIMEZONE("kraichnan_field::read_parameters");
+ this->direct_numerical_simulation::read_parameters();
+ hid_t parameter_file = H5Fopen((this->simname + ".h5").c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+ this->fftw_plan_rigor = hdf5_tools::read_string(parameter_file, "parameters/fftw_plan_rigor");
+ this->dt = hdf5_tools::read_value<double>(parameter_file, "parameters/dt");
+ 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->spectrum_slope = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_slope");
+ this->spectrum_k_cutoff = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_k_cutoff");
+ this->spectrum_coefficient = hdf5_tools::read_value<double>(parameter_file, "parameters/spectrum_coefficient");
+ this->field_random_seed = hdf5_tools::read_value<int>(parameter_file, "parameters/field_random_seed");
+ this->output_velocity = hdf5_tools::read_value<int>(parameter_file, "parameters/output_velocity");
+ this->max_velocity_estimate = hdf5_tools::read_value<double>(parameter_file, "parameters/max_velocity_estimate");
+ H5Fclose(parameter_file);
+ return EXIT_SUCCESS;
+}
+
+template <class rnumber>
+void kraichnan_field<rnumber>::update_checkpoint()
+{
+ TIMEZONE("kraichnan_field::update_checkpoint");
+ DEBUG_MSG("entered kraichan_field::update_checkpoint\n");
+ std::string fname = this->get_current_fname();
+ if (this->kk->layout->myrank == 0)
+ {
+ bool file_exists = false;
+ {
+ struct stat file_buffer;
+ file_exists = (stat(fname.c_str(), &file_buffer) == 0);
+ }
+ if (file_exists)
+ {
+ // check how many fields there are in the checkpoint file
+ // increment checkpoint if needed
+ 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,
+ &fields_stored);
+ bool dset_exists = H5Lexists(
+ group_id,
+ std::to_string(this->iteration).c_str(),
+ H5P_DEFAULT);
+ H5Gclose(group_id);
+ H5Fclose(fid);
+ if ((int(fields_stored) >= this->checkpoints_per_file) &&
+ !dset_exists)
+ this->checkpoint++;
+ }
+ else
+ {
+ // create file, create fields_stored dset
+ hid_t fid = H5Fcreate(
+ fname.c_str(),
+ H5F_ACC_EXCL,
+ H5P_DEFAULT,
+ H5P_DEFAULT);
+ hid_t gg = H5Gcreate(
+ fid,
+ "velocity",
+ H5P_DEFAULT,
+ H5P_DEFAULT,
+ H5P_DEFAULT);
+ hid_t ggg = H5Gcreate(
+ gg,
+ "real",
+ H5P_DEFAULT,
+ H5P_DEFAULT,
+ H5P_DEFAULT);
+ H5Gclose(ggg);
+ H5Gclose(gg);
+ H5Fclose(fid);
+ }
+ }
+ MPI_Bcast(&this->checkpoint, 1, MPI_INT, 0, this->kk->layout->comm);
+ DEBUG_MSG("exiting kraichan_field::update_checkpoint\n");
+}
+
+template class kraichnan_field<float>;
+template class kraichnan_field<double>;
+
diff --git a/cpp/full_code/kraichnan_field.hpp b/cpp/full_code/kraichnan_field.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d7fb215efcdfa57a6ac22c03717526839e1ab016
--- /dev/null
+++ b/cpp/full_code/kraichnan_field.hpp
@@ -0,0 +1,97 @@
+/**********************************************************************
+* *
+* Copyright 2017 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 KRAICHNAN_FIELD_HPP
+#define KRAICHNAN_FIELD_HPP
+
+
+
+#include <cstdlib>
+#include "base.hpp"
+#include "kspace.hpp"
+#include "field.hpp"
+#include "full_code/direct_numerical_simulation.hpp"
+#include "particles/particles_system_builder.hpp"
+#include "particles/particles_output_hdf5.hpp"
+#include "particles/particles_sampling.hpp"
+
+template <typename rnumber>
+class kraichnan_field: public direct_numerical_simulation
+{
+ public:
+
+ /* parameters that are read in read_parameters */
+ double dt;
+ std::string fftw_plan_rigor;
+
+ /* particle 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;
+ int output_velocity;
+
+ /* other stuff */
+ kspace<FFTW, SMOOTH> *kk;
+ field<rnumber, FFTW, THREE> *velocity;
+ double spectrum_slope;
+ double spectrum_k_cutoff;
+ double spectrum_coefficient;
+ double max_velocity_estimate;
+ int field_random_seed;
+
+ /* other stuff */
+ std::unique_ptr<abstract_particles_system<long long int, double>> ps;
+
+ particles_output_hdf5<long long int, double,3> *particles_output_writer_mpi;
+ particles_output_sampling_hdf5<long long int, double, 3> *particles_sample_writer_mpi;
+
+
+ kraichnan_field(
+ const MPI_Comm COMMUNICATOR,
+ const std::string &simulation_name):
+ direct_numerical_simulation(
+ COMMUNICATOR,
+ simulation_name){}
+ ~kraichnan_field(){}
+
+ int initialize(void);
+ int step(void);
+ int finalize(void);
+
+ virtual int read_parameters(void);
+ int write_checkpoint(void);
+ int do_stats(void);
+
+ int generate_random_velocity(void);
+ void update_checkpoint(void);
+};
+
+#endif//KRAICHNAN_FIELD_HPP
+
diff --git a/cpp/kspace.cpp b/cpp/kspace.cpp
index 06afe4a88fbe6c1b1a05040dfbbae72d73c68d18..d601345742f6c0845e15b1b0cea2b2f1770c7a71 100644
--- a/cpp/kspace.cpp
+++ b/cpp/kspace.cpp
@@ -548,9 +548,11 @@ void kspace<be, dt>::dealias(typename fftw_interface<rnumber>::complex *__restri
template <field_backend be,
kspace_dealias_type dt>
template <typename rnumber>
-void kspace<be, dt>::force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a)
+void kspace<be, dt>::project_divfree(
+ typename fftw_interface<rnumber>::complex *__restrict__ a,
+ const bool maintain_energy)
{
- TIMEZONE("kspace::force_divfree");
+ TIMEZONE("kspace::project_divfree");
this->CLOOP_K2(
[&](ptrdiff_t cindex,
ptrdiff_t xindex,
@@ -560,18 +562,46 @@ void kspace<be, dt>::force_divfree(typename fftw_interface<rnumber>::complex *__
if (k2 > 0)
{
typename fftw_interface<rnumber>::complex tval;
+ double initial_size = 1;
+ double projected_size = 1;
tval[0] = (this->kx[xindex]*((*(a + cindex*3 ))[0]) +
this->ky[yindex]*((*(a + cindex*3+1))[0]) +
this->kz[zindex]*((*(a + cindex*3+2))[0]) ) / k2;
tval[1] = (this->kx[xindex]*((*(a + cindex*3 ))[1]) +
this->ky[yindex]*((*(a + cindex*3+1))[1]) +
this->kz[zindex]*((*(a + cindex*3+2))[1]) ) / k2;
+ if (maintain_energy)
+ {
+ initial_size = sqrt(
+ ((*(a + cindex*3 ))[0])*((*(a + cindex*3 ))[0]) +
+ ((*(a + cindex*3+1))[0])*((*(a + cindex*3+1))[0]) +
+ ((*(a + cindex*3+2))[0])*((*(a + cindex*3+2))[0]) +
+ ((*(a + cindex*3 ))[1])*((*(a + cindex*3 ))[1]) +
+ ((*(a + cindex*3+1))[1])*((*(a + cindex*3+1))[1]) +
+ ((*(a + cindex*3+2))[1])*((*(a + cindex*3+2))[1]));
+ }
for (int imag_part=0; imag_part<2; imag_part++)
{
a[cindex*3 ][imag_part] -= tval[imag_part]*this->kx[xindex];
a[cindex*3+1][imag_part] -= tval[imag_part]*this->ky[yindex];
a[cindex*3+2][imag_part] -= tval[imag_part]*this->kz[zindex];
}
+ if (maintain_energy)
+ {
+ projected_size = sqrt(
+ ((*(a + cindex*3 ))[0])*((*(a + cindex*3 ))[0]) +
+ ((*(a + cindex*3+1))[0])*((*(a + cindex*3+1))[0]) +
+ ((*(a + cindex*3+2))[0])*((*(a + cindex*3+2))[0]) +
+ ((*(a + cindex*3 ))[1])*((*(a + cindex*3 ))[1]) +
+ ((*(a + cindex*3+1))[1])*((*(a + cindex*3+1))[1]) +
+ ((*(a + cindex*3+2))[1])*((*(a + cindex*3+2))[1]));
+ if (projected_size > 0)
+ for (int component=0; component<3; component++)
+ for (int imag_part=0; imag_part<2; imag_part++)
+ {
+ (*(a + cindex*3+component))[imag_part] *= initial_size / projected_size;
+ }
+ }
}
}
);
@@ -579,6 +609,72 @@ void kspace<be, dt>::force_divfree(typename fftw_interface<rnumber>::complex *__
std::fill_n((rnumber*)(a), 6, 0.0);
}
+/** \brief Rotate vector modes perpendicular to wavenumber
+ * This is different from project because it maintains the energy of the field,
+ * I want it in order to be able to generate random fields with prescribed
+ * spectra.
+ */
+template <field_backend be,
+ kspace_dealias_type dt>
+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){
+ if (k2 > 0)
+ {
+ double cosTheta;
+ double usize;
+ for (int cc=0; cc<2; cc++)
+ {
+ // compute dot product
+ cosTheta = (this->kx[xindex]*((*(a + cindex*3 ))[cc]) +
+ this->ky[yindex]*((*(a + cindex*3+1))[cc]) +
+ this->kz[zindex]*((*(a + cindex*3+2))[cc]) );
+ // now compute size of initial velocity vector
+ usize = sqrt(pow((*(a + cindex*3 ))[cc], 2) +
+ pow((*(a + cindex*3+1))[cc], 2) +
+ pow((*(a + cindex*3+2))[cc], 2));
+ // finalize computation of cos Theta
+ if (usize != 0)
+ {
+ cosTheta /= (usize * sqrt(k2));
+ // now compute cross product
+ // cross product for complex vectors is complex conjugate of regular cross product
+ double cp[3];
+ cp[0] = (*(a + cindex*3+1))[cc]*this->kz[zindex] - (*(a + cindex*3+2))[cc]*this->ky[yindex];
+ cp[1] = (*(a + cindex*3+2))[cc]*this->kx[xindex] - (*(a + cindex*3+0))[cc]*this->kz[zindex];
+ cp[2] = (*(a + cindex*3+0))[cc]*this->ky[yindex] - (*(a + cindex*3+1))[cc]*this->kx[xindex];
+ double cpsize = sqrt(cp[0]*cp[0] + cp[1]*cp[1] + cp[2]*cp[2]);
+ cp[0] /= cpsize;
+ cp[1] /= cpsize;
+ cp[2] /= cpsize;
+ double sinTheta = sqrt(1 - cosTheta*cosTheta);
+ // we are actually interested in rotating the vector with pi/2 - Theta
+ double tmpdouble = cosTheta;
+ cosTheta = sinTheta;
+ sinTheta = -tmpdouble;
+ // store initial velocity
+ double u[3];
+ u[0] = (*(a + cindex*3+0))[cc];
+ u[1] = (*(a + cindex*3+1))[cc];
+ u[2] = (*(a + cindex*3+2))[cc];
+ // store final result
+ (*(a + cindex*3+0))[cc] = u[0]*cosTheta + (cp[1]*u[2] - cp[2]*u[1])*sinTheta;
+ (*(a + cindex*3+1))[cc] = u[1]*cosTheta + (cp[2]*u[0] - cp[0]*u[2])*sinTheta;
+ (*(a + cindex*3+2))[cc] = u[2]*cosTheta + (cp[0]*u[1] - cp[1]*u[0])*sinTheta;
+ }
+ }}
+ });
+ if (this->layout->myrank == this->layout->rank[0][0])
+ std::fill_n((rnumber*)(a), 6, 0.0);
+}
+
template <field_backend be,
kspace_dealias_type dt>
template <typename rnumber,
@@ -1068,8 +1164,20 @@ template double kspace<FFTW, SMOOTH>::L2norm<double, THREE>(
template double kspace<FFTW, SMOOTH>::L2norm<double, THREExTHREE>(
const typename fftw_interface<double>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::project_divfree<float>(
+ typename fftw_interface<float>::complex *__restrict__ a,
+ const bool);
+template void kspace<FFTW, SMOOTH>::project_divfree<double>(
+ typename fftw_interface<double>::complex *__restrict__ a,
+ const bool);
+
template void kspace<FFTW, SMOOTH>::force_divfree<float>(
typename fftw_interface<float>::complex *__restrict__ a);
template void kspace<FFTW, SMOOTH>::force_divfree<double>(
typename fftw_interface<double>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::rotate_divfree<float>(
+ typename fftw_interface<float>::complex *__restrict__ a);
+template void kspace<FFTW, SMOOTH>::rotate_divfree<double>(
+ typename fftw_interface<double>::complex *__restrict__ a);
+
diff --git a/cpp/kspace.hpp b/cpp/kspace.hpp
index ea82d54be7b9f915f5ec475c3c19ba4b10161acc..80d89f2593d7462b0dd8e28c229e3c7919a9e1ae 100644
--- a/cpp/kspace.hpp
+++ b/cpp/kspace.hpp
@@ -220,7 +220,16 @@ class kspace
}
}
template <typename rnumber>
- void force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a);
+ void project_divfree(
+ typename fftw_interface<rnumber>::complex *__restrict__ a,
+ const bool maintain_energy = false);
+ // TODO: can the following be done in a cleaner way?
+ template <typename rnumber>
+ void force_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a){
+ this->template project_divfree<rnumber>(a, false);
+ }
+ template <typename rnumber>
+ void rotate_divfree(typename fftw_interface<rnumber>::complex *__restrict__ a);
};
#endif//KSPACE_HPP