D3D updates

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)
 
@@ -118,8 +118,8 @@ def _convert_time(
         provided, returns the closest matching time_index.
     """
 
-    if not isinstance(data, netCDF4.Dataset):
-        raise TypeError("data must be NetCDF4 object")
+    if not isinstance(data, (netCDF4.Dataset, xr.Dataset)):
+        raise TypeError("data must be NetCDF4 object or xarray Dataset")
 
     if not (time_index or seconds_run):
         raise ValueError("Input of time_index or seconds_run needed")
@@ -199,8 +199,8 @@ def get_layer_data(
     if not isinstance(layer_index, int):
         raise TypeError("layer_index must be an int")
 
-    if not isinstance(data, netCDF4.Dataset):
-        raise TypeError("data must be NetCDF4 object")
+    if not isinstance(data, (netCDF4.Dataset, xr.Dataset)):
+        raise TypeError("data must be NetCDF4 object or xarray Dataset")
 
     if variable not in data.variables.keys():
         raise ValueError("variable not recognized")
@@ -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 = {
@@ -534,8 +538,10 @@ def variable_interpolation(
                 f"If a string, points must be cells or faces. Got {points}"
             )
 
-    if not isinstance(data, netCDF4.Dataset):
-        raise TypeError(f"data must be netCDF4 object. Got {type(data)}")
+    if not isinstance(data, (netCDF4.Dataset, xr.Dataset)):
+        raise TypeError(
+            f"data must be netCDF4 object or xarray Dataset. Got {type(data)}"
+        )
 
     if not isinstance(to_pandas, bool):
         raise TypeError(f"to_pandas must be of type bool. Got: {type(to_pandas)}")
@@ -616,8 +622,8 @@ def get_all_data_points(
     if not isinstance(time_index, int):
         raise TypeError("time_index must be an int")
 
-    if not isinstance(data, netCDF4.Dataset):
-        raise TypeError("data must be NetCDF4 object")
+    if not isinstance(data, (netCDF4.Dataset, xr.Dataset)):
+        raise TypeError("data must be NetCDF4 object or xarray Dataset")
 
     if variable not in data.variables.keys():
         raise ValueError("variable not recognized")
@@ -637,6 +643,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"][:],
@@ -784,8 +794,8 @@ def turbulent_intensity(
             f"value of the max time index {max_time_index}"
         )
 
-    if not isinstance(data, netCDF4.Dataset):
-        raise TypeError("data must be netCDF4 object")
+    if not isinstance(data, (netCDF4.Dataset, xr.Dataset)):
+        raise TypeError("data must be netCDF4 object or xarray Dataset")
 
     for variable in ["turkin1", "ucx", "ucy", "ucz"]:
         if variable not in data.variables.keys():
@@ -878,11 +888,61 @@ def list_variables(data: Union[netCDF4.Dataset, xr.Dataset, xr.DataArray]) -> Li
     >>> print(variables)
     ['time', 'x', 'y', 'waterdepth', 'ucx', 'ucy', 'ucz', 'turkin1']
     """
-    if isinstance(data, netCDF4.Dataset):
+    if isinstance(
+        data,
+        netCDF4.Dataset,
+    ):
         return list(data.variables.keys())
     if isinstance(data, (xr.Dataset, xr.DataArray)):
         return list(data.variables.keys())
     raise TypeError(
         "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. 
+
+         NASA. (n.d.). Examining spatial (grid) convergence. Accessed Febuary 3, 2026. NASA. 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.
+    factor_of_safety: float
+        Factor of safety (default is 1.25).
+    order: int
+        Order of accuracy (default is 2).
+
+    Returns
+    -------
+    gci: float
+        Grid Convergence Index (GCI).
+    """
+
+    # Validate inputs
+    if not (np.issubdtype(refinement_ratio.dtype, np.number)):
+        raise TypeError("refinement_ratio must be a numeric values")
+    if not (np.issubdtype(factor_of_safety.dtype, np.number)):
+        raise TypeError("factor_of_safety must be a numeric values")
+    if not (np.issubdtype(order.dtype, np.number)):
+        raise TypeError("order must be a numeric values")
+    if not (np.issubdtype(fine_grid.dtype, np.number) and np.issubdtype(coarse_grid.dtype, np.number)):
+        raise TypeError("fine_grid and coarse_grid must contain numeric values")
+    if fine_grid.shape != coarse_grid.shape:
+        raise ValueError("fine_grid and coarse_grid must have the same shape")
+
+    # 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

mhkit/tests/river/test_io_d3d.py

@@ -297,6 +297,15 @@ def test_turbulent_intensity(self):
         ucx_size = np.size(ucx["ucx"])
         self.assertEqual(TI_size, ucx_size)
 
+    def test_calculate_grid_convergence_index(self):
+        fine_grid = np.array([1.0, 2.0, 3.0])
+        coarse_grid = np.array([0.5, 1.5, 2.5])
+        refinement_ratio = 2.0
+        gci = river.io.d3d.calculate_grid_convergence_index(
+            fine_grid, coarse_grid, refinement_ratio
+        )
+        assert_array_almost_equal(gci, np.array([0.2083, 0.1042, 0.0694]), decimal=3)
+
 
 if __name__ == "__main__":
     unittest.main()