Skip to content
Snippets Groups Projects
Commit 76581747 authored by Cristian Lalescu's avatar Cristian Lalescu
Browse files

fix old particle fields add method

parent 9474c6fc
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 159 additions and 140 deletions
bfps/NavierStokes.py
+ 159
− 140
View file @ 76581747
......@@ -117,10 +117,6 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
def write_fluid_stats(self):
self.fluid_includes += '#include <cmath>\n'
self.fluid_includes += '#include "fftw_tools.hpp"\n'
self.stat_src += """
//begincpp
//endcpp
"""
self.stat_src += """
//begincpp
double *velocity_moments = new double[10*4];
......@@ -399,7 +395,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
field_name = 'fs->rvelocity'):
name = quantity + '_' + name
self.fluid_includes += '#include "rFFTW_interpolator.hpp"\n'
self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(self.C_dtype, neighbours, name)
self.fluid_variables += 'rFFTW_interpolator <{0}, {1}> *{2};\n'.format(
self.C_dtype, neighbours, name)
self.parameters[name + '_type'] = interp_type
self.parameters[name + '_neighbours'] = neighbours
if interp_type == 'spline':
......@@ -434,7 +431,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
field_class = 'rFFTW_interpolator',
acc_field_name = 'rFFTW_acc'):
self.fluid_includes += '#include "{0}.hpp"\n'.format(field_class)
self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(self.C_dtype, neighbours, name)
self.fluid_variables += field_class + '<{0}, {1}> *vel_{2}, *acc_{2};\n'.format(
self.C_dtype, neighbours, name)
self.parameters[name + '_type'] = interp_type
self.parameters[name + '_neighbours'] = neighbours
if interp_type == 'spline':
......@@ -442,6 +440,7 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
beta_name = 'beta_n{0}_m{1}'.format(neighbours, smoothness)
elif interp_type == 'Lagrange':
beta_name = 'beta_Lagrange_n{0}'.format(neighbours)
if field_class == 'rFFTW_interpolator':
self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4}, fs->rvelocity);\n' +
'acc_{0} = new {1}<{2}, {3}>(fs, {4}, {5});\n').format(name,
field_class,
......@@ -449,6 +448,14 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
neighbours,
beta_name,
acc_field_name)
elif field_class == 'interpolator':
self.fluid_start += ('vel_{0} = new {1}<{2}, {3}>(fs, {4});\n' +
'acc_{0} = new {1}<{2}, {3}>(fs, {4});\n').format(name,
field_class,
self.C_dtype,
neighbours,
beta_name,
acc_field_name)
self.fluid_end += ('delete vel_{0};\n' +
'delete acc_{0};\n').format(name)
update_fields = 'fs->compute_velocity(fs->cvorticity);\n'
......@@ -474,9 +481,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
integration_steps = 2
neighbours = self.parameters[fields_name + '_neighbours']
self.parameters['tracers{0}_field'.format(self.particle_species)] = fields_name
self.parameters['tracers{0}_integration_method'.format(self.particle_species)] = integration_method
self.parameters['tracers{0}_integration_method'.format(
self.particle_species)] = integration_method
self.parameters['tracers{0}_kcut'.format(self.particle_species)] = kcut
self.parameters['tracers{0}_integration_steps'.format(self.particle_species)] = integration_steps
self.parameters['tracers{0}_integration_steps'.format(
self.particle_species)] = integration_steps
self.file_datasets_grow += """
//begincpp
temp_string = (std::string("/particles/") +
......@@ -552,7 +561,11 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'fname, fs, vel_{4},\n' +
'nparticles,\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
self.particle_species, self.C_dtype, multistep, neighbours, fields_name)
self.particle_species,
self.C_dtype,
multistep,
neighbours,
fields_name)
else:
self.particle_variables += 'rFFTW_particles<VELOCITY_TRACER, {0}, {1}> *ps{2};\n'.format(
self.C_dtype,
......@@ -562,7 +575,10 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
'fname, vel_{3},\n' +
'nparticles,\n' +
'niter_part, tracers{0}_integration_steps);\n').format(
self.particle_species, self.C_dtype, neighbours, fields_name)
self.particle_species,
self.C_dtype,
neighbours,
fields_name)
self.particle_start += ('ps{0}->dt = dt;\n' +
'ps{0}->iteration = iteration;\n' +
'ps{0}->read(stat_file);\n').format(self.particle_species)
......@@ -601,7 +617,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
with self.get_data_file() as data_file:
if 'moments' not in data_file['statistics'].keys():
return None
iter0 = min(data_file['statistics/moments/velocity'].shape[0]*self.parameters['niter_stat']-1,
iter0 = min((data_file['statistics/moments/velocity'].shape[0] *
self.parameters['niter_stat']-1),
iter0)
if type(iter1) == type(None):
iter1 = data_file['iteration'].value
......@@ -613,9 +630,9 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.statistics['kM'] = data_file['kspace/kM'].value
self.statistics['dk'] = data_file['kspace/dk'].value
if self.particle_species > 0:
self.trajectories = [
data_file['particles/' + key + '/state'][
iter0//self.parameters['niter_part']:iter1//self.parameters['niter_part']+1]
self.trajectories = [data_file['particles/' + key + '/state'][
iter0//self.parameters['niter_part'] :
iter1//self.parameters['niter_part']+1]
for key in data_file['particles'].keys()]
computation_needed = True
pp_file = h5py.File(self.get_postprocess_file_name(), 'a')
......@@ -660,7 +677,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.statistics[key + '(t)'] = (self.statistics['dk'] *
np.sum(self.statistics[key + '(t, k)'], axis = 1))
self.statistics['Uint(t)'] = np.sqrt(2*self.statistics['energy(t)'] / 3)
self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi / (2*self.statistics['Uint(t)']**2)) *
self.statistics['Lint(t)'] = ((self.statistics['dk']*np.pi /
(2*self.statistics['Uint(t)']**2)) *
np.nansum(self.statistics['energy(t, k)'] /
self.statistics['kshell'][None, :], axis = 1))
for key in ['energy',
......@@ -687,7 +705,8 @@ class NavierStokes(bfps.fluid_base.fluid_particle_base):
self.statistics['Rlambda' + suffix] = (self.statistics['Uint' + suffix] *
self.statistics['lambda' + suffix] /
self.parameters['nu'])
self.statistics['kMeta' + suffix] = self.statistics['kM']*self.statistics['etaK' + suffix]
self.statistics['kMeta' + suffix] = (self.statistics['kM'] *
self.statistics['etaK' + suffix])
if self.parameters['dealias_type'] == 1:
self.statistics['kMeta' + suffix] *= 0.8
self.statistics['Tint'] = self.statistics['Lint'] / self.statistics['Uint']
......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment