D3D Grid Convergence

examples/ADCP_Delft3D_TRTS_2021_2023_all_transects.ipynb

@@ -0,0 +1,16 @@
+{
+ "cells": [],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python_310",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.10.15"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

mhkit/river/io/d3d.py

@@ -35,8 +35,8 @@ def get_all_time(data: netCDF4.Dataset) -> NDArray:
         simulation conditions at that time.
     """
 
-    if not isinstance(data, netCDF4.Dataset):
-        raise TypeError("data must be a NetCDF4 object")
+    if not isinstance(data, (netCDF4.Dataset, xr.Dataset)):
+        raise TypeError("data must be a NetCDF4 object or xarray Dataset")
 
     seconds_run = np.ma.getdata(data.variables["time"][:], False)
 
@@ -244,6 +244,10 @@ def get_layer_data(
                 "name": "mesh2d_nLayers",
                 "coords": data.variables["mesh2d_layer_sigma"][:],
             },
+            "mesh2d_face_x mesh2d_face_y mesh2d_layer_sigma": {
+                "name": "mesh2d_nLayers",
+                "coords": data.variables["mesh2d_layer_sigma"][:],
+            },
             "mesh2d_edge_x mesh2d_edge_y": {
                 "name": "mesh2d_nInterfaces",
                 "coords": data.variables["mesh2d_interface_sigma"][:],
@@ -253,7 +257,7 @@ def get_layer_data(
             data.variables["mesh2d_waterdepth"][time_index, :], False
         )
         waterlevel = np.ma.getdata(data.variables["mesh2d_s1"][time_index, :], False)
-        coords = str(data.variables["waterdepth"].coordinates).split()
+        coords = str(data.variables["mesh2d_waterdepth"].coordinates).split()
 
     elif str(data.variables[variable].coordinates) == "FlowElem_xcc FlowElem_ycc":
         cords_to_layers = {
@@ -637,6 +641,10 @@ def get_all_data_points(
                 "name": "mesh2d_nLayers",
                 "coords": data.variables["mesh2d_layer_sigma"][:],
             },
+            "mesh2d_face_x mesh2d_face_y mesh2d_layer_sigma": {
+                "name": "mesh2d_nLayers",
+                "coords": data.variables["mesh2d_layer_sigma"][:],
+            },
             "mesh2d_edge_x mesh2d_edge_y": {
                 "name": "mesh2d_nInterfaces",
                 "coords": data.variables["mesh2d_interface_sigma"][:],
@@ -886,3 +894,33 @@ def list_variables(data: Union[netCDF4.Dataset, xr.Dataset, xr.DataArray]) -> Li
         "data must be a NetCDF4 Dataset, xarray Dataset, or "
         f"xarray DataArray. Got: {type(data)}"
     )
+
+
+def calculate_grid_convergence_index(
+    fine_grid, coarse_grid, refinement_ratio, factor_of_safety=1.25, order=2
+):
+    """
+    Calculate the Grid Convergence Index (GCI) between two grid sizes. https://www.grc.nasa.gov/WWW/wind/valid/tutorial/spatconv.html
+
+    Parameters
+    ----------
+    fine_grid: numpy.ndarray
+        Results from the finer grid.
+    coarse_grid: numpy.ndarray
+        Results from the coarser grid.
+    refinement_ratio: float
+        Refinement ratio between the grids.
+    order: int
+        Order of accuracy (default is 2).
+
+    Returns
+    -------
+    gci: float
+        Grid Convergence Index (GCI).
+    """
+    # Calculate the approximate relative error
+    error = np.abs((fine_grid - coarse_grid) / fine_grid)
+
+    # Calculate the GCI
+    gci = (factor_of_safety * error) / (refinement_ratio**order - 1)
+    return gci