CHLNDDEV / oceanmesh

Automatic coastal ocean mesh generation in Python and C++. **under development**
GNU General Public License v3.0
51 stars 15 forks source link

offset in sizing functions using dem (WL and bathy gradient) #70

Closed tomsail closed 6 months ago

tomsail commented 9 months ago

The @setter for grid.values is misleading during the calculation of the sizing functions that require a dem.

The values grid.dx and grid.dy change here for the bathy gradient and here for the wavelength

Here an example for the world model where both the wavelength and bathy gradient have an offset: om_bug

The quick fix (without touching the Grid class) is to force in both functions the dx and dy in the sizing functions.

grid.dx = dem.dx
grid.dy = dem.dy

Here is the result after using the fix: image

tomsail commented 9 months ago

to replicate :

# Note: global_stereo.shp has been generated using global_tag in pyposeidon
# https://github.com/ec-jrc/pyPoseidon/blob/9cfd3bbf5598c810004def83b1f43dc5149addd0/pyposeidon/boundary.py#L452
fname = os.path.join(os.path.dirname(__file__), "global", "global_latlon.shp")
fname2 = os.path.join(os.path.dirname(__file__), "global", "global_stereo.shp")
fnameDEM = os.path.join(os.path.dirname(__file__), "global", "dem.nc")

# it is necessary to define all the coastlines at once:
# the Shoreline class will the detect the biggest coastline (antartica and define it
# as the outside boundary)

EPSG = 4326  # EPSG:4326 or WGS84
bbox = (-180.00, 180.00, -89.00, 90.00)
extent = om.Region(extent=bbox, crs=4326)

min_edge_length = 0.5  # minimum mesh size in domain in meters
max_edge_length = 2  # maximum mesh size in domain in meters
shoreline = om.Shoreline(fname, extent.bbox, min_edge_length)
sdf = om.signed_distance_function(shoreline)
edge_length0 = om.distance_sizing_function(shoreline, rate=0.11)
edge_length1 = om.feature_sizing_function(
    shoreline,
    sdf,
    min_edge_length=min_edge_length,
    max_edge_length=max_edge_length,
    crs=EPSG,
)
dem = om.DEM(fnameDEM, bbox=extent)
edge_length2 = om.wavelength_sizing_function(
    dem,
    wl=100,
    period=12.42 * 3600,
    min_edgelength=min_edge_length,
    max_edge_length=max_edge_length,
    crs=EPSG,
)  # use the M2-tide period (in seconds)
edge_length3 = om.bathymetric_gradient_sizing_function(
    dem,
    slope_parameter=30.0,
    min_edge_length=min_edge_length,
    max_edge_length=max_edge_length,
    crs=EPSG,
)
edge_length = om.compute_minimum(
    [edge_length0, edge_length1, edge_length2, edge_length3]
)
edge_length = om.enforce_mesh_gradation(edge_length, gradation=0.11, stereo=True)

# once the size functions have been defined, wed need to mesh inside domain in
# stereographic projections. This is way we use another coastline which is
# already translated in a sterographic projection
shoreline_stereo = om.Shoreline(fname2, extent.bbox, min_edge_length, stereo=True)
domain = om.signed_distance_function(shoreline_stereo)

points, cells = om.generate_mesh(domain, edge_length, stereo=True, max_iter=100)

# remove degenerate mesh faces and other common problems in the mesh
points, cells = om.make_mesh_boundaries_traversable(points, cells)
points, cells = om.delete_faces_connected_to_one_face(points, cells)

# apply a Laplacian smoother
points, cells = om.laplacian2(points, cells, max_iter=100)

# plot
fig, ax, _ = edge_length.plot(
    holding=True,
    plot_colorbar=True,
    stereo=True,
    vmax=max_edge_length,
)

ax.triplot(points[:, 0], points[:, 1], cells, color="gray", linewidth=0.5)
shoreline_stereo.plot(ax=ax)

with dem.nc obtained from the data folder:

curl -L url https://www.dropbox.com/sh/nd2b012wrpt6qnq/AAAD7aA_qXztUhlT39YK2yBua?dl=1 > data.zip then unzip as unzip data.zip -d data/

from pyposeidon's Tutorial page

tomsail commented 6 months ago

this has been fixed with #71