CFL-TELEMAC.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import io\n",
    "import math\n",
    "import os\n",
    "import pathlib\n",
    "\n",
    "import holoviews as hv\n",
    "import hvplot.pandas\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import panel as pn\n",
    "import pandas as pd\n",
    "import pymap3d\n",
    "import xarray as xr\n",
    "\n",
    "hv.extension(\"bokeh\")\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_skews_and_base_cfls(lons, lats, depths) -> np.ndarray:\n",
    "    # The shape of each one of the input arrays needs to be (3, <no_triangles>)\n",
    "    #ell = pymap3d.Ellipsoid.from_name(\"wgs84\")\n",
    "    ell = pymap3d.Ellipsoid(6378206.4, 6378206.4, \"schism\", \"schism\")\n",
    "    local_x, local_y, _ = pymap3d.geodetic2enu(lats, lons, depths, lats[0], lons[0], depths[0], ell=ell)\n",
    "    areas = (local_x[1] * local_y[2] - local_x[2] * local_y[1]) * 0.5\n",
    "    rhos = np.sqrt(areas / np.pi)\n",
    "    max_sides = np.maximum(\n",
    "        np.sqrt(local_x[1] ** 2 + local_y[1] ** 2),\n",
    "        np.sqrt(local_x[2] ** 2 + local_y[2] ** 2),\n",
    "        np.sqrt((local_x[2] - local_x[1]) ** 2 + (local_y[2] - local_y[1]) ** 2),\n",
    "    )\n",
    "    skews = max_sides / rhos\n",
    "    base_cfls = np.sqrt(9.81 * np.maximum(0.1, depths.mean(axis=0))) / rhos / 2\n",
    "    return skews, base_cfls\n",
    "\n",
    "def get_skews_and_base_cfls_from_path(path: os.PathLike[str] | str) -> np.ndarray:\n",
    "    ds = xr.open_dataset(path, engine='selafin')\n",
    "    tri = ds.attrs['ikle2'] - 1\n",
    "    lons = ds.x.values[tri].T\n",
    "    lats = ds.y.values[tri].T\n",
    "    depths = - ds.B.isel(time=0).values[tri].T\n",
    "    skews, base_cfls = get_skews_and_base_cfls(lons=lons, lats=lats, depths=depths)\n",
    "    return skews, base_cfls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = {\"v1.2\" : [\"/home/tomsail/Documents/work/models/meshes/slf/v1.2.slf\", 50],\n",
    "         \"v2.2\" : [\"/home/tomsail/Documents/work/models/meshes/slf/v2.2.slf\", 30],\n",
    "         \"v1.5\" : [\"/home/tomsail/Documents/work/models/meshes/slf/v1.5.slf\", 50],\n",
    "         \"v2.3\" : [\"/home/tomsail/Documents/work/models/meshes/slf/v2.3.slf\", 30],\n",
    "         \"v3.1\" : [\"/home/tomsail/Documents/work/models/meshes/slf/v3.1.slf\", 20]}\n",
    "\n",
    "plots = []\n",
    "for v in files.keys():\n",
    "    file = files[v][0]\n",
    "    ideal_dt = files[v][1]\n",
    "    ds = xr.open_dataset(file, engine='selafin')\n",
    "    skews, base_cfls = get_skews_and_base_cfls_from_path(file)\n",
    "    CFL_THRESHOLD = 1\n",
    "    print(v)\n",
    "    for dt in (1, 10, 20, 30, 50, 75, 100, 120, 150, 200, 300, 400):\n",
    "        violations = (base_cfls * dt > CFL_THRESHOLD).sum()\n",
    "        print(f\"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%\")\n",
    "    _skews = pd.DataFrame({\"skew\": skews}).hvplot.hist(bins=40,bin_range=[2.5,3.5]).opts(title = v + \" skewness\")\n",
    "    _cfls = pd.DataFrame({\"cfls\": base_cfls * ideal_dt}).hvplot.hist(bins=40, bin_range = [0,1.5]).opts(title = v + \" CFL\")\n",
    "    both = (_skews + _cfls)\n",
    "    plots.append(both)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Layout(plots).cols(2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tri = ds.attrs['ikle2'] - 1\n",
    "nodes = pd.DataFrame(np.vstack((ds.x, ds.y, ds.B.isel(time=0))).T, columns=[\"lon\", \"lat\", \"depth\"])\n",
    "elements = pd.DataFrame(np.vstack( (np.ones(len(tri))* 3, tri.T)).T , columns=[\"no_nodes\", \"n1\", \"n2\", \"n3\"])\n",
    "elements = elements.assign(base_cfl=base_cfls)\n",
    "elements.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "min_cfl_per_node = pd.concat([\n",
    "    elements[[\"n1\", \"base_cfl\"]].groupby([\"n1\"]).base_cfl.min(),\n",
    "    elements[[\"n2\", \"base_cfl\"]].groupby([\"n2\"]).base_cfl.min(),\n",
    "    elements[[\"n3\", \"base_cfl\"]].groupby([\"n3\"]).base_cfl.min(),\n",
    "], axis=1).min(axis=1)\n",
    "min_cfl_per_node.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt = 30\n",
    "df = nodes.assign(\n",
    "    cfl=min_cfl_per_node * dt,\n",
    "    # CFL_violation nodes have a value of 1 if there is no violation and 4 if there is a violation. \n",
    "    # We do this in order to plot the points with a different size\n",
    "    cfl_violation=((min_cfl_per_node * dt > CFL_THRESHOLD) * 3) + 1   \n",
    ")\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot = df[df.cfl_violation == 4].hvplot.points(\n",
    "    'lon', \n",
    "    'lat',\n",
    "    c=\"depth\",\n",
    "    cmap=\"jet\",\n",
    "    geo=True,\n",
    "    tiles=\"EsriImagery\",\n",
    ").options(\n",
    "    width=1200, height=900\n",
    ")\n",
    "len(df[df.cfl_violation == 4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plots = []\n",
    "for v in files.keys():\n",
    "    file = files[v][0]\n",
    "    ideal_dt = files[v][1]\n",
    "    ds = xr.open_dataset(file, engine='selafin')\n",
    "    skews, base_cfls = get_skews_and_base_cfls_from_path(file)\n",
    "    CFL_THRESHOLD = 1\n",
    "    tri = ds.attrs['ikle2'] - 1\n",
    "    # \n",
    "    nodes = pd.DataFrame(np.vstack((ds.x, ds.y, ds.B.isel(time=0))).T, columns=[\"lon\", \"lat\", \"depth\"])\n",
    "    elements = pd.DataFrame(np.vstack( (np.ones(len(tri))* 3, tri.T)).T , columns=[\"no_nodes\", \"n1\", \"n2\", \"n3\"])\n",
    "    elements = elements.assign(cfls=base_cfls*ideal_dt, skews = skews)\n",
    "    min_cfl_per_node = pd.concat([\n",
    "        elements[[\"n1\", \"cfls\"]].groupby([\"n1\"]).cfls.min(),\n",
    "        elements[[\"n2\", \"cfls\"]].groupby([\"n2\"]).cfls.min(),\n",
    "        elements[[\"n3\", \"cfls\"]].groupby([\"n3\"]).cfls.min(),\n",
    "    ], axis=1).min(axis=1)\n",
    "    mean_skew_per_node = pd.concat([\n",
    "        elements[[\"n1\", \"skews\"]].groupby([\"n1\"]).skews.mean(),\n",
    "        elements[[\"n2\", \"skews\"]].groupby([\"n2\"]).skews.mean(),\n",
    "        elements[[\"n3\", \"skews\"]].groupby([\"n3\"]).skews.mean(),\n",
    "    ], axis=1).mean(axis=1)\n",
    "    df = nodes.assign(\n",
    "        cfl=min_cfl_per_node,\n",
    "        # CFL_violation nodes have a value of 1 if there is no violation and 4 if there is a violation. \n",
    "        # We do this in order to plot the points with a different size\n",
    "        cfl_violation=((min_cfl_per_node > CFL_THRESHOLD) * 3) + 1  \n",
    "    )\n",
    "    plot = df[df.cfl_violation == 4].hvplot.points(\n",
    "        'lon', \n",
    "        'lat',\n",
    "        c=\"depth\",\n",
    "        cmap=\"jet\",\n",
    "        geo=True,\n",
    "        tiles=\"EsriImagery\",\n",
    "    ).options(\n",
    "        width=1200, height=900, title = v\n",
    "    )\n",
    "    print(v, len(df[df.cfl_violation == 4]))\n",
    "    plots.append(plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Layout(plots).cols(1)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}