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
Data-Only-Greater · Jun 8, 2022 · 48a29fe · 48a29fe
1 parent 498e813
commit 48a29fe
Show file tree

Hide file tree

Showing 25 changed files with 370 additions and 97 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -3,7 +3,7 @@
   <head>
     <title>Redirecting to latest version</title>
     <meta charset="utf-8">
-    <meta http-equiv="refresh" content="0; url=./v0.7.0/index.html">
-    <link rel="canonical" href="https://data-only-greater.github.io/SNL-Delft3D-CEC-Verify/v0.7.0/index.html">
+    <meta http-equiv="refresh" content="0; url=./v0.8.0/index.html">
+    <link rel="canonical" href="https://data-only-greater.github.io/SNL-Delft3D-CEC-Verify/v0.8.0/index.html">
   </head>
 </html>
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.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
@@ -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
@@ -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
@@ -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
@@ -1,6 +1,6 @@
 [metadata]
 name = SNL-Delft3D-CEC-Verify
-version = 0.7.0
+version = 0.8.0
 author = Mathew Topper
 author_email = [email protected]
 description = Automated verification of SNL-Delft3D-CEC based on the 2014 Mycek experiment