Mesh density as a function of CFL?

[ivicajan]

Dear all,
if anyone is using SCHISM running in 3D baroclinic mode, then you probably know about desired CFL range (basically not too small and not too big i.e. range of CFL between 0.6 and 6) for a given time step.

Is there an option or plan to have CFL optimised grid? It would be an additional constraint for a given time step.

Thanks in advance,
Ivica

[SorooshMani-NOAA]

@ivicajan there's this function to limit the Courant number in the mesh-size function, based on some approximations.

OCSMesh/ocsmesh/hfun/raster.py

Lines 684 to 690 in 18d15ac

	def add_courant_num_constraint(
	self,
	upper_bound: float = 0.9,
	lower_bound: Optional[float] = None,
	timestep: float = 150,
	wave_amplitude: float = 2
	) -> None:

Is this what you're looking for or do you want more control over any aspect of it?

[ivicajan]

Thanks for a swift reply!

I am looking for something to tweak my old grid, and not sure if it is possible. In your post it is defined on the raster only? do you have a few lines example how to add the CFL lower/upper constrain?
For example:
mesh_obj = ocsmesh.Mesh.open('hgrid.gr3', crs=4326)
geom_obj = ocsmesh.Geom(mesh_obj_2)

Ultimate goal is to cleanup and optimize the grid for different time step (i.e. make it work for larger time steps -> coarsen the grid).

Thanks in advance.

[SorooshMani-NOAA]

Modifying mesh is not OCSMesh's strongest suit, but you can definitely remesh while specifying initial points (your old mesh nodes). For remesh I haven't yet implemented a clean API in OCSMesh, but you can still use OCSMesh to prepare inputs for Jigsaw (the mesh engine OCSMesh uses).

The scripts below adds a CFL constraint on the mesh size, but doesn't enforce old edges and vertices to be fixed. It just approximates the old mesh size where the CFL is already met.

from copy import deepcopy

import numpy as np
from ocsmesh import Mesh, Hfun, Geom, utils, JigsawDriver

# import other stuff...


old_mesh_crs = ...
desired_timestep = ...
high_cfl_value = ...
low_cfl_value = ...
nu = ...

old_mesh = Mesh.open('/path/to/your/old/mesh', crs=old_mesh_crs)

# Get topography from the old mesh (assuming it's the best available)
topo = old_mesh.msh_t.value


# Use old mesh's edge sizes for the new hfun
hfun = Hfun(old_mesh)
hfun.size_from_mesh()

# Apply CFL constraints
hfun_msht = deepcopy(hfun.msh_t())
vel = utils.estimate_particle_velocity_from_depth(
    depth=topo, wave_amplitude=nu
)
size_cfl_hi = utils.get_element_size_courant(
    characteristic_velocity_magnitude=vel
    timestep=desired_timestep
    target_courant=high_cfl_value
)
size_cfl_lo = utils.get_element_size_courant(
    characteristic_velocity_magnitude=vel
    timestep=desired_timestep
    target_courant=low_cfl_value
)
## Sanity check (at each location hi Courant # is larger than lo)
assert np.all(size_cfl_hi > size_cfl_lo)

## What is the Courant number for current mesh sizes
old_cfl = utils.approximate_courant_number_for_depth(
    depth=topo,
    timestep=desired_timestep,
    element_size=hfun_msht.value, # Old element size from mesh
    wave_amplitude=nu
)

## Adjust the size values where it falls out of bound
hi_mask = size_cfl_hi < old_cfl
hfun_msht.value[hi_mask] = size_cfl_hi[hi_mask]

lo_mask = size_cfl_lo > old_cfl
hfun_msht.value[lo_mask] = size_cfl_lo[lo_mask]


geom = Geom(utils.get_mesh_polygons(old_mesh.msh_t))
adj_hfun = Hfun(hfun_msht)

driver = JigsawDriver(geom, adj_hfun)
new_mesh = driver.run()

new_mesh.write('output/path/file.2dm', format='2dm', overwrite=True)

Please let me know if you find this useful.

[ivicajan]

Thanks,
It helped a lot.

I was getting my head around sizes (msh_t.values) and couldn't get it right unitl I figured that vel (in your example is not barotropic wave speed that is used to define CFL). When I added that part all went in the right direction.

hfun_msht = deepcopy(hfun_msht_old_wgs84)
vel = utils.estimate_particle_velocity_from_depth(
depth=topo, wave_amplitude=nu
)
size_cfl_hi = utils.get_element_size_courant(
characteristic_velocity_magnitude=vel + np.sqrt(-topo*9.81),
timestep=desired_timestep,
target_courant=high_cfl_value
)
size_cfl_lo = utils.get_element_size_courant(
characteristic_velocity_magnitude=vel + np.sqrt(-topo*9.81),
timestep=desired_timestep,
target_courant=low_cfl_value
)

Then I had to apply the new sizes:
bad_hi = np.where(old_cfl>=high_cfl_value)[0][:]
bad_lo = np.where(old_cfl<=low_cfl_value)[0][:]

and then:
hfun_msht.value[bad_hi] = size_cfl_hi[bad_hi]
hfun_msht.value[bad_lo] = size_cfl_lo[bad_lo]

When specifying Geom there is a requirement to add csr, something along the lines:
geom = Geom(utils.get_mesh_polygons(hfun_msht_old), crs = CRS.from_user_input('epsg:4326'))

In the end, I didn't know how to interpolate depths from the old mesh onto new.
Then I ended writing a small part to interpolate from old grid to the new (there must be a better way).
If you want I can send you the whole notebook.

Thanks once again !
Ivica

[SorooshMani-NOAA]

No problem, I'm glad that it almost worked. I didn't really test the script, I was just trying to put together something as a guide. Also sorry that I missed some calculations for velocity from

OCSMesh/ocsmesh/features/constraint.py

Lines 296 to 304 in 18d15ac

	u_mag = utils.estimate_particle_velocity_from_depth(
	ref_values, self._nu
	)
	depth_mask = utils.can_velocity_be_approximated_by_linear_wave_theory(
	ref_values, self._nu
	)
	char_vel = u_mag + np.sqrt(constants.g * np.abs(ref_values))
	# For overland where h < nu the characteristic velocity is 2 * sqrt(g*h)
	char_vel[~depth_mask] = 2 * u_mag[~depth_mask]

where the same thing is applied for raster data.

For interpolation, when you have a mesh you can use

OCSMesh/ocsmesh/utils.py

Lines 1217 to 1226 in 18d15ac

	def interpolate(src: jigsaw_msh_t, dst: jigsaw_msh_t, **kwargs):
	if src.mshID == 'euclidean-grid' and dst.mshID == 'euclidean-mesh':
	interpolate_euclidean_grid_to_euclidean_mesh(src, dst, **kwargs)
	elif src.mshID == 'euclidean-mesh' and dst.mshID == 'euclidean-mesh':
	interpolate_euclidean_mesh_to_euclidean_mesh(src, dst, **kwargs)
	else:
	raise NotImplementedError(
	f'Not implemented type combination: source={src.mshID}, '
	f'dest={dst.mshID}')

You need to get the new_mesh.msh_t and pass that to the interpolation function.

[ivicajan]

Possibly you have it in your fingers and I would need ages to figure, how to add shape file (such as coastline - poly lines) and define depth for them to be say -1m and add that to the current raster? In my case I have DEM that is quite coarse and adding "poly lines" (as shape file) with some predefined depths would help a lot. Something along:
coast = gpd.read_file('coastline.shp')
hfun.add_feature(coast.unary_union,
expansion_rate = 0.01,
target_size = 100,
depth = -1,
)

[SorooshMani-NOAA]

My understanding is that you'd like to adjust the DEM data before using it for size function definition, right? If that's the case then take a look at:

OCSMesh/ocsmesh/raster.py

Lines 1221 to 1227 in 18d15ac

	def adjust(
	self,
	geom: Union[None, Polygon, MultiPolygon] = None,
	inside_min: float = -np.inf,
	outside_max: float = np.inf,
	cond: Optional[Callable[[npt.NDArray[float]], npt.NDArray[bool]]] = None
	) -> None:

This requires an input (Multi)Polygon. If you just want to add a line the distance from which defines the size value, then you can use the add_feature or add_patch of the size function.

Does that answer your question?

[ivicajan]

Thanks
it is along the line you wrote, need to fix DEM.
The reason for that is I have bad channels inside the mesh.
My DEM is not that good to idenfiy channels (as in your manual using something similar to: hfun.add_channel(level=0, width=1500, target_size=200, expansion_rate=0.01))
I would rather use shape file as a input to the channel detection rather than raster.

Cheers,
Ivica

[SorooshMani-NOAA]

I see. So maybe take a look at Example 4 in
https://github.com/SorooshMani-NOAA/OCSMeshTutorial/blob/main/UFS_Nov_2022/ocsmesh_basics.ipynb
as well. I guess this is similar to what you want to do for channels.