diff --git a/.conda/recipe/meta.yaml b/.conda/recipe/meta.yaml index e1f8ac0..ad46a4f 100644 --- 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 index 50f5aa5..565354c 100644 --- 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 index a6ef3e4..f43c4cb 100644 --- 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 index 1ae369c..48756f9 100644 --- a/examples/basic.py +++ b/examples/basic.py @@ -4,6 +4,7 @@ import platform import tempfile from pathlib import Path +from contextlib import contextmanager from collections import defaultdict import pandas as pd @@ -11,23 +12,8 @@ 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 index e014f4e..609f60e 100644 --- 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 = mathew.topper@dataonlygreater.com 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 index 09b4812..382238c 100644 --- 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 index 4b4a13c..4b39f1a 100644 --- 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 index 5822bcd..948cf9f 100644 --- 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 index 5641347..b1f27a8 100644 --- 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 diff --git a/src/snl_d3d_cec_verify/templates/structured/D3D.mdf b/src/snl_d3d_cec_verify/templates/structured/D3D.mdf index 25845a3..c68739e 100644 --- a/src/snl_d3d_cec_verify/templates/structured/D3D.mdf +++ b/src/snl_d3d_cec_verify/templates/structured/D3D.mdf @@ -37,8 +37,8 @@ Fclou = 0.0000000e+000 Sarea = 0.0000000e+000 Temint = #Y# Roumet = #M# -Ccofu = 2.3000000e-002 -Ccofv = 2.3000000e-002 +Ccofu = {{ '{:<13.7E}'.format(bed_roughness) }} +Ccofv = {{ '{:<13.7E}'.format(bed_roughness) }} Xlo = 0.0000000e+000 Htur2d = #N# Vicouv = {{ '{:<13.7E}'.format(horizontal_eddy_viscosity) }} diff --git a/src/snl_d3d_cec_verify/templates/structured/turbines.ini b/src/snl_d3d_cec_verify/templates/structured/turbines.ini index 54329e2..097e02b 100644 --- a/src/snl_d3d_cec_verify/templates/structured/turbines.ini +++ b/src/snl_d3d_cec_verify/templates/structured/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 diff --git a/tests/test_cases.py b/tests/test_cases.py index eaf7146..01df072 100644 --- a/tests/test_cases.py +++ b/tests/test_cases.py @@ -34,11 +34,17 @@ def test_casestudy_fields(): 'turb_pos_y', 'turb_pos_z', 'discharge', + 'bed_roughness', 'horizontal_eddy_viscosity', 'horizontal_eddy_diffusivity', 'vertical_eddy_viscosity', 'vertical_eddy_diffusivity', 'simulate_turbines', + 'turbine_turbulence_model', + 'beta_p', + 'beta_d', + 'c_epp', + 'c_epd', 'horizontal_momentum_filter', 'stats_interval', 'restart_interval'] @@ -71,11 +77,17 @@ def test_casestudy_values(cases): 3, -1, 6.0574, + 0.023, 1e-06, 1e-06, 1e-06, 1e-06, True, + 'delft', + 1., + 1.84, + 0.77, + 0.13, True, None, 0] diff --git a/tests/test_template.py b/tests/test_template.py index 159dba7..0a99d70 100644 --- a/tests/test_template.py +++ b/tests/test_template.py @@ -12,6 +12,7 @@ from snl_d3d_cec_verify.template import (_FMTemplateExtras, _StructuredTemplateExtras, _package_template_path, + _get_turbine_turbulence_code, Template) @@ -131,6 +132,21 @@ def test_package_template_path(): assert str(expected) in str(result) +@pytest.mark.parametrize("turbine_turbulence_model, expected", [ + ("delft", 0), + ("canopy", 1)]) +def test_get_turbine_turbulence_code(turbine_turbulence_model, expected): + assert _get_turbine_turbulence_code(turbine_turbulence_model) == expected + + +def test_get_turbine_turbulence_code_bad_model(): + + with pytest.raises(ValueError) as excinfo: + _get_turbine_turbulence_code("mock") + + assert "Unrecognised turbine turbulence model 'mock'" in str(excinfo) + + def test_template_bad_type(): with pytest.raises(ValueError) as excinfo: @@ -161,7 +177,7 @@ def test_template_fm_call(mocker, tmp_path): case = CaseStudy() project_path = "mock_template" - excluded_fields = ["dx", "dy"] + excluded_fields = ["dx", "dy", "turbine_turbulence_model"] expected_fields = [field for field in case.fields if field not in excluded_fields] @@ -203,7 +219,7 @@ def test_template_structured_call(mocker, tmp_path): case = CaseStudy() project_path = "mock_template" - excluded_fields = ["dx", "dy"] + excluded_fields = ["dx", "dy", "turbine_turbulence_model"] expected_fields = [field for field in case.fields if field not in excluded_fields]