diff --git a/mhkit/river/io/d3d.py b/mhkit/river/io/d3d.py index 7295d7e11..8a760ee78 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) @@ -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,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()) @@ -886,3 +899,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 diff --git a/mhkit/tests/river/test_io_d3d.py b/mhkit/tests/river/test_io_d3d.py index f41ba4962..e2e159082 100644 --- a/mhkit/tests/river/test_io_d3d.py +++ b/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()