Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

max_courant_time seems to have no effect on refinement #73

Closed
veenstrajelmer opened this issue Jul 5, 2023 · 1 comment · Fixed by Deltares/MeshKernel#226
Closed
Labels
bug Something isn't working

Comments

@veenstrajelmer
Copy link
Collaborator

Describe the bug
When refining a grid based on bathymetry samples, there is a parameter called max_courant_time. I would expect it to do something similar as min_edge_size, but instead it does not seem to do anything.

To Reproduce

import meshkernel
import xarray as xr
import matplotlib.pyplot as plt
plt.close('all')

#general settings
lon_min,lon_max = -6,2
lat_min,lat_max = 48.5,51.2
lon_res,lat_res = 0.2,0.2
is_geographic = True

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters(origin_x=lon_min,
                                                     origin_y=lat_min,
                                                     upper_right_x=lon_max,
                                                     upper_right_y=lat_max,
                                                     block_size_x=lon_res,
                                                     block_size_y=lat_res)
mk = meshkernel.MeshKernel(is_geographic=is_geographic)
mk.curvilinear_make_uniform_on_extension(make_grid_parameters)
mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d


#select bathy data and convert to GriddedSamples
data_bathy = xr.open_dataset(r'p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc')
data_bathy_sel = data_bathy.sel(lon=slice(lon_min-1,lon_max+1),lat=slice(lat_min-1,lat_max+1))
lon_np = data_bathy_sel.lon.to_numpy()
lat_np = data_bathy_sel.lat.to_numpy()
values_np = data_bathy_sel.elevation.to_numpy().flatten().astype('float')
gridded_samples = meshkernel.GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np)


#refinement
mesh_refinement_parameters = meshkernel.MeshRefinementParameters(min_edge_size=5000, #always in meters
                                                                 refinement_type=meshkernel.RefinementType(1), #Wavecourant/1,
                                                                 connect_hanging_nodes=False,
                                                                 smoothing_iterations=0,
                                                                 max_courant_time=200) #TODO: altering this does not seem to have effect

mk.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
                                          mesh_refinement_params=mesh_refinement_parameters)

mesh2d_refinedgrid = mk.mesh2d_get()
fig, ax = plt.subplots()
mesh2d_refinedgrid.plot_edges(ax,linewidth=0.8)

Expected behavior
A different result when supplying different values for max_courant_time

Version info (please complete the following information):

@veenstrajelmer
Copy link
Collaborator Author

When using this code with the meshkernel 3.0.0 pre-release, the issue is solved. max_courant_time 100 and 200 result in different grids:

import meshkernel
import xarray as xr
import matplotlib.pyplot as plt
plt.close('all')

#general settings
lon_min,lon_max = -6,2
lat_min,lat_max = 48.5,51.2
lon_res,lat_res = 0.2,0.2

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters(origin_x=lon_min,
                                                     origin_y=lat_min,
                                                     upper_right_x=lon_max,
                                                     upper_right_y=lat_max,
                                                     block_size_x=lon_res,
                                                     block_size_y=lat_res)
mk = meshkernel.MeshKernel(projection=True)
mk.get_projection()
mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d


#select bathy data and convert to GriddedSamples
data_bathy = xr.open_dataset(r'p:\metocean-data\open\GEBCO\2021\GEBCO_2021.nc')
data_bathy_sel = data_bathy.sel(lon=slice(lon_min-1,lon_max+1),lat=slice(lat_min-1,lat_max+1))
lon_np = data_bathy_sel.lon.to_numpy()
lat_np = data_bathy_sel.lat.to_numpy()
values_np = data_bathy_sel.elevation.to_numpy().flatten().astype('float')
gridded_samples = meshkernel.GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np)


#refinement
mesh_refinement_parameters = meshkernel.MeshRefinementParameters(min_edge_size=1000, #always in meters
                                                                 refinement_type=meshkernel.RefinementType(1), #Wavecourant/1,
                                                                 connect_hanging_nodes=False,
                                                                 smoothing_iterations=0,
                                                                 max_courant_time=100)

mk.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
                                          mesh_refinement_params=mesh_refinement_parameters)

mesh2d_refinedgrid = mk.mesh2d_get()
fig, ax = plt.subplots()
mesh2d_refinedgrid.plot_edges(ax,linewidth=0.8)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

Successfully merging a pull request may close this issue.

1 participant