+
+
+
+
+
+
+In [1]:
+
+
+
+
+
+import io
+import math
+import os
+import pathlib
+
+import holoviews as hv
+import hvplot.pandas
+import numpy as np
+import pandas as pd
+import pymap3d
+import pyposeidon.mesh as pmesh
+
+hv.extension("bokeh")
+np.set_printoptions(suppress=True)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [2]:
+
+
+
+
+
+def parse_hgrid_nodes(path: os.PathLike[str] | str) -> pd.DataFrame:
+ with open(path, "rb") as fd:
+ _ = fd.readline()
+ _, no_points = map(int, fd.readline().strip().split(b" "))
+ content = io.BytesIO(b''.join(next(fd) for _ in range(no_points)))
+ nodes = pd.read_csv(
+ content,
+ engine="pyarrow",
+ sep="\t",
+ header=None,
+ names=["lon", "lat", "depth"],
+ index_col=0
+ )
+ nodes = nodes.reset_index(drop=True)
+ return nodes
+
+def parse_hgrid_elements3(path: os.PathLike[str] | str) -> pd.DataFrame:
+ with open(path, "rb") as fd:
+ _ = fd.readline()
+ no_elements, no_points = map(int, fd.readline().strip().split(b" "))
+ for _ in range(no_points):
+ next(fd)
+ content = io.BytesIO(b''.join(next(fd) for _ in range(no_elements)))
+ elements = pd.read_csv(
+ content,
+ engine="pyarrow",
+ sep="\t",
+ header=None,
+ names=["no_nodes", "n1", "n2", "n3"],
+ index_col=0
+ )
+ elements = elements.assign(
+ n1=elements.n1 - 1,
+ n2=elements.n2 - 1,
+ n3=elements.n3 - 1,
+ ).reset_index(drop=True)
+ return elements
+
+
+
+
+
+
+
+In [3]:
+
+
+
+
+
+def get_skews_and_base_cfls(lons, lats, depths) -> np.ndarray:
+ # The shape of each one of the input arrays needs to be (3, <no_triangles>)
+ #ell = pymap3d.Ellipsoid.from_name("wgs84")
+ ell = pymap3d.Ellipsoid(6378206.4, 6378206.4, "schism", "schism")
+ local_x, local_y, _ = pymap3d.geodetic2enu(lats, lons, depths, lats[0], lons[0], depths[0], ell=ell)
+ areas = (local_x[1] * local_y[2] - local_x[2] * local_y[1]) * 0.5
+ rhos = np.sqrt(areas / np.pi)
+ max_sides = np.maximum(
+ np.sqrt(local_x[1] ** 2 + local_y[1] ** 2),
+ np.sqrt(local_x[2] ** 2 + local_y[2] ** 2),
+ np.sqrt((local_x[2] - local_x[1]) ** 2 + (local_y[2] - local_y[1]) ** 2),
+ )
+ skews = max_sides / rhos
+ base_cfls = np.sqrt(9.81 * np.maximum(0.1, depths.mean(axis=0))) / rhos / 2
+ return skews, base_cfls
+
+def get_skews_and_base_cfls_from_path(path: os.PathLike[str] | str) -> np.ndarray:
+ nodes = parse_hgrid_nodes(path)
+ elements = parse_hgrid_elements3(path)
+ tri = elements[["n1", "n2", "n3"]].values
+ lons = nodes.lon.values[tri].T
+ lats = nodes.lat.values[tri].T
+ depths = nodes.depth.values[tri].T
+ skews, base_cfls = get_skews_and_base_cfls(lons=lons, lats=lats, depths=depths)
+ return skews, base_cfls
+
+
+
+
+
+
+
+In [4]:
+
+
+
+
+
+path = "/home/panos/Prog/poseidon/seareport_meshes/meshes/global-v0.1.gr3"
+path = "/home/panos/Prog/poseidon/seareport_meshes/meshes/global-v0.gr3"
+path = "/home/panos/Prog/git/schism/src/Utility/Grid_Scripts/hgrid.gr3"
+path = "/home/panos/Prog/poseidon/seareport_meshes/meshes/global-v0.2.gr3"
+skews, base_cfls = get_skews_and_base_cfls_from_path(path)
+
+
+
+
+
+
+
+In [5]:
+
+
+
+
+
+CFL_THRESHOLD = 0.4
+for dt in (1, 50, 75, 100, 120, 150, 200, 300, 400, 600, 900, 1200, 1800, 3600):
+ violations = (base_cfls * dt < CFL_THRESHOLD).sum()
+ print(f"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%")
+
+
+
+
+
+
+
+
+
+
+
+1 716169 100.00% + 50 438858 61.28% + 75 173842 24.27% + 100 122625 17.12% + 120 107230 14.97% + 150 87782 12.26% + 200 67136 9.37% + 300 39294 5.49% + 400 21114 2.95% + 600 564 0.08% + 900 0 0.00% +1200 0 0.00% +1800 0 0.00% +3600 0 0.00% ++
+
+
+
+
+
+
+In [6]:
+
+
+
+
+
+pd.DataFrame({"skew": skews}).describe()
+
+
+
+
+
+
+
+
+Out[6]:
+
+
+
+
+
+
+
+
+| + | skew | +
|---|---|
| count | +716169.000000 | +
| mean | +2.819506 | +
| std | +0.144765 | +
| min | +2.507587 | +
| 25% | +2.721086 | +
| 50% | +2.780963 | +
| 75% | +2.885701 | +
| max | +6.666141 | +
+
+
+
+
+
+
+In [7]:
+
+
+
+
+
+df = pd.DataFrame({"cfl": base_cfls * 400})
+df[df.cfl < 0.4].describe()
+df[df.cfl < 0.4].hvplot.hist()
+
+
+
+
+
+
+
+
+Out[7]:
+
+
+
+
+
+
+
+
+| + | cfl | +
|---|---|
| count | +21114.000000 | +
| mean | +0.330807 | +
| std | +0.036809 | +
| min | +0.199615 | +
| 25% | +0.302369 | +
| 50% | +0.327259 | +
| 75% | +0.360951 | +
| max | +0.400000 | +
+
+
+
+
+Out[7]:
+
+
+
+
+
+
+
+
+
+
+
+
+
+In [8]:
+
+
+
+
+
+nodes = parse_hgrid_nodes(path)
+elements = parse_hgrid_elements3(path)
+elements = elements.assign(base_cfl=base_cfls)
+elements.head()
+
+
+
+
+
+
+
+
+Out[8]:
+
+
+
+
+
+
+
+
+| + | no_nodes | +n1 | +n2 | +n3 | +base_cfl | +
|---|---|---|---|---|---|
| 0 | +3 | +189444 | +0 | +244205 | +0.003675 | +
| 1 | +3 | +0 | +189444 | +335357 | +0.004649 | +
| 2 | +3 | +0 | +207657 | +218301 | +0.007291 | +
| 3 | +3 | +207657 | +0 | +319429 | +0.006633 | +
| 4 | +3 | +0 | +218301 | +361609 | +0.006373 | +
+
+
+
+
+
+
+In [9]:
+
+
+
+
+
+min_cfl_per_node = pd.concat([
+ elements[["n1", "base_cfl"]].groupby(["n1"]).base_cfl.min(),
+ elements[["n2", "base_cfl"]].groupby(["n2"]).base_cfl.min(),
+ elements[["n3", "base_cfl"]].groupby(["n3"]).base_cfl.min(),
+], axis=1).min(axis=1)
+min_cfl_per_node.head()
+
+
+
+
+
+
+
+
+Out[9]:
+
+
+0 0.003675 +1 0.007158 +2 0.009419 +3 0.002360 +4 0.007654 +dtype: float64+
+
+
+
+
+
+
+In [15]:
+
+
+
+
+
+dt = 600
+df = nodes.assign(
+ cfl=min_cfl_per_node * dt,
+ # CFL_violation nodes have a value of 1 if there is no violation and 4 if there is a violation.
+ # We do this in order to plot the points with a different size
+ cfl_violation=((min_cfl_per_node * dt < CFL_THRESHOLD) * 3) + 1
+)
+df.head()
+
+
+
+
+
+
+
+
+Out[15]:
+
+
+
+
+
+
+
+
+| + | lon | +lat | +depth | +cfl | +cfl_violation | +
|---|---|---|---|---|---|
| 0 | +0.000000 | +90.000000 | +4228.0 | +2.205227 | +1 | +
| 1 | +-23.447414 | +-19.913775 | +4603.0 | +4.294965 | +1 | +
| 2 | +118.483781 | +-15.726201 | +5477.0 | +5.651691 | +1 | +
| 3 | +15.495643 | +-28.939963 | +187.0 | +1.416196 | +1 | +
| 4 | +94.479055 | +-14.068985 | +5710.0 | +4.592690 | +1 | +
+
+
+
+
+
+
+In [16]:
+
+
+
+
+
+df.hvplot.points(
+ 'lon',
+ 'lat',
+ c="cfl_violation",
+ cmap="colorblind",
+ geo=True,
+ tiles=True,
+).options(
+ width=900, height=600
+).opts(
+ hv.opts.Points(size=hv.dim('cfl_violation'))
+)
+
+
+
+
+
+
+
+
+
Out[16]:
+
+
+
+
+