From 0041e720c17f022f7384a80fe14aec5f76c6b670 Mon Sep 17 00:00:00 2001 From: Panos Mavrogiorgos Date: Mon, 6 May 2024 12:44:07 +0300 Subject: [PATCH] Add CFL-schism --- CFL-schism.html | 9380 ++++++++++++++++++++++++++++++++++++++++++++++ CFL-schism.ipynb | 251 ++ index.html | 1 + 3 files changed, 9632 insertions(+) create mode 100644 CFL-schism.html create mode 100644 CFL-schism.ipynb diff --git a/CFL-schism.html b/CFL-schism.html new file mode 100644 index 0000000..9e8a12a --- /dev/null +++ b/CFL-schism.html @@ -0,0 +1,9380 @@ + + + + + +CFL-schism + + + + + + + + + + + + +
+
+ + diff --git a/CFL-schism.ipynb b/CFL-schism.ipynb new file mode 100644 index 0000000..3a40291 --- /dev/null +++ b/CFL-schism.ipynb @@ -0,0 +1,251 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "00e18507-11e9-4ef3-a165-d9beb2cdd290", + "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 numpy as np\n", + "import pandas as pd\n", + "import pymap3d\n", + "import pyposeidon.mesh as pmesh\n", + "\n", + "hv.extension(\"bokeh\")\n", + "np.set_printoptions(suppress=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca4352df-6db1-4087-bd8a-64b9346e83fd", + "metadata": {}, + "outputs": [], + "source": [ + "def parse_hgrid_nodes(path: os.PathLike[str] | str) -> pd.DataFrame:\n", + " with open(path, \"rb\") as fd:\n", + " _ = fd.readline()\n", + " _, no_points = map(int, fd.readline().strip().split(b\" \"))\n", + " content = io.BytesIO(b''.join(next(fd) for _ in range(no_points)))\n", + " nodes = pd.read_csv(\n", + " content,\n", + " engine=\"pyarrow\",\n", + " sep=\"\\t\",\n", + " header=None,\n", + " names=[\"lon\", \"lat\", \"depth\"],\n", + " index_col=0\n", + " )\n", + " nodes = nodes.reset_index(drop=True)\n", + " return nodes\n", + " \n", + "def parse_hgrid_elements3(path: os.PathLike[str] | str) -> pd.DataFrame:\n", + " with open(path, \"rb\") as fd:\n", + " _ = fd.readline()\n", + " no_elements, no_points = map(int, fd.readline().strip().split(b\" \"))\n", + " for _ in range(no_points):\n", + " next(fd) \n", + " content = io.BytesIO(b''.join(next(fd) for _ in range(no_elements)))\n", + " elements = pd.read_csv(\n", + " content,\n", + " engine=\"pyarrow\",\n", + " sep=\"\\t\",\n", + " header=None,\n", + " names=[\"no_nodes\", \"n1\", \"n2\", \"n3\"],\n", + " index_col=0\n", + " )\n", + " elements = elements.assign(\n", + " n1=elements.n1 - 1,\n", + " n2=elements.n2 - 1,\n", + " n3=elements.n3 - 1,\n", + " ).reset_index(drop=True)\n", + " return elements" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d628ca1-7925-44b3-84de-90484e55cb0f", + "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, )\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", + " nodes = parse_hgrid_nodes(path)\n", + " elements = parse_hgrid_elements3(path)\n", + " tri = elements[[\"n1\", \"n2\", \"n3\"]].values\n", + " lons = nodes.lon.values[tri].T\n", + " lats = nodes.lat.values[tri].T\n", + " depths = nodes.depth.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, + "id": "322443c2-14aa-435e-b0c6-5d6a42dfd1db", + "metadata": {}, + "outputs": [], + "source": [ + "path = \"/home/panos/Prog/poseidon/seareport_meshes/meshes/global-v0.1.gr3\"\n", + "path = \"/home/panos/Prog/poseidon/seareport_meshes/meshes/global-v0.gr3\"\n", + "path = \"/home/panos/Prog/git/schism/src/Utility/Grid_Scripts/hgrid.gr3\"\n", + "path = \"/home/panos/Prog/poseidon/seareport_meshes/meshes/global-v0.2.gr3\"\n", + "skews, base_cfls = get_skews_and_base_cfls_from_path(path)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d25aa8d7-61b3-452e-963d-7d918ed0c273", + "metadata": {}, + "outputs": [], + "source": [ + "CFL_THRESHOLD = 0.4\n", + "for dt in (1, 50, 75, 100, 120, 150, 200, 300, 400, 600, 900, 1200, 1800, 3600):\n", + " violations = (base_cfls * dt < CFL_THRESHOLD).sum()\n", + " print(f\"{dt:>4d} {violations:>12d} {violations / len(base_cfls) * 100:>8.2f}%\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8016a0c5-bb91-44a0-abcf-579e19c5bc22", + "metadata": {}, + "outputs": [], + "source": [ + "pd.DataFrame({\"skew\": skews}).describe()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d03b21c-af6c-4e31-8040-63ec52affab2", + "metadata": {}, + "outputs": [], + "source": [ + "df = pd.DataFrame({\"cfl\": base_cfls * 400})\n", + "df[df.cfl < 0.4].describe()\n", + "df[df.cfl < 0.4].hvplot.hist()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8988ce4e-7607-4f65-8f93-2c7f0b9cac60", + "metadata": {}, + "outputs": [], + "source": [ + "nodes = parse_hgrid_nodes(path)\n", + "elements = parse_hgrid_elements3(path)\n", + "elements = elements.assign(base_cfl=base_cfls)\n", + "elements.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad6461d1-4a49-4b78-9703-e2dc989c1f2b", + "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, + "id": "8def9317-af88-49a6-b867-0c935ca6ecc7", + "metadata": {}, + "outputs": [], + "source": [ + "dt = 600\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, + "id": "ee158896-1dd1-463f-a623-4b8b4ca9425d", + "metadata": {}, + "outputs": [], + "source": [ + "df.hvplot.points(\n", + " 'lon', \n", + " 'lat',\n", + " c=\"cfl_violation\",\n", + " cmap=\"colorblind\",\n", + " geo=True,\n", + " tiles=True,\n", + ").options(\n", + " width=900, height=600\n", + ").opts(\n", + " hv.opts.Points(size=hv.dim('cfl_violation'))\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bb6f5592-81be-4acb-994b-0d2142ecd34f", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "pos_dev", + "language": "python", + "name": "pos_dev" + }, + "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": 5 +} diff --git a/index.html b/index.html index 03f86ce..02001fb 100644 --- a/index.html +++ b/index.html @@ -45,6 +45,7 @@