Merge branch 'feature/parameter_overwrite' into 'master'

tweaks parameter priority

See merge request !126

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(