In [59]:
import io
import math
import os
import pathlib
import holoviews as hv
import hvplot.pandas
import matplotlib.pyplot as plt
import numpy as np
import panel as pn
import pandas as pd
import pymap3d
import xarray as xr
hv.extension("bokeh")
np.set_printoptions(suppress=True)
In [2]:
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:
ds = xr.open_dataset(path, engine='selafin')
tri = ds.attrs['ikle2'] - 1
lons = ds.x.values[tri].T
lats = ds.y.values[tri].T
depths = - ds.B.isel(time=0).values[tri].T
skews, base_cfls = get_skews_and_base_cfls(lons=lons, lats=lats, depths=depths)
return skews, base_cfls
In [109]:
files = {"v1.2" : ["/home/tomsail/Documents/work/models/meshes/slf/v1.2.slf", 50],
"v2.2" : ["/home/tomsail/Documents/work/models/meshes/slf/v2.2.slf", 30],
"v1.5" : ["/home/tomsail/Documents/work/models/meshes/slf/v1.5.slf", 50],
"v2.3" : ["/home/tomsail/Documents/work/models/meshes/slf/v2.3.slf", 30],
"v3.1" : ["/home/tomsail/Documents/work/models/meshes/slf/v3.1.slf", 20]}
plots = []
for v in files.keys():
file = files[v][0]
ideal_dt = files[v][1]
ds = xr.open_dataset(file, engine='selafin')
skews, base_cfls = get_skews_and_base_cfls_from_path(file)
CFL_THRESHOLD = 1
print(v)
for dt in (1, 10, 20, 30, 50, 75, 100, 120, 150, 200, 300, 400):
violations = (base_cfls * dt > CFL_THRESHOLD).sum()
print(f"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%")
_skews = pd.DataFrame({"skew": skews}).hvplot.hist(bins=40,bin_range=[2.5,3.5]).opts(title = v + " skewness")
_cfls = pd.DataFrame({"cfls": base_cfls * ideal_dt}).hvplot.hist(bins=40, bin_range = [0,1.5]).opts(title = v + " CFL")
both = (_skews + _cfls)
plots.append(both)
v1.2 1 0 0.00% 10 0 0.00% 20 0 0.00% 30 26970 0.85% 50 698959 22.16% 75 1522191 48.25% 100 1988032 63.02% 120 2200208 69.74% 150 2375695 75.31% 200 2529722 80.19% 300 2706719 85.80% 400 2960200 93.83% v2.2 1 0 0.00% 10 1 0.00% 20 101192 1.72% 30 795500 13.51% 50 2272648 38.60% 75 3343353 56.79% 100 3958923 67.25% 120 4261793 72.39% 150 4550877 77.30% 200 4856354 82.49% 300 5222952 88.72% 400 5625413 95.55% v1.5 1 0 0.00% 10 0 0.00% 20 23 0.00% 30 24176 0.74% 50 761333 23.44% 75 1877009 57.80% 100 2345059 72.21% 120 2492877 76.76% 150 2604729 80.20% 200 2712201 83.51% 300 2872349 88.45% 400 2994566 92.21% v2.3 1 0 0.00% 10 4 0.00% 20 20468 0.37% 30 291651 5.34% 50 1468831 26.88% 75 2822569 51.65% 100 3460263 63.32% 120 3745588 68.55% 150 4058524 74.27% 200 4451012 81.46% 300 4972399 91.00% 400 5257787 96.22% v3.1 1 0 0.00% 10 408 0.00% 20 454600 4.56% 30 1269435 12.73% 50 3102337 31.12% 75 5070873 50.86% 100 6291124 63.10% 120 6977276 69.99% 150 7785133 78.09% 200 8892359 89.19% 300 9703458 97.33% 400 9893952 99.24%
In [110]:
hv.Layout(plots).cols(2)
Out[110]:
In [75]:
tri = ds.attrs['ikle2'] - 1
nodes = pd.DataFrame(np.vstack((ds.x, ds.y, ds.B.isel(time=0))).T, columns=["lon", "lat", "depth"])
elements = pd.DataFrame(np.vstack( (np.ones(len(tri))* 3, tri.T)).T , columns=["no_nodes", "n1", "n2", "n3"])
elements = elements.assign(base_cfl=base_cfls)
elements.head()
Out[75]:
no_nodes | n1 | n2 | n3 | base_cfl | |
---|---|---|---|---|---|
0 | 3.0 | 65536.0 | 575051.0 | 1405784.0 | 0.014614 |
1 | 3.0 | 575051.0 | 65536.0 | 968366.0 | 0.017156 |
2 | 3.0 | 221589.0 | 65536.0 | 1405784.0 | 0.014376 |
3 | 3.0 | 65536.0 | 221589.0 | 1584364.0 | 0.014610 |
4 | 3.0 | 65536.0 | 499883.0 | 968366.0 | 0.017411 |
In [76]:
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[76]:
1.0 0.019817 2.0 0.008588 3.0 0.014766 4.0 0.007269 5.0 0.002010 dtype: float32
In [86]:
dt = 30
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[86]:
lon | lat | depth | cfl | cfl_violation | |
---|---|---|---|---|---|
0 | -180.000000 | 90.000000 | -4228.0 | NaN | NaN |
1 | 85.697014 | -32.376286 | -3436.0 | 0.594523 | 1.0 |
2 | -57.604584 | -45.496311 | -3668.0 | 0.257643 | 1.0 |
3 | -90.156075 | 27.119200 | -2249.0 | 0.442983 | 1.0 |
4 | -139.083511 | -48.631084 | -4250.0 | 0.218060 | 1.0 |
In [87]:
plot = df[df.cfl_violation == 4].hvplot.points(
'lon',
'lat',
c="depth",
cmap="jet",
geo=True,
tiles="EsriImagery",
).options(
width=1200, height=900
)
len(df[df.cfl_violation == 4])
Out[87]:
2600
In [88]:
plot
Out[88]:
In [107]:
plots = []
for v in files.keys():
file = files[v][0]
ideal_dt = files[v][1]
ds = xr.open_dataset(file, engine='selafin')
skews, base_cfls = get_skews_and_base_cfls_from_path(file)
CFL_THRESHOLD = 1
tri = ds.attrs['ikle2'] - 1
#
nodes = pd.DataFrame(np.vstack((ds.x, ds.y, ds.B.isel(time=0))).T, columns=["lon", "lat", "depth"])
elements = pd.DataFrame(np.vstack( (np.ones(len(tri))* 3, tri.T)).T , columns=["no_nodes", "n1", "n2", "n3"])
elements = elements.assign(cfls=base_cfls*ideal_dt, skews = skews)
min_cfl_per_node = pd.concat([
elements[["n1", "cfls"]].groupby(["n1"]).cfls.min(),
elements[["n2", "cfls"]].groupby(["n2"]).cfls.min(),
elements[["n3", "cfls"]].groupby(["n3"]).cfls.min(),
], axis=1).min(axis=1)
mean_skew_per_node = pd.concat([
elements[["n1", "skews"]].groupby(["n1"]).skews.mean(),
elements[["n2", "skews"]].groupby(["n2"]).skews.mean(),
elements[["n3", "skews"]].groupby(["n3"]).skews.mean(),
], axis=1).mean(axis=1)
df = nodes.assign(
cfl=min_cfl_per_node,
# 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 > CFL_THRESHOLD) * 3) + 1
)
plot = df[df.cfl_violation == 4].hvplot.points(
'lon',
'lat',
c="depth",
cmap="jet",
geo=True,
tiles="EsriImagery",
).options(
width=1200, height=900, title = v
)
print(v, len(df[df.cfl_violation == 4]))
plots.append(plot)
v1.2 251006 v2.2 275439 v1.5 293080 v2.3 102597 v3.1 171305
In [108]:
hv.Layout(plots).cols(1)
Out[108]: