Commit d4f77830 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

Merge branch 'feature/Stokes-drag' into 'develop'

Feature/stokes drag

See merge request clalescu/turtle!26
parents 6ddcee26 3ad8cf2d
Pipeline #68268 passed with stage
in 7 minutes
...@@ -3,20 +3,20 @@ ...@@ -3,20 +3,20 @@
# Copyright 2019 Max Planck Institute # # Copyright 2019 Max Planck Institute #
# for Dynamics and Self-Organization # # 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 # # it under the terms of the GNU General Public License as published #
# by the Free Software Foundation, either version 3 of the License, # # by the Free Software Foundation, either version 3 of the License, #
# or (at your option) any later version. # # 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 # # but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. # # GNU General Public License for more details. #
# # # #
# You should have received a copy of the GNU General Public License # # 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 # # Contact: Cristian.Lalescu@ds.mpg.de #
# # # #
...@@ -300,6 +300,7 @@ set(cpp_for_lib ...@@ -300,6 +300,7 @@ set(cpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.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/full_code/NSVEp_extra_sampling.cpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp ${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.cpp
${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.cpp ${PROJECT_SOURCE_DIR}/cpp/full_code/ornstein_uhlenbeck_process.cpp
...@@ -348,8 +349,11 @@ set(hpp_for_lib ...@@ -348,8 +349,11 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.hpp ${PROJECT_SOURCE_DIR}/cpp/full_code/test_interpolation.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.hpp ${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEparticles.hpp
${PROJECT_SOURCE_DIR}/cpp/full_code/NSVEcomplex_particles.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/full_code/NSVEp_extra_sampling.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_inner_computer.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_input.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/abstract_particles_output.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.hpp
...@@ -367,7 +371,6 @@ set(hpp_for_lib ...@@ -367,7 +371,6 @@ set(hpp_for_lib
${PROJECT_SOURCE_DIR}/cpp/particles/particles_distr_mpi.hpp ${PROJECT_SOURCE_DIR}/cpp/particles/particles_distr_mpi.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_field_computer.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_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_input_hdf5.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_output_hdf5.hpp ${PROJECT_SOURCE_DIR}/cpp/particles/particles_output_hdf5.hpp
${PROJECT_SOURCE_DIR}/cpp/particles/particles_output_mpiio.hpp ${PROJECT_SOURCE_DIR}/cpp/particles/particles_output_mpiio.hpp
...@@ -419,3 +422,66 @@ else() ...@@ -419,3 +422,66 @@ 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)") 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() endif()
install(CODE "execute_process(COMMAND python3 ${PROJECT_SOURCE_DIR}/setup.py install --force --prefix=${CMAKE_INSTALL_PREFIX} WORKING_DIRECTORY ${PROJECT_BINARY_DIR}/python/)") 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
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_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 --nparticles 10000
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)
...@@ -439,7 +439,13 @@ class DNS(_code): ...@@ -439,7 +439,13 @@ class DNS(_code):
assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0) assert (self.parameters['niter_todo'] % self.parameters['niter_stat'] == 0)
assert (self.parameters['niter_todo'] % self.parameters['niter_out'] == 0) assert (self.parameters['niter_todo'] % self.parameters['niter_out'] == 0)
assert (self.parameters['niter_out'] % self.parameters['niter_stat'] == 0) assert (self.parameters['niter_out'] % self.parameters['niter_stat'] == 0)
if self.dns_type in ['NSVEparticles_no_output', 'NSVEcomplex_particles', 'NSVEparticles', 'static_field', 'kraichnan_field']: if self.dns_type in [
'NSVEparticles_no_output',
'NSVEcomplex_particles',
'NSVE_Stokes_particles',
'NSVEparticles',
'static_field',
'kraichnan_field']:
assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0) assert (self.parameters['niter_todo'] % self.parameters['niter_part'] == 0)
assert (self.parameters['niter_out'] % self.parameters['niter_part'] == 0) assert (self.parameters['niter_out'] % self.parameters['niter_part'] == 0)
_code.write_par(self, iter0 = iter0) _code.write_par(self, iter0 = iter0)
...@@ -648,6 +654,10 @@ class DNS(_code): ...@@ -648,6 +654,10 @@ class DNS(_code):
'NSVEparticles', 'NSVEparticles',
help = 'plain Navier-Stokes vorticity formulation, with basic fluid tracers') 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( parser_NSVEp2p = subparsers.add_parser(
'NSVEcomplex_particles', 'NSVEcomplex_particles',
help = 'plain Navier-Stokes vorticity formulation, with oriented active particles') help = 'plain Navier-Stokes vorticity formulation, with oriented active particles')
...@@ -660,6 +670,7 @@ class DNS(_code): ...@@ -660,6 +670,7 @@ class DNS(_code):
'NSVEparticles_no_output', 'NSVEparticles_no_output',
'NSVEp2', 'NSVEp2',
'NSVEp2p', 'NSVEp2p',
'NSVE_Stokes_particles',
'NSVEp_extra', 'NSVEp_extra',
'static_field', 'static_field',
'kraichnan_field']: 'kraichnan_field']:
...@@ -673,6 +684,7 @@ class DNS(_code): ...@@ -673,6 +684,7 @@ class DNS(_code):
'NSVEparticles_no_output', 'NSVEparticles_no_output',
'NSVEp2', 'NSVEp2',
'NSVEp2p', 'NSVEp2p',
'NSVE_Stokes_particles',
'NSVEp_extra', 'NSVEp_extra',
'static_field', 'static_field',
'kraichnan_field']: 'kraichnan_field']:
...@@ -691,6 +703,10 @@ class DNS(_code): ...@@ -691,6 +703,10 @@ class DNS(_code):
pars['spectrum_slope'] = float(-5./3) pars['spectrum_slope'] = float(-5./3)
pars['spectrum_k_cutoff'] = float(16) pars['spectrum_k_cutoff'] = float(16)
pars['spectrum_coefficient'] = float(1) pars['spectrum_coefficient'] = float(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 return pars
def prepare_launch( def prepare_launch(
self, self,
...@@ -726,7 +742,14 @@ class DNS(_code): ...@@ -726,7 +742,14 @@ class DNS(_code):
self.dns_type = opt.DNS_class self.dns_type = opt.DNS_class
self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__ self.name = self.dns_type + '-' + self.fluid_precision + '-v' + TurTLE.__version__
# merge parameters if needed # merge parameters if needed
if self.dns_type in ['NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling', 'static_field', 'kraichnan_field']: if self.dns_type in [
'NSVEparticles',
'NSVEcomplex_particles',
'NSVE_Stokes_particles',
'NSVEparticles_no_output',
'NSVEp_extra_sampling',
'static_field',
'kraichnan_field']:
for k in self.NSVEp_extra_parameters.keys(): for k in self.NSVEp_extra_parameters.keys():
self.parameters[k] = self.NSVEp_extra_parameters[k] self.parameters[k] = self.NSVEp_extra_parameters[k]
if type(extra_parameters) != type(None): if type(extra_parameters) != type(None):
...@@ -797,7 +820,14 @@ class DNS(_code): ...@@ -797,7 +820,14 @@ class DNS(_code):
# hardcoded FFTW complex representation size # hardcoded FFTW complex representation size
field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize field_size = 3*(opt.nx+2)*opt.ny*opt.nz*self.fluid_dtype.itemsize
checkpoint_size = field_size checkpoint_size = field_size
if self.dns_type in ['kraichnan_field', 'static_field', 'NSVEparticles', 'NSVEcomplex_particles', 'NSVEparticles_no_output', 'NSVEp_extra_sampling']: if self.dns_type in [
'kraichnan_field',
'static_field',
'NSVEparticles',
'NSVEcomplex_particles',
'NSVE_Stokes_particles',
'NSVEparticles_no_output',
'NSVEp_extra_sampling']:
rhs_size = self.parameters['tracers0_integration_steps'] rhs_size = self.parameters['tracers0_integration_steps']
if type(opt.tracers0_integration_steps) != type(None): if type(opt.tracers0_integration_steps) != type(None):
rhs_size = opt.tracers0_integration_steps rhs_size = opt.tracers0_integration_steps
...@@ -837,7 +867,9 @@ class DNS(_code): ...@@ -837,7 +867,9 @@ class DNS(_code):
integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps'] integration_steps = self.NSVEp_extra_parameters['tracers0_integration_steps']
if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys(): if 'tracers{0}_integration_steps'.format(species) in self.parameters.keys():
integration_steps = self.parameters['tracers{0}_integration_steps'.format(species)] 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 ncomponents = 6
with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file: with h5py.File(self.get_checkpoint_0_fname(), 'a') as data_file:
nn = self.parameters['nparticles'] nn = self.parameters['nparticles']
...@@ -869,12 +901,18 @@ class DNS(_code): ...@@ -869,12 +901,18 @@ class DNS(_code):
if nn > batch_size: if nn > batch_size:
dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size) dset[cc*batch_size:(cc+1)*batch_size, :3] = get_random_phases(batch_size)
if dset.shape[1] == 6: 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 nn -= batch_size
else: else:
dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn) dset[cc*batch_size:cc*batch_size+nn, :3] = get_random_phases(nn)
if dset.shape[1] == 6: 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 nn = 0
cc += 1 cc += 1
except Exception as e: except Exception as e:
...@@ -1021,6 +1059,8 @@ class DNS(_code): ...@@ -1021,6 +1059,8 @@ class DNS(_code):
if self.dns_type in ['NSVEcomplex_particles']: if self.dns_type in ['NSVEcomplex_particles']:
particle_file.create_group('tracers0/orientation') particle_file.create_group('tracers0/orientation')
particle_file.create_group('tracers0/velocity_gradient') 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']: if self.dns_type in ['NSVEp_extra_sampling']:
particle_file.create_group('tracers0/velocity_gradient') particle_file.create_group('tracers0/velocity_gradient')
particle_file.create_group('tracers0/pressure') particle_file.create_group('tracers0/pressure')
...@@ -1041,6 +1081,7 @@ class DNS(_code): ...@@ -1041,6 +1081,7 @@ class DNS(_code):
'static_field', 'static_field',
'NSVEparticles', 'NSVEparticles',
'NSVEcomplex_particles', 'NSVEcomplex_particles',
'NSVE_Stokes_particles',
'NSVEparticles_no_output', 'NSVEparticles_no_output',
'NSVEp_extra_sampling']: 'NSVEp_extra_sampling']:
if not os.path.exists(self.get_checkpoint_0_fname()): if not os.path.exists(self.get_checkpoint_0_fname()):
...@@ -1076,10 +1117,16 @@ class DNS(_code): ...@@ -1076,10 +1117,16 @@ class DNS(_code):
f, f,
'vorticity/complex/{0}'.format(0)) 'vorticity/complex/{0}'.format(0))
else: 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, write_to_file = False,
spectra_slope = 2.0, spectra_slope = 2.0,
amplitude = 0.05) amplitude = 0.05)
f['vorticity/complex/{0}'.format(0)] = data f['vorticity/complex/{0}'.format(0)] = data
f.close() f.close()
if self.dns_type == 'kraichnan_field': if self.dns_type == 'kraichnan_field':
...@@ -1093,6 +1140,7 @@ class DNS(_code): ...@@ -1093,6 +1140,7 @@ class DNS(_code):
'static_field', 'static_field',
'NSVEparticles', 'NSVEparticles',
'NSVEcomplex_particles', 'NSVEcomplex_particles',
'NSVE_Stokes_particles',
'NSVEparticles_no_output', 'NSVEparticles_no_output',
'NSVEp_extra_sampling']: 'NSVEp_extra_sampling']:
self.generate_particle_data(opt = opt) self.generate_particle_data(opt = opt)
......
#! /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 main():
"""
Run test where the trajectory of initially moving particles in a
quiescent flow is compared to the analytical exponential
"""
drag_coefficients = ['30.','20.','10.','1.','0.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',
'--checkpoints_per_file', '{0}'.format(3),
'--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)
def theory_exponential(t):
return f2['parameters/initial_particle_vel']*np.exp(-t*np.float(drag_coeff))
numerics, = plt.plot(time_values, particle_momentum, alpha = 0.3, lw=2)
plt.plot(time_values, theory_exponential(time_values), color=numerics.get_color(), ls = '--')
plt.plot(0,0,color=numerics.get_color(),label=r'$\frac{1}{St}=$'+drag_coeff)
plt.legend(loc='lower left')
print(f['tracers0/momentum/0'][:])
print(f['tracers0/momentum/100'][:]-f['tracers0/momentum/0'][:])
print(f['tracers0/position/100'][:]-f['tracers0/position/0'][:])
plt.savefig('Stokes_test1.pdf')
j+=1
return None
if __name__ == '__main__':
main()
...@@ -4,20 +4,20 @@ ...@@ -4,20 +4,20 @@
# Copyright 2019 Max Planck Institute # # Copyright 2019 Max Planck Institute #
# for Dynamics and Self-Organization # # 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 # # it under the terms of the GNU General Public License as published #
# by the Free Software Foundation, either version 3 of the License, # # by the Free Software Foundation, either version 3 of the License, #
# or (at your option) any later version. # # 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 # # but WITHOUT ANY WARRANTY; without even the implied warranty of #
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the #
# GNU General Public License for more details. # # GNU General Public License for more details. #
# # # #
# You should have received a copy of the GNU General Public License # # 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 # # Contact: Cristian.Lalescu@ds.mpg.de #
# # # #
......
...@@ -799,8 +799,8 @@ int field<rnumber, be, fc>::write_filtered( ...@@ -799,8 +799,8 @@ int field<rnumber, be, fc>::write_filtered(
} }
} }
} }
DEBUG_MSG("count[0] = %ld, offset[0] = %ld\n", //DEBUG_MSG("count[0] = %ld, offset[0] = %ld\n",
count[0], offset[0]); // count[0], offset[0]);
if (nz>=2) if (nz>=2)
{ {
assert(nz%2==0); assert(nz%2==0);
...@@ -812,8 +812,8 @@ int field<rnumber, be, fc>::write_filtered( ...@@ -812,8 +812,8 @@ int field<rnumber, be, fc>::write_filtered(
offset[1] = cz*nz/2; offset[1] = cz*nz/2;
memshape [1] = this->clayout->sizes[1]; memshape [1] = this->clayout->sizes[1];
memoffset[1] = cz*(this->clayout->sizes[1] - nz/2); memoffset[1] = cz*(this->clayout->sizes[1] - nz/2);
DEBUG_MSG("cz = %d, count[1] + offset[1] = %ld\n", //DEBUG_MSG("cz = %d, count[1] + offset[1] = %ld\n",
cz, count[1] + offset[1]); // cz, count[1] + offset[1]);
//now write data //now write data
mspace = H5Screate_simple(ndim(fc), memshape, NULL); mspace = H5Screate_simple(ndim(fc), memshape, NULL);
...@@ -830,7 +830,7 @@ int field<rnumber, be, fc>::write_filtered( ...@@ -830,7 +830,7 @@ int field<rnumber, be, fc>::write_filtered(
offset[1] = 0; offset[1] = 0;
memshape [1] = this->clayout->sizes[1]; memshape [1] = this->clayout->sizes[1];
memoffset[1] = 0; memoffset[1] = 0;
DEBUG_MSG("filtered write nz=1\n"); //DEBUG_MSG("filtered write nz=1\n");
//now write data //now write data
mspace = H5Screate_simple(ndim(fc), memshape, NULL); mspace = H5Screate_simple(ndim(fc), memshape, NULL);
...@@ -921,9 +921,9 @@ void field<rnumber, be, fc>::compute_rspace_stats( ...@@ -921,9 +921,9 @@ void field<rnumber, be, fc>::compute_rspace_stats(
// and possible component indices that are NOT averaged over // and possible component indices that are NOT averaged over
// HDF5 dataset has 2 extra indices for time and order of moment, // HDF5 dataset has 2 extra indices for time and order of moment,
// and possible component indices // and possible component indices
DEBUG_MSG("ndims = %d, ndim(fc) = %d, dsetname = %s\n", //DEBUG_MSG("ndims = %d, ndim(fc) = %d, dsetname = %s\n",
ndims, ndim(fc), // ndims, ndim(fc),
dset_name.c_str()); // dset_name.c_str());
assert(ndims == int(ndim(fc))-1); assert(ndims == int(ndim(fc))-1);
assert(dims[1] == nmoments); assert(dims[1] == nmoments);
switch(ndims) switch(ndims)
...@@ -2130,12 +2130,12 @@ int make_gaussian_random_field( ...@@ -2130,12 +2130,12 @@ int make_gaussian_random_field(
rseed*omp_get_max_threads()*output_field->clayout->nprocs + rseed*omp_get_max_threads()*output_field->clayout->nprocs +
output_field->clayout->myrank*omp_get_max_threads() + output_field->clayout->myrank*omp_get_max_threads() +
thread_id); 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); rgen[thread_id].seed(current_seed);
} }
output_field->real_space_representation = false; output_field->real_space_representation = false;
*output_field = 0.0; *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 // inside loop use only thread-local random number generator
kk->CLOOP_K2([&]( kk->CLOOP_K2([&](
ptrdiff_t cindex, ptrdiff_t cindex,
......
/**********************************************************************
* *
* 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 *