Update templates (#45)

This PR adds new variables which can be configured in the templates using the `CaseStudy` class. These are: * `bed_roughness`: set the uniform bed roughness as a manning number * `turbine_turbulence_model`: choose between the default turbulence model for the turbine wake or a canopy model. * `beta_p`: set the turbine turbulence canopy model "production" coefficient for the canopy model * `beta_d`: set the turbine turbulence canopy model "dissipation" coefficient for the canopy model * `c_epp`: set the turbine turbulence canopy model "production" closure coefficient for the canopy model * `c_epd`: set the turbine turbulence canopy model "dissipation" closure coefficient for the canopy model Additionally, the model files created in the `basic.py` example can now be kept for inspection by passing the `--persistent` flag as a command line argument. Bump to version 0.8.2
Data-Only-Greater · Jul 4, 2022 · 974a53d · 974a53d
1 parent f242804
commit 974a53d
Show file tree

Hide file tree

Showing 13 changed files with 184 additions and 52 deletions.
diff --git a/.conda/recipe/meta.yaml b/.conda/recipe/meta.yaml
@@ -1,4 +1,4 @@
-{% set version = "0.8.1" %}
+{% set version = "0.8.2" %}
 
 package:
   name: snl-delft3d-cec-verify

diff --git a/README.md b/README.md
@@ -33,7 +33,7 @@ From a conda prompt create a named environment in which to install the
 for future updates:
 
 ```
-(base) > conda create -y -n snld3d --override-channels -c conda-forge -c dataonlygreater snl-delft3d-cec-verify=0.8.1
+(base) > conda create -y -n snld3d --override-channels -c conda-forge -c dataonlygreater snl-delft3d-cec-verify=0.8.2
 (base) > conda activate snld3d
 (snld3d) > conda config --env --add channels conda-forge --add channels dataonlygreater
 (snld3d) > conda config --env --set channel_priority strict
@@ -159,6 +159,14 @@ If successful, the report files (and images) will be placed into a
 sub-directory based on the model type. For the flexible mesh model, this is 
 `fm/basic_report`.
 
+By default, the temporary directories in which the models are run are deleted
+upon completion. To keep these directories, use the `--persistent` flag. For 
+example:
+
+```
+(snld3d) > python basic.py structured --persistent
+```
+
 ### Validation Example
 
 Required files:

diff --git a/docs/conf.py b/docs/conf.py
@@ -25,7 +25,7 @@
 author = 'Mathew Topper'
 
 # The full version, including alpha/beta/rc tags
-release = '0.8.1'
+release = '0.8.2'
 
 
 # -- General configuration ---------------------------------------------------
@@ -57,7 +57,7 @@
 smv_remote_whitelist = r'^(origin)$'
 smv_tag_whitelist = r'^v(\d+\.\d+\.\d+)$' # r'^v(?!0.4.0|0.4.1|0.4.2)\d+\.\d+\.\d+$'
 smv_released_pattern = r'^refs/tags/.*$'
-smv_latest_version = 'v0.8.1'
+smv_latest_version = 'v0.8.2'
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ['_templates']

diff --git a/examples/basic.py b/examples/basic.py
@@ -4,30 +4,16 @@
 import platform
 import tempfile
 from pathlib import Path
+from contextlib import contextmanager
 from collections import defaultdict
 
 import pandas as pd
 import matplotlib.pyplot as plt
 
 from snl_d3d_cec_verify import CaseStudy, Report, Result, Runner, Template
 
-def get_d3d_bin_path():
-
-    env = dict(os.environ)
-
-    if 'D3D_BIN' in env:
-        root = Path(env['D3D_BIN'].replace('"', ''))
-        print('D3D_BIN found')
-    else:
-        root = Path("..") / "src" / "bin"
-        print('D3D_BIN not found')
-
-    print(f'Setting bin folder path to {root.resolve()}')
-
-    return root.resolve()
 
-
-def main(template_type):
+def main(template_type, persistent=False):
 
     template = Template(template_type)
     runner = Runner(get_d3d_bin_path())
@@ -40,7 +26,7 @@ def main(template_type):
 
     for i, case in enumerate(cases):
 
-        with tempfile.TemporaryDirectory() as tmpdirname:
+        with mkdtemp_persistent(persistent=persistent) as tmpdirname:
 
             # Create the project and then run it
             template(case, tmpdirname)
@@ -84,6 +70,9 @@ def main(template_type):
                                      "Turbulent kinetic energy, $k$ "
                                      "(m^2^/s^2^)",
                                      width="5.7in")
+
+            if persistent:
+                print(f"Template path: {tmpdirname}")
 
     df = pd.DataFrame(data)
 
@@ -118,6 +107,32 @@ def main(template_type):
         print(report)
 
 
+def get_d3d_bin_path():
+
+    env = dict(os.environ)
+
+    if 'D3D_BIN' in env:
+        root = Path(env['D3D_BIN'].replace('"', ''))
+        print('D3D_BIN found')
+    else:
+        root = Path("..") / "src" / "bin"
+        print('D3D_BIN not found')
+
+    print(f'Setting bin folder path to {root.resolve()}')
+
+    return root.resolve()
+
+
+def mkdtemp_persistent(*args, persistent=True, **kwargs):
+    if persistent:
+        @contextmanager
+        def normal_mkdtemp():
+            yield tempfile.mkdtemp(*args, **kwargs)
+        return normal_mkdtemp()
+    else:
+        return tempfile.TemporaryDirectory(*args, **kwargs)
+
+
 if __name__ == "__main__":
 
     import argparse
@@ -127,6 +142,9 @@ def main(template_type):
                         choices=['fm', 'structured'],
                         help=("the type of model to be exectuted - choose "
                               "'fm' or 'structured'"))
+    parser.add_argument('--persistent',
+                        action='store_true',
+                        help=("make model directories persistent"))
 
     args = parser.parse_args()
-    main(args.MODEL)
+    main(args.MODEL, args.persistent)
diff --git a/setup.cfg b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = SNL-Delft3D-CEC-Verify
-version = 0.8.1
+version = 0.8.2
 author = Mathew Topper
 author_email = [email protected]
 description = Automated verification of SNL-Delft3D-CEC based on the 2014 Mycek experiment

diff --git a/src/snl_d3d_cec_verify/cases.py b/src/snl_d3d_cec_verify/cases.py
@@ -65,6 +65,8 @@ class CaseStudy:
     :param turb_pos_z: turbine z-position, in meters. Defaults to {turb_pos_z}
     :param discharge: inlet boundary discharge, in cubic meters per second.
         Defaults to {discharge}
+    :param bed_roughness: uniform bed roughness coefficient, as a Manning 
+        number, in seconds over metres cube rooted. Defaults to {bed_roughness}
     :param horizontal_eddy_viscosity: uniform horizontal eddy viscosity, in
         metres squared per second. Defaults to {horizontal_eddy_viscosity}
     :param horizontal_eddy_diffusivity: uniform horizontal eddy diffusivity,
@@ -75,6 +77,21 @@ class CaseStudy:
         in metres squared per second. Defaults to {vertical_eddy_diffusivity}
     :param simulate_turbines: simulate turbines, defaults to
         {simulate_turbines}
+    :param turbine_turbulence_model: turbine turbulence model. Set to 
+        ``'delft'`` for the default model or ``'canopy'`` to use the 
+        Réthoré (2009) canopy model. Defaults to {turbine_turbulence_model}
+    :param beta_p: turbine turbulence canopy model "production" coefficient, 
+        :math:`\\beta_p`. Applies to the ``'canopy'`` turbine turbulence model 
+        only. Defaults to {beta_p}.
+    :param beta_d: turbine turbulence canopy model "dissipation" coefficient, 
+        :math:`\\beta_d`. Applies to the ``'canopy'`` turbine turbulence model 
+        only. Defaults to {beta_d}.
+    :param c_epp: turbine turbulence canopy model "production" closure 
+        coefficient, :math:`C_{{\\varepsilon p}}`. Applies to the 
+        ``'canopy'`` turbine turbulence model only. Defaults to {c_epp}.
+    :param c_epd: turbine turbulence canopy model "dissipation" closure 
+        coefficient, :math:`C_{{\\varepsilon d}}`. Applies to the 
+        ``'canopy'`` turbine turbulence model only. Defaults to {c_epd}.
     :param horizontal_momentum_filter: use high-order horizontal momentum 
         filter. Applies to the ``'fm'`` model only. Defaults to 
         {horizontal_momentum_filter}
@@ -102,7 +119,7 @@ class CaseStudy:
     dt_max: OneOrMany[Num] = 1
 
     #: initial time step, in seconds. For the ``'structured'`` model, this 
-    # is the fixed time step
+    #: is the fixed time step
     dt_init: OneOrMany[Num] = 1
 
     turb_pos_x: OneOrMany[Num] = 6 #: turbine x-position, in meters
@@ -112,6 +129,10 @@ class CaseStudy:
     #: inlet boundary discharge, in cubic meters per second
     discharge: OneOrMany[Num] = 6.0574
 
+    #: uniform bed roughness coefficient, as a Manning number, in seconds
+    #: over metres cube rooted
+    bed_roughness: OneOrMany[Num] = 0.023
+
     #: uniform horizontal eddy viscosity, in metres squared per second
     horizontal_eddy_viscosity: OneOrMany[Num] = 1e-06
 
@@ -127,12 +148,36 @@ class CaseStudy:
     #: simulate turbines
     simulate_turbines: OneOrMany[bool] = True
 
+    #: turbine turbulence model. Set to ``'delft'`` for the default model or
+    #: ``'canopy'`` to use the Réthoré (2009) canopy model.
+    turbine_turbulence_model: OneOrMany[str] = 'delft'
+
+    #: turbine turbulence canopy model "production" coefficient, 
+    #: :math:`\beta_p`. Applies to the ``'canopy'`` turbine turbulence model 
+    #: only.
+    beta_p: OneOrMany[Num] = 1.
+
+    #: turbine turbulence canopy model "dissipation" coefficient, 
+    #: :math:`\beta_d`. Applies to the ``'canopy'`` turbine turbulence model 
+    #: only.
+    beta_d: OneOrMany[Num] = 1.84
+
+    #: turbine turbulence canopy model "production" closure coefficient, 
+    #: :math:`C_{\varepsilon p}`. Applies to the ``'canopy'`` turbine 
+    #: turbulence model only.
+    c_epp: OneOrMany[Num] =  0.77
+
+    #: turbine turbulence canopy model "dissipation" closure coefficient, 
+    #: :math:`C_{\varepsilon d}`. Applies to the ``'canopy'`` turbine 
+    #: turbulence model only.
+    c_epd: OneOrMany[Num] = 0.13
+
     #: use high-order horizontal momentum filter. Applies to the ``'fm'`` 
-    # model only
+    #: model only
     horizontal_momentum_filter: OneOrMany[bool] = True
 
     #: interval for simulation progress output, in seconds. Applies to the 
-    # ``'fm'`` model only
+    #: ``'fm'`` model only
     stats_interval: OneOrManyOptional[Num] = None
 
     #:interval for restart file output, in seconds
@@ -141,7 +186,7 @@ class CaseStudy:
     def __post_init__(self):
 
         mutli_values = {n: v for n, v in zip(self.fields, self.values)
-                                                if isinstance(v, Sequence)}
+                                                        if is_sequence(v)}
 
         # Unpack single length sequences
         for name, value in mutli_values.copy().items():
@@ -189,7 +234,7 @@ def get_case(self, index: int = 0) -> CaseStudy:
             return CaseStudy(*self.values)
 
         self._multi_index_check(index)
-        case_values = [value[index] if isinstance(value, Sequence) else value
+        case_values = [value[index] if is_sequence(value) else value
                                                    for value in self.values]
 
         return CaseStudy(*case_values)
@@ -248,7 +293,7 @@ def is_equal(self, other: object,
 
             other_v = other_dict[f]
 
-            if isinstance(v, Sequence):
+            if is_sequence(v):
                 if tuple(v) != tuple(other_v): return False
             else:
                 if v != other_v: return False
@@ -263,6 +308,8 @@ def _multi_index_check(self, index: int):
         if not (-1 * length <= index <= length - 1):
             raise IndexError("index out of range")
 
+
+
     def __eq__(self, other: object) -> bool:
         return self.is_equal(other)
 
@@ -271,12 +318,16 @@ def __getitem__(self, item: int) -> CaseStudy:
 
     def __len__(self) -> int:
 
-        mutli_values = [v for v in self.values if isinstance(v, Sequence)]
+        mutli_values = [v for v in self.values if is_sequence(v)]
 
         if not mutli_values: return 1
         return len(mutli_values[0])
 
 
+def is_sequence(x):
+    return isinstance(x, Sequence) and not isinstance(x, str)
+
+
 @docstringtemplate
 @dataclass(eq=False, frozen=True)
 class MycekStudy(CaseStudy):

diff --git a/src/snl_d3d_cec_verify/template.py b/src/snl_d3d_cec_verify/template.py
@@ -88,9 +88,10 @@ class Template:
     #: variables to ignore in the given :class:`.CaseStudy` objects when
     #: filling templates
     no_template: List[str] = field(
-                                default_factory=lambda: ["dx",
-                                                         "dy",
-                                                         "simulate_turbines"])
+                        default_factory=lambda: ["dx",
+                                                 "dy",
+                                                 "simulate_turbines",
+                                                 "turbine_turbulence_model"])
 
     _template_tmp: tempfile.TemporaryDirectory = field(init=False, repr=False)
     _extras: _BaseTemplateExtras = field(init=False, repr=False)
@@ -218,6 +219,11 @@ def data_hook(self, case: CaseStudy,
             simulate_turbines = ""
 
         data["simulate_turbines"] = simulate_turbines
+
+        assert isinstance(case.turbine_turbulence_model, str)
+        turbine_turbulence_code = _get_turbine_turbulence_code(
+                                            case.turbine_turbulence_model)
+        data["turbine_turbulence_code"] = turbine_turbulence_code
 
     def write_grid(self, project_path: StrOrPath,
                          dx: Num,
@@ -239,12 +245,21 @@ def data_hook(self, case: CaseStudy,
         data["version"] = f"{version('SNL-Delft3D-CEC-Verify')}"
         data["os"] = f"{platform.system()}"
 
+        if case.simulate_turbines:
+            simulate_turbines = "Filtrb = #turbines.ini#"
+        else:
+            simulate_turbines = ""
+
+        data["simulate_turbines"] = simulate_turbines
+
+        assert isinstance(case.turbine_turbulence_model, str)
+        turbine_turbulence_code = _get_turbine_turbulence_code(
+                                            case.turbine_turbulence_model)
+        data["turbine_turbulence_code"] = turbine_turbulence_code
+
         if not case.simulate_turbines:
-            data["simulate_turbines"] = ""
             return
 
-        data["simulate_turbines"] = "Filtrb = #turbines.ini#"
-
         # Inform the type checker that we have Num for single value cases
         dx = cast(Num, case.dx)
         dy = cast(Num, case.dy)
@@ -285,3 +300,15 @@ def write_grid(self, project_path: StrOrPath,
 def _package_template_path(template_type) -> Path:
     this_dir = os.path.dirname(os.path.realpath(__file__))
     return Path(this_dir) / "templates" / template_type
+
+
+def _get_turbine_turbulence_code(turbine_turbulence_model: str) -> int:
+
+    if turbine_turbulence_model == "delft":
+        return 0
+    elif turbine_turbulence_model == "canopy":
+        return 1
+
+    err_msg = ("Unrecognised turbine turbulence model "
+               f"'{turbine_turbulence_model}'")
+    raise ValueError(err_msg)
diff --git a/src/snl_d3d_cec_verify/templates/fm/input/FlowFM.mdu b/src/snl_d3d_cec_verify/templates/fm/input/FlowFM.mdu
@@ -78,9 +78,9 @@ Epshu                             = 0.0001          # Threshold water depth for
 HorizontalMomentumFilter          = {{ '{:<15}'.format(horizontal_momentum_filter) }} # Use high-order horizontal momentum filter  (1: yes, 0: no)
 
 [physics]
-UnifFrictCoef                     = 0.023           # Uniform friction coefficient (0: no friction)
+UnifFrictCoef                     = {{ '{:<15.9G}'.format(bed_roughness) }} # Uniform friction coefficient (0: no friction)
 UnifFrictType                     = 1               # Uniform friction type (0: Chezy, 1: Manning, 2: White-Colebrook)
-UnifFrictCoef1D                   = 0.023           # Uniform friction coefficient in 1D links (0: no friction)
+UnifFrictCoef1D                   = {{ '{:<15.9G}'.format(bed_roughness) }} # Uniform friction coefficient in 1D links (0: no friction)
 UnifFrictCoefLin                  = 0               # Uniform linear friction coefficient for ocean models (0: no friction)
 Umodlin                           = 0               # Linear friction umod, for ifrctyp=4,5,6
 Vicouv                            = {{ '{:<15.9G}'.format(horizontal_eddy_viscosity) }} # Uniform horizontal eddy viscosity

diff --git a/src/snl_d3d_cec_verify/templates/fm/input/turbines.ini b/src/snl_d3d_cec_verify/templates/fm/input/turbines.ini
@@ -12,10 +12,10 @@
     ThrustCurve  = #Turbine Type 1#    # Reference to C_T in curves.trb file
     PowerCurve   = #Turbine Type 1#    # Reference to C_P in curves.trb file
     NDiaDist4Vel = 1.000               # If TurbineModel = 0, reference velocity for momentum sink (turbine) 
-    Beta_p       = 0.0                 # Rethore model tke production
-    Beta_d       = 0.0                 # Rethore model tke dissipation
-    Cep4         = 0.0                 # Rethore model epsilon production
-    Cep5         = 0.0                 # Rethore model epsilon dissipation
+    Beta_p       = {{ '{:<19}'.format(beta_p) }} # Rethore model tke production
+    Beta_d       = {{ '{:<19}'.format(beta_d) }} # Rethore model tke dissipation
+    Cep4         = {{ '{:<19}'.format(c_epp) }} # Rethore model epsilon production
+    Cep5         = {{ '{:<19}'.format(c_epd) }} # Rethore model epsilon dissipation
     TurbineModel    = 1                # Set 0 = Delft-TurbineModel Set 1 = SNL-TurbineModel
-    Turbulencemodel = 0                # Set 0 = Default Delft turbulence model, Set 1 = Rethore turbulence model
+    Turbulencemodel = {{ '{:<15}'.format(turbine_turbulence_code) }}  # Set 0 = Default Delft turbulence model, Set 1 = Rethore turbulence model
     CurvesFil    = #curves.trb#        # location of thrust and power coefficients