diff --git a/TurTLE/DNS.py b/TurTLE/DNS.py
index b04c8b360307632e759eddbe81a03d555473ba96..606014b09a87b490f1dda781fc6bc6e501d12eb8 100644
--- a/TurTLE/DNS.py
+++ b/TurTLE/DNS.py
@@ -93,7 +93,7 @@ class DNS(_code):
return None
def generate_default_parameters(self):
# these parameters are relevant for all DNS classes
- self.parameters['fftw_plan_rigor'] = 'FFTW_ESTIMATE'
+ self.parameters['fftw_plan_rigor'] = str('FFTW_ESTIMATE')
self.parameter_description['fftw_plan_rigor'] = 'FFTW plan rigor to use. One of `FFTW_ESTIMATE`, `FFTW_MEASURE`, `FFTW_PATIENT`. Please see FFTW documentation.'
self.parameters['dealias_type'] = int(1)
self.parameter_description['dealias_type'] = 'Dealiasing mehtod to use, integer. Options are: two-thirds (0) or smooth (1).'
@@ -125,7 +125,7 @@ class DNS(_code):
self.parameter_description['fk0'] = 'Forcing parameter: if forcing acts on wavenumber band, this is the smallest wavenumber where it acts (in code units).'
self.parameters['fk1'] = float(4.0)
self.parameter_description['fk1'] = 'Forcing parameter: if forcing acts on wavenumber band, this is the largest wavenumber where it acts (in code units).'
- self.parameters['forcing_type'] = 'fixed_energy_injection_rate'
+ self.parameters['forcing_type'] = str('fixed_energy_injection_rate')
self.parameter_description['forcing_type'] = 'Forcing parameter: what type of force to use.'
self.parameters['histogram_bins'] = int(256)
self.parameter_description['histogram_bins'] = 'During statistics, histograms of real-valued fields are computed using a number of `HISTOGRAM_BINS` bins.'
@@ -324,22 +324,39 @@ class DNS(_code):
'--src-wd',
type = str,
dest = 'src_work_dir',
+ help = ('Empty by default.'
+ + ' Location of simulation from which to source initial conditions.'
+ + ' See also `--src-simname`.'),
default = '')
parser.add_argument(
'--src-simname',
type = str,
dest = 'src_simname',
+ help = ('Empty by default.'
+ + ' Name of simulation from which to source initial conditions.'
+ + ' If the parameter file of the source simulation is present, '
+ + 'the new DNS will use the same physical parameters '
+ + '(forcing parameters and viscosity) and the same dealiasing option.'
+ + ' Any other parameters (e.g. "nx") have to be specified separately.'
+ + ' See also `--overwrite-src-parameters`.'),
default = '')
parser.add_argument(
'--src-iteration',
type = int,
dest = 'src_iteration',
+ help = ('0 by default.'
+ + ' Iteration of source simulation to be used as initial condition.'
+ + ' See also `--src-simname`.'),
default = 0)
parser.add_argument(
'--overwrite-src-parameters',
action = 'store_true',
dest = 'overwrite_src_parameters',
- help = 'False by default. If using a source simulation, do we keep its physical parameters? Note: if this option is switched on, default values are used for parameters that are not set explicitly by the command --- NOT the values of the source simulation.')
+ help = ('False by default.'
+ + ' If using a source simulation, do we keep its physical parameters?'
+ + ' Note: if this option is switched on, default values are used for '
+ + 'parameters that are not set explicitly by the command --- '
+ + 'NOT the values of the source simulation.'))
parser.add_argument(
'--kMeta',
type = float,
@@ -472,7 +489,23 @@ class DNS(_code):
self,
args = [],
extra_parameters = None):
- """Set up reasonable parameters.
+ """Assign DNS parameters.
+
+ There are three ways to specify parameters:
+
+ * Directly via the command line.
+ * If the simulation takes its initial condition from another simulation
+ (via the :code:`--src-dir` and :code:`--src-simname` parameters) and
+ some physical parameter is not provided on the command line, then this
+ method will attempt to read that parameter from the parameter file of
+ the source simulation.
+ * If a parameter is not specified via the command line, and a source
+ simulation is not provided, or the source parameter file is missing, or
+ the :code:`--overwrite-src-parameters` option is turned on, this method
+ will generate a "reasonable" value for that parameter.
+ "Reasonable" refers to either some fixed values for the forcing
+ parameters, or an empirical extrapolation of the viscosity to satisfy a
+ given value of the :code:`--kMeta` command line parameter.
With the default Lundgren forcing applied in the band [2, 4],
we can estimate the dissipation, therefore we can estimate
@@ -520,10 +553,32 @@ class DNS(_code):
for k in self.extra_parameters[self.dns_type].keys():
self.parameters[k] = self.extra_parameters[self.dns_type][k]
+ # for parameters there are three options, in order of preference:
+ # * value given from command line
+ # * value used by source simulation
+ # * default value
+
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:
- opt.src_work_dir = os.path.realpath(opt.work_dir)
+
+ if len(opt.src_simname) > 0:
+ if len(opt.src_work_dir) == 0:
+ opt.src_work_dir = os.path.realpath(opt.work_dir)
+ src_file_name = os.path.join(
+ os.path.realpath(opt.src_work_dir),
+ opt.src_simname + '.h5')
+ if os.path.exists(src_file_name):
+ # open source file for use later
+ src_file = h5py.File(src_file_name, 'r')
+ else:
+ print('TurTLE WARNING: unable to copy parameters of source simulation.'
+ + 'File {0} is missing.\n'.format(src_file_name)
+ + 'DNS will use default values for physical parameters.\n')
+ src_file = None
+ else:
+ src_file = None
+ if opt.overwrite_src_parameters:
+ src_file = None
if type(opt.dkx) == type(None):
opt.dkx = 2. / opt.Lx
if type(opt.dky) == type(None):
@@ -536,16 +591,60 @@ class DNS(_code):
opt.ny = opt.n
if type(opt.nz) == type(None):
opt.nz = opt.n
+ if type(opt.friction_coefficient) == type(None):
+ if type(src_file) != type(None):
+ self.parameters['friction_coefficient'] = type(self.parameters['friction_coefficient'])(
+ src_file['parameters/friction_coefficient'][()])
+ else:
+ opt.friction_coefficient = self.parameters['friction_coefficient']
+ if type(opt.fmode) == type(None):
+ if type(src_file) != type(None):
+ self.parameters['fmode'] = type(self.parameters['fmode'])(
+ src_file['parameters/fmode'][()])
+ else:
+ opt.fmode = self.parameters['fmode']
+ if type(opt.famplitude) == type(None):
+ if type(src_file) != type(None):
+ self.parameters['famplitude'] = type(self.parameters['famplitude'])(
+ src_file['parameters/famplitude'][()])
+ else:
+ opt.famplitude = self.parameters['famplitude']
+ if type(opt.energy) == type(None):
+ if type(src_file) != type(None):
+ self.parameters['energy'] = type(self.parameters['energy'])(
+ src_file['parameters/energy'][()])
+ else:
+ opt.energy = self.parameters['energy']
if type(opt.fk0) == type(None):
- opt.fk0 = self.parameters['fk0']
+ if type(src_file) != type(None):
+ self.parameters['fk0'] = type(self.parameters['fk0'])(
+ src_file['parameters/fk0'][()])
+ else:
+ opt.fk0 = self.parameters['fk0']
if type(opt.fk1) == type(None):
- opt.fk1 = self.parameters['fk1']
+ if type(src_file) != type(None):
+ self.parameters['fk1'] = type(self.parameters['fk1'])(
+ src_file['parameters/fk1'][()])
+ else:
+ opt.fk1 = self.parameters['fk1']
if type(opt.forcing_type) == type(None):
- opt.forcing_type = self.parameters['forcing_type']
+ if type(src_file) != type(None):
+ self.parameters['forcing_type'] = str(
+ src_file['parameters/forcing_type'][()].decode())
+ else:
+ opt.forcing_type = self.parameters['forcing_type']
if type(opt.injection_rate) == type(None):
- opt.injection_rate = self.parameters['injection_rate']
+ if type(src_file) != type(None):
+ self.parameters['injection_rate'] = type(self.parameters['injection_rate'])(
+ src_file['parameters/injection_rate'][()])
+ else:
+ opt.injection_rate = self.parameters['injection_rate']
if type(opt.dealias_type) == type(None):
- opt.dealias_type = self.parameters['dealias_type']
+ if type(src_file) != type(None):
+ self.parameters['dealias_type'] = type(self.parameters['dealias_type'])(
+ src_file['parameters/dealias_type'][()])
+ else:
+ opt.dealias_type = self.parameters['dealias_type']
if (opt.nx > opt.n or
opt.ny > opt.n or
opt.nz > opt.n):
@@ -561,34 +660,33 @@ class DNS(_code):
kM = opt.n * 0.5
if opt.dealias_type == 1:
kM *= 0.8
- # tweak forcing/viscosity based on forcing type
- if opt.forcing_type == 'linear':
- # custom famplitude for 288 and 576
- if opt.n == 288:
- self.parameters['famplitude'] = 0.45
- elif opt.n == 576:
- self.parameters['famplitude'] = 0.47
- elif opt.forcing_type == 'fixed_energy_injection_rate':
- # use the fact that mean dissipation rate is equal to injection rate
- if type(opt.nu) == type(None):
- opt.nu = (
- opt.injection_rate *
- (opt.kMeta / kM)**4)**(1./3)
- if opt.injection_rate == 0:
- opt.nu = (
- (opt.kMeta / kM)**(4./3) *
- (np.pi / 1.5)**(1./3) *
- (2*self.parameters['energy'] / 3)**0.5)
- elif opt.forcing_type == 'fixed_energy':
- kf = 1. / (1./opt.fk0 +
- 1./opt.fk1)
- if type(opt.nu) == type(None):
- opt.nu = (
- (opt.kMeta / kM)**(4./3) *
- (np.pi / kf)**(1./3) *
- (2*self.parameters['energy'] / 3)**0.5)
+
if type(opt.nu) == type(None):
- opt.nu = (opt.kMeta * 2 / opt.n)**(4./3)
+ # check if we are using viscosity of source simulation
+ if type(src_file) != type(None):
+ opt.nu = type(self.parameters['nu'])(
+ src_file['parameters/nu'][()])
+ # empirical viscosity based on forcing type
+ else:
+ if opt.forcing_type == 'fixed_energy_injection_rate':
+ # use the fact that mean dissipation rate is equal to injection rate
+ opt.nu = (
+ opt.injection_rate *
+ (opt.kMeta / kM)**4)**(1./3)
+ if opt.injection_rate == 0:
+ opt.nu = (
+ (opt.kMeta / kM)**(4./3) *
+ (np.pi / 1.5)**(1./3) *
+ (2*self.parameters['energy'] / 3)**0.5)
+ elif opt.forcing_type == 'fixed_energy':
+ kf = 1. / (1./opt.fk0 +
+ 1./opt.fk1)
+ opt.nu = (
+ (opt.kMeta / kM)**(4./3) *
+ (np.pi / kf)**(1./3) *
+ (2*self.parameters['energy'] / 3)**0.5)
+ else:
+ opt.nu = (opt.kMeta * 2 / opt.n)**(4./3)
self.parameters['nu'] = opt.nu
if type(opt.checkpoints_per_file) == type(None):
# hardcoded FFTW complex representation size
@@ -898,35 +996,6 @@ class DNS(_code):
checkpoint_field + '/complex/{0}'.format(opt.src_iteration),
f,
checkpoint_field + '/complex/{0}'.format(0))
- if not opt.overwrite_src_parameters:
- print('DNS will use physical parameters from source simulation.\n'
- 'This means `dealias_type`, `energy`, `famplitude`, `fk0`, `fk1`, `fmode`, `friction_coefficient`, `injection_rate`, `nu` and `forcing_type`.')
- # timestep left out on purpose
- src_file = os.path.join(
- os.path.realpath(opt.src_work_dir),
- opt.src_simname + '.h5')
- if os.path.exists(src_file):
- f0 = h5py.File(src_file, 'r')
- for kk in ['dealias_type',
- 'energy',
- 'famplitude',
- 'fk0',
- 'fk1',
- 'fmode',
- 'friction_coefficient',
- 'injection_rate',
- 'nu']:
- if (kk in f0['parameters'].keys()) and (kk in self.parameters.keys()):
- self.parameters[kk] = type(self.parameters[kk])(f0['parameters/' + kk][()])
- self.parameters['forcing_type'] = str(f0['parameters/forcing_type'][()], 'ascii')
- f0.close()
- else:
- print('TurTLE WARNING: unable to copy parameters of source simulation.'
- + 'File {0} is missing.\n'.format(src_file)
- + 'DNS will use default values for physical parameters.\n')
- else:
- print('TurTLE WARNING: DNS will use default values for physical parameters.\n'
- 'This means `dealias_type`, `energy`, `famplitude`, `fk0`, `fk1`, `fmode`, `friction_coefficient`, `injection_rate`, `nu` and `forcing_type` may differ from values used by the source DNS.')
else:
if self.dns_type == 'NSVE_Stokes_particles':
data = self.generate_vector_field(