From bbbdf679ce00c537b2072dca6401909b1df5fe11 Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Tue, 18 Nov 2025 13:29:47 -0900 Subject: [PATCH 1/7] begin to add xarray capability and add new variables in coords --- mhkit/river/io/d3d.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/mhkit/river/io/d3d.py b/mhkit/river/io/d3d.py index 7295d7e11..b4c854264 100644 --- a/mhkit/river/io/d3d.py +++ b/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"][:], @@ -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"][:], From bb27980efaa51deaed38c850db6b0186ba3cea21 Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Wed, 19 Nov 2025 10:18:15 -0900 Subject: [PATCH 2/7] fix waterdepth name --- mhkit/river/io/d3d.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mhkit/river/io/d3d.py b/mhkit/river/io/d3d.py index b4c854264..f0491e69d 100644 --- a/mhkit/river/io/d3d.py +++ b/mhkit/river/io/d3d.py @@ -257,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 = { From 2e78821536e87a3f317df1e05653df4c39734b52 Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Wed, 28 Jan 2026 09:28:58 -0900 Subject: [PATCH 3/7] calculate_grid_convergence_index --- mhkit/river/io/d3d.py | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/mhkit/river/io/d3d.py b/mhkit/river/io/d3d.py index f0491e69d..79b0399e5 100644 --- a/mhkit/river/io/d3d.py +++ b/mhkit/river/io/d3d.py @@ -894,3 +894,30 @@ 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 \ No newline at end of file From a43446ca94456830054d655fe5c5bdcd5ae6ace4 Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Wed, 28 Jan 2026 10:05:15 -0900 Subject: [PATCH 4/7] add test for calculate_grid_convergence_index --- mhkit/tests/river/test_io_d3d.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/mhkit/tests/river/test_io_d3d.py b/mhkit/tests/river/test_io_d3d.py index f41ba4962..86caf8bd1 100644 --- a/mhkit/tests/river/test_io_d3d.py +++ b/mhkit/tests/river/test_io_d3d.py @@ -296,7 +296,15 @@ def test_turbulent_intensity(self): ucx = river.io.d3d.get_all_data_points(data, "ucx", time_index) 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 + ) + self.assertAlmostEqual(gci, np.array([0.2083, 0.1042, 0.0694])) if __name__ == "__main__": unittest.main() From 0f5e015bf48707854c674818a2d8932dd6fdc547 Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Wed, 28 Jan 2026 10:28:18 -0900 Subject: [PATCH 5/7] black an pylint formating --- mhkit/river/io/d3d.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/mhkit/river/io/d3d.py b/mhkit/river/io/d3d.py index 79b0399e5..847795a82 100644 --- a/mhkit/river/io/d3d.py +++ b/mhkit/river/io/d3d.py @@ -895,7 +895,10 @@ def list_variables(data: Union[netCDF4.Dataset, xr.Dataset, xr.DataArray]) -> Li f"xarray DataArray. Got: {type(data)}" ) -def calculate_grid_convergence_index(fine_grid, coarse_grid, refinement_ratio,factor_of_safety=1.25, order=2): + +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 @@ -905,7 +908,7 @@ def calculate_grid_convergence_index(fine_grid, coarse_grid, refinement_ratio,fa Results from the finer grid. coarse_grid: numpy.ndarray Results from the coarser grid. - refinement_ratio: float + refinement_ratio: float Refinement ratio between the grids. order: int Order of accuracy (default is 2). @@ -920,4 +923,4 @@ def calculate_grid_convergence_index(fine_grid, coarse_grid, refinement_ratio,fa # Calculate the GCI gci = (factor_of_safety * error) / (refinement_ratio**order - 1) - return gci \ No newline at end of file + return gci From d2155a8be08e7776767db6293f790f3f8abf6433 Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Wed, 28 Jan 2026 10:28:39 -0900 Subject: [PATCH 6/7] black and pulint formatting --- mhkit/tests/river/test_io_d3d.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/mhkit/tests/river/test_io_d3d.py b/mhkit/tests/river/test_io_d3d.py index 86caf8bd1..e2e159082 100644 --- a/mhkit/tests/river/test_io_d3d.py +++ b/mhkit/tests/river/test_io_d3d.py @@ -296,7 +296,7 @@ def test_turbulent_intensity(self): ucx = river.io.d3d.get_all_data_points(data, "ucx", time_index) 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]) @@ -304,7 +304,8 @@ def test_calculate_grid_convergence_index(self): gci = river.io.d3d.calculate_grid_convergence_index( fine_grid, coarse_grid, refinement_ratio ) - self.assertAlmostEqual(gci, np.array([0.2083, 0.1042, 0.0694])) + assert_array_almost_equal(gci, np.array([0.2083, 0.1042, 0.0694]), decimal=3) + if __name__ == "__main__": unittest.main() From 3d115173f213212ed9f779b50bcbaa216db205eb Mon Sep 17 00:00:00 2001 From: Emily Browning Date: Wed, 28 Jan 2026 10:38:58 -0900 Subject: [PATCH 7/7] add xarray input to error checks --- mhkit/river/io/d3d.py | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/mhkit/river/io/d3d.py b/mhkit/river/io/d3d.py index 847795a82..8a760ee78 100644 --- a/mhkit/river/io/d3d.py +++ b/mhkit/river/io/d3d.py @@ -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") @@ -538,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)}") @@ -620,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") @@ -792,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(): @@ -886,7 +888,10 @@ 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())