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)
No description has been provided for this image No description has been provided for this image
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]: