From 48a29fe60b69966b0101c9bfd303ea24620458c7 Mon Sep 17 00:00:00 2001 From: Mathew Topper Date: Wed, 8 Jun 2022 18:30:48 +0100 Subject: [PATCH] Add Turbulence to Basic and Validation examples (#42) This PR adds turbulent quantities to the `basic.py` and `validation.py` examples. For `basic.py`, plots of the turbulent kinetic energy are added. For `validation.py` comparison is made with newly extracted turbulence intensity values from the Mycek experiment. Other changes are as follows: * The interpolation of the turbulent kinetic energy from the edges to the face centres was improved. Note, however, that the results differ between Linux and Windows platforms. * The Transect class now requires an integer id as its first argument. For the Validate class, all transects stored within must have a unique ID. This change allows consistent "naming" of the stored transects when new ones are added - before they may have been renumbered based on the alphabetical order. * The workflows using setup-miniconda were modified to use native shells where possible as cygwin is not working on Windows, currently. * Bump to version 0.8.0 --- .conda/recipe/meta.yaml | 2 +- .github/workflows/docs.yml | 8 +- .github/workflows/static_type_checks.yml | 8 +- .github/workflows/unit_tests.yml | 8 +- README.md | 4 +- docs/_assets/gh-pages-redirect.html | 4 +- docs/conf.py | 4 +- examples/basic.py | 32 ++++- examples/grid_convergence.py | 13 +- examples/validation.py | 87 +++++++++++-- setup.cfg | 2 +- src/snl_d3d_cec_verify/result/__init__.py | 64 +++++++--- src/snl_d3d_cec_verify/result/faces.py | 58 +++++++-- .../result/mycek2014/11a.yaml | 1 + .../result/mycek2014/11b.yaml | 1 + .../result/mycek2014/11c.yaml | 16 +++ .../result/mycek2014/11d.yaml | 16 +++ .../result/mycek2014/A12a_xstar5.yaml | 1 + .../result/mycek2014/A12b_xstar5.yaml | 1 + test_data/transects/one.yaml | 1 + test_data/transects/two.yaml | 1 + test_data/transects_bad/one.yaml | 8 ++ test_data/transects_bad/two.yaml | 8 ++ tests/test_result.py | 114 ++++++++++++------ tests/test_result_faces.py | 5 +- 25 files changed, 370 insertions(+), 97 deletions(-) create mode 100644 src/snl_d3d_cec_verify/result/mycek2014/11c.yaml create mode 100644 src/snl_d3d_cec_verify/result/mycek2014/11d.yaml create mode 100644 test_data/transects_bad/one.yaml create mode 100644 test_data/transects_bad/two.yaml diff --git a/.conda/recipe/meta.yaml b/.conda/recipe/meta.yaml index da63b9c..5a36d69 100644 --- a/.conda/recipe/meta.yaml +++ b/.conda/recipe/meta.yaml @@ -1,4 +1,4 @@ -{% set version = "0.7.0" %} +{% set version = "0.8.0" %} package: name: snl-delft3d-cec-verify diff --git a/.github/workflows/docs.yml b/.github/workflows/docs.yml index 805736e..d44f40a 100644 --- a/.github/workflows/docs.yml +++ b/.github/workflows/docs.yml @@ -4,12 +4,18 @@ jobs: doctest: runs-on: ${{ matrix.os }} strategy: + fail-fast: false matrix: os: [windows-latest, ubuntu-latest] python-version: ['3.9', '3.10'] + include: + - os: windows-latest + shell: pwsh + - os: ubuntu-latest + shell: bash -l {0} defaults: run: - shell: bash -l {0} + shell: ${{ matrix.shell }} steps: - uses: actions/checkout@v2 with: diff --git a/.github/workflows/static_type_checks.yml b/.github/workflows/static_type_checks.yml index 846b291..f408c44 100644 --- a/.github/workflows/static_type_checks.yml +++ b/.github/workflows/static_type_checks.yml @@ -4,12 +4,18 @@ jobs: mypy: runs-on: ${{ matrix.os }} strategy: + fail-fast: false matrix: os: [windows-latest, ubuntu-latest] python-version: ['3.9', '3.10'] + include: + - os: windows-latest + shell: pwsh + - os: ubuntu-latest + shell: bash -l {0} defaults: run: - shell: bash -l {0} + shell: ${{ matrix.shell }} steps: - uses: actions/checkout@v2 with: diff --git a/.github/workflows/unit_tests.yml b/.github/workflows/unit_tests.yml index d956729..a1d6d13 100644 --- a/.github/workflows/unit_tests.yml +++ b/.github/workflows/unit_tests.yml @@ -4,12 +4,18 @@ jobs: pytest: runs-on: ${{ matrix.os }} strategy: + fail-fast: false matrix: os: [windows-latest, ubuntu-latest] python-version: ['3.9', '3.10'] + include: + - os: windows-latest + shell: pwsh + - os: ubuntu-latest + shell: bash -l {0} defaults: run: - shell: bash -l {0} + shell: ${{ matrix.shell }} steps: - uses: actions/checkout@v2 with: diff --git a/README.md b/README.md index 36ba58b..3412f67 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.7.0 +(base) > conda create -y -n snld3d --override-channels -c conda-forge -c dataonlygreater snl-delft3d-cec-verify=0.8.0 (base) > conda activate snld3d (snld3d) > conda config --env --add channels conda-forge --add channels dataonlygreater (snld3d) > conda config --env --set channel_priority strict @@ -119,7 +119,7 @@ can be optionally converted to Word format if the `pypandoc` package is installed. To install it, type: ``` -(snld3d) > conda install -y pypandoc pandoc=2.12 +(snld3d) > conda install -y pypandoc pandoc ``` Currently, a compiled copy of SNL-Delft3D-CEC or SNL-Delft3D-FM-CEC must be diff --git a/docs/_assets/gh-pages-redirect.html b/docs/_assets/gh-pages-redirect.html index 8a24cf1..162aed0 100644 --- a/docs/_assets/gh-pages-redirect.html +++ b/docs/_assets/gh-pages-redirect.html @@ -3,7 +3,7 @@ Redirecting to latest version - - + + diff --git a/docs/conf.py b/docs/conf.py index ae0ff30..71ddcaf 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.7.0' +release = '0.8.0' # -- 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.7.0' +smv_latest_version = 'v0.8.0' # 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 a445d50..1ae369c 100644 --- a/examples/basic.py +++ b/examples/basic.py @@ -50,30 +50,49 @@ def main(template_type): result = Result(tmpdirname) turb_ds = result.faces.extract_turbine_centre(-1, case) turb_u = turb_ds["$u$"].values.take(0) + turb_k = turb_ds["$k$"].values.take(0) # Record data for table data["discharge"].append(case.discharge) data["u"].append(turb_u) + data["k"].append(turb_k) # Add report section with plot report.content.add_heading( f"Discharge: {case.discharge} (cubic meters per second)", level=2) - fig, ax = plt.subplots() turbz = result.faces.extract_turbine_z(-1, case) + + fig, ax = plt.subplots(figsize=(6, 3)) turbz["$u$"].plot(ax=ax, x="$x$", y="$y$") - plot_name = f"discharge_case_{i}.png" + plot_name = f"discharge_case_{i}_u.png" + plot_path = report_dir / plot_name + plt.savefig(plot_path, bbox_inches="tight") + + report.content.add_image(plot_name, + "Velocity, $u$ (m/s)", + width="5.4in") + + fig, ax = plt.subplots(figsize=(6, 3)) + turbz["$k$"].plot(ax=ax, x="$x$", y="$y$") + plot_name = f"discharge_case_{i}_k.png" plot_path = report_dir / plot_name - plt.savefig(plot_path) + plt.savefig(plot_path, bbox_inches="tight") - report.content.add_image(plot_name, "u-velocity (m/s)") + report.content.add_image(plot_name, + "Turbulent kinetic energy, $k$ " + "(m^2^/s^2^)", + width="5.7in") df = pd.DataFrame(data) + report.content.add_heading("Results", level=2) + report.content.add_table(df, index=False, - caption="Turbine centre velocity per discharge level") + caption="Turbine centre velocity and TKE per " + "discharge level") os_name = platform.system() report.title = f"Basic Example ({os_name})" @@ -91,7 +110,8 @@ def main(template_type): 'docx', outputfile=f"{report_dir / 'report.docx'}", extra_args=[f'--resource-path={report_dir}', - '--reference-doc=reference.docx']) + '--reference-doc=reference.docx'], + sandbox=False) except ImportError: diff --git a/examples/grid_convergence.py b/examples/grid_convergence.py index d95e1a4..61aff44 100644 --- a/examples/grid_convergence.py +++ b/examples/grid_convergence.py @@ -248,6 +248,8 @@ def main(template_type, max_experiments, omp_num_threads): for i, transect in enumerate(global_validate): + if transect.name not in ["$u$", "$u_0$"]: continue + description = transect.attrs['description'] transect_df = transect_grouped.get_group(description).drop("Transect", axis=1) @@ -362,6 +364,8 @@ def main(template_type, max_experiments, omp_num_threads): for i, transect in enumerate(global_validate): + if transect.name not in ["$u$", "$u_0$"]: continue + description = transect.attrs['description'] report.content.add_heading(description, level=3) @@ -446,7 +450,8 @@ def main(template_type, max_experiments, omp_num_threads): extra_args=['-C', f'--resource-path={report_dir}', '--bibliography=examples.bib', - '--reference-doc=reference.docx']) + '--reference-doc=reference.docx'], + sandbox=False) except ImportError: @@ -533,6 +538,8 @@ def plot_transects(case, for i, transect in enumerate(validate): + if transect.name not in ["$u$", "$u_0$"]: continue + transect_true = transect.to_xarray() # Compare transect @@ -569,9 +576,11 @@ def get_rmse(estimated, observed): def get_transect_error(case, validate, result, factor, data): - + for i, transect in enumerate(validate): + if transect.name not in ["$u$", "$u_0$"]: continue + transect_true = transect.to_xarray() # Compare transect diff --git a/examples/validation.py b/examples/validation.py index d374096..6d472ba 100644 --- a/examples/validation.py +++ b/examples/validation.py @@ -40,7 +40,7 @@ def get_d3d_bin_path(): return root.resolve() -def get_u0(da, case, transect): +def get_normalised(da, case, transect): da = get_reset_origin(da, (case.turb_pos_x, case.turb_pos_y, case.turb_pos_z)) @@ -61,6 +61,11 @@ def get_gamma0(da, case, transect): return da +def get_TI(ds): + velmag = np.sqrt(ds["$u$"]**2 + ds["$v$"]**2 + ds["$w$"]**2) + return (100 / velmag) * np.sqrt(2 * ds["$k$"] / 3) + + def get_rmse(estimated, observed): return np.sqrt(((estimated - observed) ** 2).mean()) @@ -75,7 +80,9 @@ def main(template_type): report_dir = Path(template_type) / "validation_report" report_dir.mkdir(exist_ok=True, parents=True) - data = defaultdict(list) + + data_u0 = defaultdict(list) + data_I0 = defaultdict(list) with tempfile.TemporaryDirectory() as tmpdirname: @@ -88,6 +95,8 @@ def main(template_type): for i, transect in enumerate(validate): + if transect.name not in ["$u$", "$u_0$"]: continue + # Compare transect transect_sim = result.faces.extract_z(-1, **transect) transect_true = transect.to_xarray() @@ -100,8 +109,12 @@ def main(template_type): major_axis = f"${transect.attrs['major_axis']}^*$" # Create and save a u0 figure - transect_sim_u0 = get_u0(transect_sim["$u$"], case, transect_true) - transect_true_u0 = get_u0(transect_true, case, transect_true) + transect_sim_u0 = get_normalised(transect_sim["$u$"], + case, + transect_true) + transect_true_u0 = get_normalised(transect_true, + case, + transect_true) fig, ax = plt.subplots(figsize=(4, 2.75), dpi=300) transect_sim_u0.plot(ax=ax, x=major_axis, label='Delft3D') @@ -122,8 +135,8 @@ def main(template_type): # Calculate RMS error and store rmse = get_rmse(transect_sim_u0.values, transect_true_u0.values) - data["Transect"].append(transect.attrs['description']) - data["RMSE"].append(rmse) + data_u0["Transect"].append(transect.attrs['description']) + data_u0["RMSE (m/s)"].append(rmse) # Create and save a gamma0 figure transect_sim_gamma0 = get_gamma0(transect_sim["$u$"], @@ -149,11 +162,64 @@ def main(template_type): "engineered from [@mycek2014, fig. " f"{transect.attrs['figure']}].") report.content.add_image(plot_name, caption, width="4in") + + for i, transect in enumerate(validate): + + if transect.name != "$I_0$": continue + + # Compare transect + transect_sim = result.faces.extract_z(-1, **transect) + transect_sim = transect_sim.assign({"$I$": get_TI}) + transect_true = transect.to_xarray() + + # Add report section with plot + report.content.add_heading(transect.attrs['description'], + level=2) + + # Determine plot x-axis + major_axis = f"${transect.attrs['major_axis']}^*$" + + # Create and save a u0 figure + transect_sim_I0 = get_normalised(transect_sim["$I$"], + case, + transect_true) + transect_true_I0 = get_normalised(transect_true, + case, + transect_true) + + fig, ax = plt.subplots(figsize=(4, 2.75), dpi=300) + transect_sim_I0.plot(ax=ax, x=major_axis, label='Delft3D') + transect_true_I0.plot(ax=ax, x=major_axis, label='Experiment') + ax.legend() + ax.grid() + ax.set_title("") + + plot_name = f"transect_I0_{i}.png" + plot_path = report_dir / plot_name + plt.savefig(plot_path, bbox_inches='tight') + + # Add figure with caption + caption = ("$I_0$ comparison (%). Experimental data " + "reverse engineered from [@mycek2014, fig. " + f"{transect.attrs['figure']}].") + report.content.add_image(plot_name, caption, width="4in") + + # Calculate RMS error and store + rmse = get_rmse(transect_sim_I0.values, transect_true_I0.values) + data_I0["Transect"].append(transect.attrs['description']) + data_I0["RMSE (%)"].append(rmse) - # Add table for errors - df = pd.DataFrame(data) - caption = "Root-mean-square errors in $u_0$." + # Add tables for errors report.content.add_heading("Errors", level=2) + + df = pd.DataFrame(data_u0) + caption = "Root-mean-square errors in $u_0$." + report.content.add_table(df, + index=False, + caption=caption) + + df = pd.DataFrame(data_I0) + caption = "Root-mean-square errors in $I_0$." report.content.add_table(df, index=False, caption=caption) @@ -182,7 +248,8 @@ def main(template_type): extra_args=['-C', f'--resource-path={report_dir}', '--bibliography=examples.bib', - '--reference-doc=reference.docx']) + '--reference-doc=reference.docx'], + sandbox=False) except ImportError: diff --git a/setup.cfg b/setup.cfg index b2f1ddc..3712bc2 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = SNL-Delft3D-CEC-Verify -version = 0.7.0 +version = 0.8.0 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/result/__init__.py b/src/snl_d3d_cec_verify/result/__init__.py index 4ac43f1..c58fe6e 100644 --- a/src/snl_d3d_cec_verify/result/__init__.py +++ b/src/snl_d3d_cec_verify/result/__init__.py @@ -7,6 +7,7 @@ import itertools from abc import abstractmethod from typing import (Any, + Dict, Hashable, List, Mapping, @@ -289,7 +290,9 @@ class Validate(): Validate(0: Centreline velocity (3\\% TI) 1: Centreline velocity (15\\% TI) 2: Axial velocity at $x^*=5$ (3\\% TI) - 3: Axial velocity at $x^*=5$ (15\\% TI)) + 3: Axial velocity at $x^*=5$ (15\\% TI) + 4: Centreline turbulence intensity (3\\% TI) + 5: Centreline turbulence intensity (15\\% TI)) >>> validate[0].to_xarray() #doctest: +ELLIPSIS @@ -326,7 +329,8 @@ class Validate(): """ - _transects: List[Transect] = field(default_factory=list, init=False) + _transects: Dict[Any, Transect] = field(default_factory=dict, init=False) + _last_key_idx: int = field(default=-1, init=False) case: InitVar[Optional[CaseStudy]] = None #: :meta private: data_dir: InitVar[Optional[StrOrPath]] = None #: :meta private: @@ -354,7 +358,7 @@ def __post_init__(self, case: Optional[CaseStudy], turb_pos_z = case.turb_pos_z translation = (turb_pos_x, turb_pos_y, turb_pos_z) - transects = [] + transects = {} for item in sorted(data_dir.iterdir()): @@ -362,13 +366,33 @@ def __post_init__(self, case: Optional[CaseStudy], if not item.suffix == '.yaml': continue transect = Transect.from_yaml(item, translation) - transects.append(transect) + + if transect.id in transects: + err_msg = (f"Transect ID '{transect.id}' given in {str(item)} " + "is already used") + raise RuntimeError(err_msg) + + transects[transect.id] = transect object.__setattr__(self, '_transects', transects) def __getitem__(self, item: int) -> Transect: return self._transects[item] + def __iter__(self): + return self + + def __next__(self): + + keys = sorted(self._transects.keys()) + object.__setattr__(self, '_last_key_idx', self._last_key_idx + 1) + + if self._last_key_idx == len(keys): + object.__setattr__(self, '_last_key_idx', -1) + raise StopIteration + else: + return self._transects[keys[self._last_key_idx]] + def __len__(self) -> int: return len(self._transects) @@ -379,12 +403,15 @@ def __repr__(self): if self._transects: - transect = self._transects[0] - msg += f"0: {transect.attrs['description']}" + keys = sorted(self._transects.keys()) - for i, transect in enumerate(self._transects[1:]): + transect = self._transects[keys[0]] + msg += f"{transect.id}: {transect.attrs['description']}" + + for key in keys[1:]: + transect = self._transects[key] msg += "\n" + " " * indent - msg += f"{i + 1}: {transect.attrs['description']}" + msg += f"{transect.id}: {transect.attrs['description']}" msg += ")" @@ -409,7 +436,7 @@ class Transect(): Data is stored for each x and y pair, in order. For example: - >>> x = Transect(-1, [1, 2, 3, 4], [2, 2, 2, 2], [5, 4, 3, 2]) + >>> x = Transect(0, -1, [1, 2, 3, 4], [2, 2, 2, 2], [5, 4, 3, 2]) >>> x.to_xarray() #doctest: +ELLIPSIS array([5, 4, 3, 2]) @@ -439,8 +466,9 @@ class Transect(): $u$ (dim_0) float64 0.7793 0.7776 0.7766 0.7757 $v$ (dim_0) float64 1.193e-17 4.679e-17 2.729e-17 -2.519e-17 $w$ (dim_0) float64 -0.001658 0.0001347 -0.00114 0.0002256 - $k$ (dim_0) float64 0.004756 0.004067 0.003406 0.003925 + $k$ (dim_0) float64 0.0047... 0.0046... 0.004... 0.0045... + :param id: integer identifier for the transect :param z: z-level of the transect, in meters :param x: x-coordinates of the transect, in meters :param y: y-coordinates of the transect, in meters @@ -455,6 +483,9 @@ class Transect(): """ + #: id for the transect + id: int + #: z-level of the transect, in meters z: Num = field(repr=False) @@ -500,6 +531,7 @@ def __post_init__(self, translation: Vector): @docstringtemplate @classmethod def from_csv(cls: Type[T], path: StrOrPath, + id: int, name: Optional[str] = None, attrs: Optional[dict[str, str]] = None, translation: Vector = (0, 0, 0)) -> T: @@ -514,6 +546,7 @@ def from_csv(cls: Type[T], path: StrOrPath, 9, 3, 0, 3 :param path: path to the CSV file to load + :param id: integer identifier :param name: name of data in the resulting :class:`.Transect` object :param attrs: attributes for the resulting :class:`.Transect` object :param translation: translation of the transect origin, defaults to @@ -546,7 +579,8 @@ def from_csv(cls: Type[T], path: StrOrPath, if attrs is None: attrs = {} attrs["path"] = str(path) - return cls(z[0], + return cls(id, + z[0], np.array(cols["x"]), np.array(cols["y"]), data=data, @@ -560,8 +594,9 @@ def from_yaml(cls: Type[T], path: StrOrPath, translation: Vector = (0, 0, 0)) -> T: """Create a new :class:`.Transect` object from a YAML file. - The YAML file must have ``z`` , ``x`` and ``y``, and keys where ``z`` - is a single value and ``x`` and ``y`` are arrays. + The YAML file must have ``id``, ``z`` , ``x`` and ``y``, and keys + where ``id`` is an integer identifier, ``z`` is a single value and + ``x`` and ``y`` are arrays. Optionally, a ``data`` key can be given as an array, ``name`` key as a single value (treated as a string) and an ``attrs`` key with a nested @@ -598,7 +633,8 @@ def from_yaml(cls: Type[T], path: StrOrPath, if "attrs" in raw: attrs = raw["attrs"] attrs["path"] = str(path) - return cls(raw["z"], + return cls(raw["id"], + raw["z"], x=np.array(raw["x"]), y=np.array(raw["y"]), data=data, diff --git a/src/snl_d3d_cec_verify/result/faces.py b/src/snl_d3d_cec_verify/result/faces.py index c21414a..dbfe6e9 100644 --- a/src/snl_d3d_cec_verify/result/faces.py +++ b/src/snl_d3d_cec_verify/result/faces.py @@ -2,6 +2,8 @@ from __future__ import annotations +import warnings +import itertools import collections from abc import ABC, abstractmethod from typing import (cast, @@ -91,7 +93,7 @@ class Faces(ABC, _FacesDataClassMixin): $u$ ($x$, $y$) float64 0.781 0.781 0.781 ... 0.7763 0.7763 0.7763 $v$ ($x$, $y$) float64 -3.237e-18 1.423e-17 ... -8.598e-17 -4.824e-17 $w$ ($x$, $y$) float64 -0.01472 -0.01472 ... 0.001343 0.001343 - $k$ ($x$, $y$) float64 0.004796 0.004796 ... 0.003683 0.003683 + $k$ ($x$, $y$) float64 0.004796 0.004796 ... 0.00368... 0.00368... :param nc_path: path to the ``.nc`` file containing results :param n_steps: number of time steps in the simulation @@ -135,7 +137,7 @@ def extract_turbine_centre(self, t_step: int, $u$ (dim_0) float64 0.7748 $v$ (dim_0) float64 -2.942e-17 $w$ (dim_0) float64 0.0002786 - $k$ (dim_0) float64 0.004307 + $k$ (dim_0) float64 0.004... The position extracted can also be shifted using the ``offset_x``, ``offset_y`` and ``offset_z`` parameters. @@ -209,7 +211,7 @@ def extract_turbine_centreline(self, t_step: int, $u$ (dim_0) float64 0.7748 0.7747 0.7745 0.7745 ... 0.7759 0.7762 nan $v$ (dim_0) float64 -2.942e-17 4.192e-17 9.126e-17 ... -8.523e-17 nan $w$ (dim_0) float64 0.0002786 -0.0004764 0.0003097 ... -7.294e-05 nan - $k$ (dim_0) float64 0.004307 0.004222 0.004162 ... 0.003691 nan + $k$ (dim_0) float64 0.004... 0.0042... 0.00417... 0.00370... nan The position extracted can also be shifted using the ``offset_x``, ``offset_y`` and ``offset_z`` parameters. @@ -279,7 +281,7 @@ def extract_turbine_z(self, t_step: int, $u$ ($x$, $y$) float64 0.781 0.781 0.781 ... 0.7763 0.7763 0.7763 $v$ ($x$, $y$) float64 -3.237e-18 1.423e-17 ... -8.598e-17 -4.824e-17 $w$ ($x$, $y$) float64 -0.01472 -0.01472 ... 0.001343 0.001343 - $k$ ($x$, $y$) float64 0.004796 0.004796 ... 0.003683 0.003683 + $k$ ($x$, $y$) float64 0.004796 0.004796 ... 0.00368... 0.00368... The z-plane can be shifted using the ``offset_z`` parameter. @@ -339,7 +341,7 @@ def extract_z(self, t_step: int, $u$ (dim_0) float64 0.7748 0.7747 0.7745 0.7745 0.7746 $v$ (dim_0) float64 -3.877e-18 4.267e-17 5.452e-17 5.001e-17 8.011e-17 $w$ (dim_0) float64 0.0002786 -0.0004764 ... -0.0002754 0.0003252 - $k$ (dim_0) float64 0.003764 0.003686 0.004169 0.004107 0.003526 + $k$ (dim_0) float64 0.004... 0.0042... 0.00417... 0.004... 0.00402... If ``x`` and ``y`` are not given, then the results are returned at the face centres. @@ -357,7 +359,7 @@ def extract_z(self, t_step: int, $u$ ($x$, $y$) float64 0.781 0.781 0.781 ... 0.7763 0.7763 0.7763 $v$ ($x$, $y$) float64 -3.237e-18 1.423e-17 ... -8.598e-17 -4.824e-17 $w$ ($x$, $y$) float64 -0.01472 -0.01472 ... 0.001343 0.001343 - $k$ ($x$, $y$) float64 0.004796 0.004796 ... 0.003683 0.003683 + $k$ ($x$, $y$) float64 0.004796 0.004796 ... 0.00368... 0.00368... :param t_step: Time step index :param z: z-level at which to extract data @@ -415,7 +417,7 @@ def extract_sigma(self, t_step: int, $u$ (dim_0) float64 0.7747 0.7746 0.7744 0.7745 0.7745 $v$ (dim_0) float64 -3.88e-18 4.267e-17 5.452e-17 5.002e-17 8.013e-17 $w$ (dim_0) float64 0.0002791 -0.0004769 ... -0.0002756 0.0003256 - $k$ (dim_0) float64 0.003767 0.003689 0.004172 0.004109 0.003528 + $k$ (dim_0) float64 0.004... 0.0042... 0.0041... 0.004... 0.0040... If ``x`` and ``y`` are not given, then the results are returned at the face centres. @@ -433,7 +435,7 @@ def extract_sigma(self, t_step: int, $u$ ($x$, $y$) float64 0.7809 0.7809 0.7809 ... 0.7763 0.7763 0.7763 $v$ ($x$, $y$) float64 -3.29e-18 1.419e-17 ... -8.598e-17 -4.824e-17 $w$ ($x$, $y$) float64 -0.01473 -0.01473 ... 0.001343 0.001343 - $k$ ($x$, $y$) float64 0.004803 0.004803 ... 0.003683 0.003683 + $k$ ($x$, $y$) float64 0.004803 0.004803 ... 0.00368... 0.00368... :param t_step: Time step index :param sigma: sigma-level at which to extract data @@ -635,7 +637,6 @@ def _map_to_faces_frame_with_tke(map_path: StrOrPath, edgest["y"] = edgest['geometry'].apply( lambda line: np.array(line.centroid.coords[0])[1]) edgesdf = pd.DataFrame(edgest[["x", "y", "sigma", "u1", "turkin1"]]) - edgesdf['turkin1'] = edgesdf['turkin1'].fillna(0) facest = facest.set_index(["x", "y", "sigma"]) facest = facest.sort_index() @@ -648,13 +649,44 @@ def _map_to_faces_frame_with_tke(map_path: StrOrPath, for sigma, group in edgesdf.groupby(by="sigma"): + # Fill missing values + groupna = group[pd.isna(group["turkin1"])] + group = group[~pd.isna(group["turkin1"])] + + points = np.array(list(zip(group.x, group.y))) + values = group.turkin1.values + + group_x = sorted(groupna.x.unique()) + group_y = sorted(groupna.y.unique()) + group_grid_x, group_grid_y = np.meshgrid(group_x, group_y) + + group_grid_z = interpolate.griddata(points, + values, + (group_grid_x, group_grid_y), + method='nearest') + + turkin1 = [] + + for i, j in itertools.product(range(len(group_x)), + range(len(group_y))): + turkin1.append(group_grid_z[j, i]) + + groupna = groupna.set_index(["x", "y"]) + groupna = groupna.sort_index() + groupna["turkin1"] = turkin1 + group = pd.concat([group, groupna.reset_index()]) + + # Interpolate onto faces grid points = np.array(list(zip(group.x, group.y))) values = group.turkin1.values - grid_z = interpolate.griddata(points, - values, - (grid_x, grid_y), - method='linear') + # Can warn about s being too low + with warnings.catch_warnings(): + warnings.simplefilter('ignore', RuntimeWarning) + f = interpolate.interp2d(points[:, 0], + points[:, 1], + values) + grid_z = f(x, y) data = {"x": grid_x.ravel(), "y": grid_y.ravel(), diff --git a/src/snl_d3d_cec_verify/result/mycek2014/11a.yaml b/src/snl_d3d_cec_verify/result/mycek2014/11a.yaml index 9586774..547e2b9 100644 --- a/src/snl_d3d_cec_verify/result/mycek2014/11a.yaml +++ b/src/snl_d3d_cec_verify/result/mycek2014/11a.yaml @@ -1,3 +1,4 @@ +id: 0 z: 0 x: [ 0.84, 1.4 , 2.1 , 2.8 , 3.5 , 4.2 , 4.9 , 5.6 , 6.3 , 7. ] y: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] diff --git a/src/snl_d3d_cec_verify/result/mycek2014/11b.yaml b/src/snl_d3d_cec_verify/result/mycek2014/11b.yaml index 7661ec1..5d6d1f9 100644 --- a/src/snl_d3d_cec_verify/result/mycek2014/11b.yaml +++ b/src/snl_d3d_cec_verify/result/mycek2014/11b.yaml @@ -1,3 +1,4 @@ +id: 1 z: 0 x: [ 0.84, 1.4 , 2.1 , 2.8 , 3.5 , 4.2 , 4.9 , 5.6 , 6.3 , 7. ] y: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] diff --git a/src/snl_d3d_cec_verify/result/mycek2014/11c.yaml b/src/snl_d3d_cec_verify/result/mycek2014/11c.yaml new file mode 100644 index 0000000..2723d2c --- /dev/null +++ b/src/snl_d3d_cec_verify/result/mycek2014/11c.yaml @@ -0,0 +1,16 @@ +id: 4 +z: 0 +x: [ 0.84, 1.4 , 2.1 , 2.8 , 3.5 , 4.2 , 4.9 , 5.6 , 6.3 , 7. ] +y: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] +data: [21.89514033, 16.84147933, 13.96882743, 17.59540208, 19.66209161, + 18.97314855, 15.73663828, 14.16381307, 12.38310869, 11.53818643] +name: $I_0$ +attrs: + description: Centreline turbulence intensity (3\% TI) + $D$: 0.7 + $U_\infty$: 0.8 + $I_\infty$: 3\% + turbine_position: [6, 2, -1] + doi: 10.1016/j.renene.2013.12.036 + figure: 11c + major_axis: x diff --git a/src/snl_d3d_cec_verify/result/mycek2014/11d.yaml b/src/snl_d3d_cec_verify/result/mycek2014/11d.yaml new file mode 100644 index 0000000..3688060 --- /dev/null +++ b/src/snl_d3d_cec_verify/result/mycek2014/11d.yaml @@ -0,0 +1,16 @@ +id: 5 +z: 0 +x: [ 0.84, 1.4 , 2.1 , 2.8 , 3.5 , 4.2 , 4.9 , 5.6 , 6.3 , 7. ] +y: [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.] +data: [32.08633094, 28.41726619, 21.51079137, 18.92086331, 16.65467626, + 17.0323741, 13.57913669, 13.36330935, 13.09352518, 13.30935252] +name: $I_0$ +attrs: + description: Centreline turbulence intensity (15\% TI) + $D$: 0.7 + $U_\infty$: 0.8 + $I_\infty$: 15\% + turbine_position: [6, 2, -1] + doi: 10.1016/j.renene.2013.12.036 + figure: 11d + major_axis: x diff --git a/src/snl_d3d_cec_verify/result/mycek2014/A12a_xstar5.yaml b/src/snl_d3d_cec_verify/result/mycek2014/A12a_xstar5.yaml index 98570b1..515c1f0 100644 --- a/src/snl_d3d_cec_verify/result/mycek2014/A12a_xstar5.yaml +++ b/src/snl_d3d_cec_verify/result/mycek2014/A12a_xstar5.yaml @@ -1,3 +1,4 @@ +id: 2 z: 0 x: [3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5] diff --git a/src/snl_d3d_cec_verify/result/mycek2014/A12b_xstar5.yaml b/src/snl_d3d_cec_verify/result/mycek2014/A12b_xstar5.yaml index a7471f7..77854ad 100644 --- a/src/snl_d3d_cec_verify/result/mycek2014/A12b_xstar5.yaml +++ b/src/snl_d3d_cec_verify/result/mycek2014/A12b_xstar5.yaml @@ -1,3 +1,4 @@ +id: 3 z: 0 x: [3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5] diff --git a/test_data/transects/one.yaml b/test_data/transects/one.yaml index a011f18..71823e7 100644 --- a/test_data/transects/one.yaml +++ b/test_data/transects/one.yaml @@ -1,3 +1,4 @@ +id: 0 z: 0 x: [ 1, 2, 3 ] y: [ 1, 1, 1 ] diff --git a/test_data/transects/two.yaml b/test_data/transects/two.yaml index 98579db..33e1f82 100644 --- a/test_data/transects/two.yaml +++ b/test_data/transects/two.yaml @@ -1,3 +1,4 @@ +id: 1 z: 1 x: [ 1, 2, 3 ] y: [ 1, 1, 1 ] diff --git a/test_data/transects_bad/one.yaml b/test_data/transects_bad/one.yaml new file mode 100644 index 0000000..71823e7 --- /dev/null +++ b/test_data/transects_bad/one.yaml @@ -0,0 +1,8 @@ +id: 0 +z: 0 +x: [ 1, 2, 3 ] +y: [ 1, 1, 1 ] +data: [0, 0, 1] +attrs: + description: mock 1 + diff --git a/test_data/transects_bad/two.yaml b/test_data/transects_bad/two.yaml new file mode 100644 index 0000000..72c1edc --- /dev/null +++ b/test_data/transects_bad/two.yaml @@ -0,0 +1,8 @@ +id: 0 +z: 1 +x: [ 1, 2, 3 ] +y: [ 1, 1, 1 ] +data: [0, 0, 2] +attrs: + description: mock 2 + diff --git a/tests/test_result.py b/tests/test_result.py index 97f5cae..e04fd1e 100644 --- a/tests/test_result.py +++ b/tests/test_result.py @@ -200,7 +200,7 @@ def test_result__repr__(result, data_dir): def test_transect_xy_length_mismatch(): with pytest.raises(ValueError) as excinfo: - Transect(z=1, x=[1, 2, 3], y=[1, 1]) + Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1]) assert "Length of x and y must match" in str(excinfo) @@ -208,59 +208,68 @@ def test_transect_xy_length_mismatch(): def test_transect_data_length_mismatch(): with pytest.raises(ValueError) as excinfo: - Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], data=[2]) + Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], data=[2]) assert "Length of data must match x and y" in str(excinfo) @pytest.mark.parametrize("other, expected", [ (1, False), - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1]), True), - (Transect(z=0, x=[1, 2, 3], y=[1, 1, 1]), False), - (Transect(z=1, x=[0, 2, 3], y=[1, 1, 1]), False), - (Transect(z=1, x=[1, 2, 3], y=[0, 1, 1]), False), - (Transect(z=1, + (Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1]), True), + (Transect(id=0, z=0, x=[1, 2, 3], y=[1, 1, 1]), False), + (Transect(id=0, z=1, x=[0, 2, 3], y=[1, 1, 1]), False), + (Transect(id=0, z=1, x=[1, 2, 3], y=[0, 1, 1]), False), + (Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 0]), False)]) def test_transect_eq(other, expected): - test = Transect(z=1, x=[1, 2, 3], y=[1, 1, 1]) + test = Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1]) assert (test == other) is expected @pytest.mark.parametrize("other, expected", [ - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 0]), True), - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], data=[1, 0, 0]), False)]) + (Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 0]), True), + (Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], data=[1, 0, 0]), False)]) def test_transect_eq_data(other, expected): - test = Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 0]) + test = Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 0]) assert (test == other) is expected @pytest.mark.parametrize("other, expected", [ - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], name="mock"), True), - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], name="not mock"), False), - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1]), False)]) + (Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], name="mock"), True), + (Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], name="not mock"), False), + (Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1]), False)]) def test_transect_eq_name(other, expected): - test = Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], name="mock") + test = Transect(id=0, z=1, x=[1, 2, 3], y=[1, 1, 1], name="mock") assert (test == other) is expected @pytest.mark.parametrize("other, expected", [ - (Transect(z=1, + (Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1], attrs={"mock": "mock"}), True), - (Transect(z=1, + (Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1], attrs={"mock": "not mock"}), False), - (Transect(z=1, + (Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1], attrs={"not mock": "mock"}), False), - (Transect(z=1, x=[1, 2, 3], y=[1, 1, 1]), False)]) + (Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1]), False)]) def test_transect_eq_attrs(other, expected): - test = Transect(z=1, x=[1, 2, 3], y=[1, 1, 1], attrs={"mock": "mock"}) + test = Transect(id=0, + z=1, + x=[1, 2, 3], + y=[1, 1, 1], + attrs={"mock": "mock"}) assert (test == other) is expected @@ -277,8 +286,9 @@ def test_transect_from_csv(tmp_path, mocker): mocker.patch('snl_d3d_cec_verify.result.open', mocker.mock_open(read_data=csv)) - test = Transect.from_csv(d, name="mock") - expected = Transect(z=0, + test = Transect.from_csv(d, 0, name="mock") + expected = Transect(id=0, + z=0, x=[7, 8, 9], y=[3, 3, 3], name="mock", @@ -301,10 +311,12 @@ def test_transect_from_csv_attrs(tmp_path, mocker): mocker.mock_open(read_data=csv)) test = Transect.from_csv(d, + 0, name="mock", attrs={"mock": "mock", "path": "mock"}) - expected = Transect(z=0, + expected = Transect(id=0, + z=0, x=[7, 8, 9], y=[3, 3, 3], name="mock", @@ -327,8 +339,9 @@ def test_transect_from_csv_translation(tmp_path, mocker): mocker.patch('snl_d3d_cec_verify.result.open', mocker.mock_open(read_data=csv)) - test = Transect.from_csv(d, translation=(6, 2, -1)) - expected = Transect(z=0, + test = Transect.from_csv(d, 0, translation=(6, 2, -1)) + expected = Transect(id=0, + z=0, x=[7, 8, 9], y=[3, 3, 3], attrs={"path": str(d.resolve())}) @@ -349,8 +362,9 @@ def test_transect_from_csv_with_data(tmp_path, mocker): mocker.patch('snl_d3d_cec_verify.result.open', mocker.mock_open(read_data=csv)) - test = Transect.from_csv(d) - expected = Transect(z=0, + test = Transect.from_csv(d, 0) + expected = Transect(id=0, + z=0, x=[7, 8, 9], y=[3, 3, 3], data=[1, 2, 3], @@ -373,14 +387,15 @@ def test_transect_from_csv_multi_z_error(tmp_path, mocker): mocker.mock_open(read_data=csv)) with pytest.raises(ValueError) as excinfo: - Transect.from_csv(d) + Transect.from_csv(d, 0) assert "only supports fixed z-value" in str(excinfo) def test_transect_from_yaml(tmp_path, mocker): - text = ("z: -1.0\n" + text = ("id: 0\n" + "z: -1.0\n" "x: [7, 8, 9]\n" "y: [3, 3, 3]\n") @@ -391,7 +406,8 @@ def test_transect_from_yaml(tmp_path, mocker): mocker.mock_open(read_data=text)) test = Transect.from_yaml(d) - expected = Transect(z=-1, + expected = Transect(id=0, + z=-1, x=[7, 8, 9], y=[3, 3, 3], attrs={"path": str(d.resolve())}) @@ -401,7 +417,8 @@ def test_transect_from_yaml(tmp_path, mocker): def test_transect_from_yaml_optional(tmp_path, mocker): - text = ("z: -1.0\n" + text = ("id: 0\n" + "z: -1.0\n" "x: [7, 8, 9]\n" "y: [3, 3, 3]\n" "data: [1, 2, 3]\n" @@ -417,7 +434,8 @@ def test_transect_from_yaml_optional(tmp_path, mocker): mocker.mock_open(read_data=text)) test = Transect.from_yaml(d) - expected = Transect(z=-1, + expected = Transect(id=0, + z=-1, x=[7, 8, 9], y=[3, 3, 3], data=[1, 2, 3], @@ -425,14 +443,12 @@ def test_transect_from_yaml_optional(tmp_path, mocker): attrs={"path": str(d.resolve()), "mock": "mock"}) - print(test) - print(expected) - assert test == expected @pytest.fixture def transect(): - return Transect(z=1, + return Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 1], @@ -471,7 +487,8 @@ def test_trasect_to_xarray(dataarray, transect): def test_trasect_to_xarray_no_data(): - transect = Transect(z=1, + transect = Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1]) result = transect.to_xarray() @@ -496,6 +513,17 @@ def test_validate_empty(tmp_path): assert repr(test) == "Validate()" +def test_validate_id_clash(data_dir): + + transects_path = data_dir / "transects_bad" + + with pytest.raises(RuntimeError) as excinfo: + Validate(data_dir=transects_path) + + assert "Transect ID '0'" in str(excinfo) + assert "is already used" in str(excinfo) + + def test_validate_mycek(): test = Validate() assert len(test) >= 0 @@ -516,7 +544,8 @@ def test_validate_repr(validate): def test_validate_get_item(validate): transect = validate[0] - expected = Transect(z=-1, + expected = Transect(id=0, + z=-1, x=[7, 8, 9], y=[4, 4, 4], data=[0, 0, 1], @@ -526,6 +555,12 @@ def test_validate_get_item(validate): assert transect == expected +def test_validate_iterate(validate): + test = [i for i, _ in enumerate(validate)] + expected = [0, 1] + assert test == expected + + def test_get_reset_origin(dataarray): result = get_reset_origin(dataarray, (1, 1, 1)) @@ -558,7 +593,8 @@ def test_get_normalised_data(dataarray): def test_get_normalised_data_latex(): - transect = Transect(z=1, + transect = Transect(id=0, + z=1, x=[1, 2, 3], y=[1, 1, 1], data=[0, 0, 1], diff --git a/tests/test_result_faces.py b/tests/test_result_faces.py index b2bcc16..850421e 100644 --- a/tests/test_result_faces.py +++ b/tests/test_result_faces.py @@ -440,7 +440,7 @@ def test_map_to_faces_frame_with_tke(data_dir): assert faces_frame["w"].min() > -0.02 assert faces_frame["w"].max() < 0.02 assert faces_frame["tke"].min() > 0 - assert faces_frame["tke"].max() < 0.0085 + assert faces_frame["tke"].max() < 0.0089 sigma_slice = _faces_frame_to_slice(faces_frame, pd.Timestamp('2001-01-01 01:00:00'), @@ -448,6 +448,7 @@ def test_map_to_faces_frame_with_tke(data_dir): -0.75) assert np.isclose(sigma_slice["$z$"].values.mean(), -1.5009617997833038) + assert round(sigma_slice["$k$"].values.mean(), 6) == 0.006289 def test_map_to_faces_frame_with_tke_none(data_dir): @@ -496,7 +497,7 @@ def test_map_to_faces_frame_with_tke_none(data_dir): assert faces_frame["w"].min() > -0.02 assert faces_frame["w"].max() < 0.02 assert faces_frame["tke"].min() > 0 - assert faces_frame["tke"].max() < 0.0085 + assert faces_frame["tke"].max() < 0.0089 def test_map_to_faces_frame(data_dir):