GraysHarbor and NY tidal examples

GraysHarborBC/3Methods.ipynb

@@ -0,0 +1,175 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "from IPython.display import Image"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### The project follows from https://github.com/rjleveque/geoclaw_tides/tree/main/GraysHarbor\n",
+    "In Dr. Randall Leveque's example, the tides are forced inside the domain, the boundary condition is implemented as zero-order extrapolation.\n",
+    "\n",
+    "The maximal tidal forcing is implemented from longitude -124.7 to -124.3 and a linear tapered forcing is implemented from longitude -124.3 to -124.19.\n",
+    "\n",
+    "We call this method ocean forcing."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Simple and Riemann\n",
+    "\n",
+    "We explore two ways of implementing a tidal signal at the boundary, we call the two methods Simple and Riemann\n",
+    "\n",
+    "Geoclaw breaks the domain into smaller rectangle cells, label them $ (i,j)$  in terms of locations, define \n",
+    "\n",
+    "$ s(t): $ tidal signal \n",
+    "\n",
+    "val$ (1, i, j):$  water height on $ (i,j)$ cell\n",
+    "\n",
+    "val$ (2, i, j):$  velocity in x direction on $ (i,j)$ cell\n",
+    "\n",
+    "val$ (3, i, j):$  velocity in y direction on $ (i,j)$ cell\n",
+    "\n",
+    "$\\Delta t:$ next jump in time\n",
+    "\n",
+    "### Simple\n",
+    "(i, j) is the boundary ghost cell\n",
+    "\n",
+    "jump_h$=  (s(t + \\Delta t) - s(t))$.\n",
+    "\n",
+    "#### Code:\n",
+    "\n",
+    "val(1, i, j) = val(1, i + 1, j) + jump_h\n",
+    "\n",
+    "val(2, i, j) = val(2, i + 1, j)\n",
+    "\n",
+    "val(3, i, j) = 0\n",
+    "\n",
+    "### Riemann\n",
+    "\n",
+    "Geoclaw uses shallow water Riemann solver, so this method sets $h_m - h_r$ in the Riemann problem as the tidal signal.\n",
+    "\n",
+    "Ignore the velocity in y direction (since tidal signal is traveling not in the y direction), 1D shallow water equations are \n",
+    "\n",
+    "\n",
+    "$\\begin{aligned} h_{t}+(h u)_{x} &=0 \\\\(h u)_{t}+\\left(h u^{2}+\\frac{1}{2} g h^{2}\\right)_{x} &=0 . \\end{aligned}$\n",
+    "\n",
+    "\n",
+    "\n",
+    "\n",
+    "The Riemann problem for 1D shallow water equations involve giving $q_{l}$ and $q_{r}$, solve for $q_{m}$\n",
+    "\n",
+    "\n",
+    "$q_{l}=\\left[\\begin{array}{c}h_{l} \\\\ u_{l}h_{l} \\end{array}\\right],\\quad q_{m}=\\left[\\begin{array}{c}h_{m} \\\\  u_{m}h_{m}\\end{array}\\right], \\quad q_{r}=\\left[\\begin{array}{c}h_{r} \\\\ u_{r}h_{r} \\end{array}\\right]$\n",
+    "\n",
+    "In our example, $q_{l}$ is the values of the ghost cell that we want to implement, $q_{r}$ is the information given inside the physical domain next to the ghost cell. \n",
+    "\n",
+    "We assume $h_m - h_r = $ jump_h $=  (s(t + \\Delta t) - s(t))$.\n",
+    "\n",
+    "So in our example, the question is, given $q_r$, set $q_l$ such that $h_m = h_r +$ jump_h, we also want the 1-wave to not travel to the right. \n",
+    "\n",
+    "The solution is not unique, the easiest way is to find the unique $u_m$, then set $q_l = q_m$. \n",
+    "\n",
+    "A website is attached below which is an interactive view of the Riemann problem, it lets you set $q_l$, $q_r$ and see $q_m$\n",
+    "\n",
+    "http://www.clawpack.org/riemann_book/phase_plane/shallow_water_small.html\n",
+    "\n",
+    "\n",
+    "##### Solve $u_m$\n",
+    "Case 1: $h_{m}>h_{r}$\n",
+    "\n",
+    "By Lax entropy condition, we have a 2-shock if $h_{m}>h_{r}$\n",
+    "\n",
+    "For a shock wave moving at speed $s$, Rankine--Hugoniot jump conditions gives\n",
+    "\n",
+    "$$\\begin{aligned} s\\left(h_{r}-h_{m}\\right) &=h_{r} u_{r}-h_{m} u_{m} \\\\ s\\left(h_{r} u_{r}-h_{m} u_{m}\\right) &=h_{r} u_{r}^{2}-h_{m} u_{m}^{2}+\\frac{g}{2}\\left(h_{r}^{2}-h^{2}\\right) \\end{aligned}$$\n",
+    "\n",
+    "\n",
+    "Combine the two equations, let $\\alpha = h_m - h_r$\n",
+    "\n",
+    "$$u_{m}=\\frac{h_{r} u_{r}+\\alpha\\left[u_{r}-\\sqrt{g h_{r}\\left(1+\\frac{\\alpha}{h_{r}}\\right)\\left(1+\\frac{\\alpha}{2 h_{r}}\\right)}\\right].}{h_{m}}$$\n",
+    "\n",
+    "\n",
+    "Case 2: $h_{m} < h_{r}$\n",
+    "\n",
+    "It's 2-rarefaction, since the variation within the rarefaction wave is at all points proportional to the corresponding eigenvector $r_p$, the solution can be found by solving $\\tilde{q}^{\\prime}(\\xi)= r^{p}(\\tilde{q}(\\xi))$, where $\\tilde{q}(\\xi)$ is a parameterization of the solution\n",
+    "\n",
+    "The eigenvectors are $r^{1}=\\left[\\begin{array}{c}1 \\\\ u-\\sqrt{g h}\\end{array}\\right], \\quad r^{2}=\\left[\\begin{array}{c}1 \\\\ u+\\sqrt{g h}\\end{array}\\right]$\n",
+    "\n",
+    "Consider $r^{1}$, then $\\tilde{q}^{\\prime}(\\xi)= r^{p}(\\tilde{q}(\\xi))$ is\n",
+    "\n",
+    "\n",
+    "\n",
+    "$$\\begin{aligned} h^{\\prime}(\\xi) &=q_{1}^{\\prime}(\\xi)=1 \\\\(h u)^{\\prime}(\\xi) &=q_{2}^{\\prime}(\\xi)=u \\pm \\sqrt{g h}=\\tilde{q}_{2} / \\tilde{q}_{1}-\\sqrt{g \\tilde{q}_{1}} \\end{aligned}$$\n",
+    "\n",
+    "\n",
+    "Fixing $(u_r,u_r h_r)$,\n",
+    "\n",
+    "$$ u_{m}= \\frac{h_{m} u_{r}-2 h_{r}\\left(\\sqrt{g h_{r}}-\\sqrt{g h_{m}}\\right)}{h_{m}}$$\n",
+    "\n",
+    "\n",
+    "So the logic of setting ghost cells goes like this,\n",
+    "\n",
+    "#### Code:\n",
+    "\n",
+    "jump_h$=  (s(t + \\Delta t) - s(t))$.\n",
+    "\n",
+    "val(3, i, j) = 0\n",
+    "\n",
+    "val(1, i, j) = val(1, nxl + 1, j) + jump_h\n",
+    "\n",
+    "h_r = val(1, nxl + 1, j)\n",
+    "\n",
+    "h_m = val(1, i, j)\n",
+    "\n",
+    "u_r = val(2, nxl + 1, j)\n",
+    "\n",
+    "if (h_r < h_m)\n",
+    "\n",
+    "    val(2, i, j) = (h_r*u_r + jump_h*(u_r - sqrt(9.81*h_r*(1+jump_h/h_r)*(1+jump_h/(2*h_r)))))/h_m\n",
+    "                      \n",
+    "else\n",
+    "\n",
+    "    val(2, i, j) = (h_m*u_r - 2*h_r*(sqrt(9.81*h_r)-sqrt(9.81*h_m)))/h_m"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

GraysHarborBC/README.md

@@ -0,0 +1,31 @@
+## This directory simulate tides in Grays Harbor, WA with 3 methods. 
+
+### First, a sine tidal signal is implemented to get 
+
+1.The high tide time lag between the sine wave and NOAA station 1102 at Westport   
+
+2.The high tide amplitude at NOAA station 1102
+
+This is done in ./sine
+
+### Second, simulates the King Tide event from December, 2015
+
+The tidal condition is implemented using the existing NOAA 1102 tides prediction data, but the data is shifted by the time lag and divided by the amplitude observed in the sine example.
+
+This is done in ./kingtide2015
+
+## 3Methods.ipynb
+
+It explains "Simple", "Riemann" and ocean forcing method.
+
+## `topo` directory
+
+It creates Grays Harbor topography. 
+
+
+## `comparison` directory
+
+It compares tidal predictions of the 3 methods (Simple, Riemann and ocean forcing) used for King Tide event from December, 2015
+
+
+Feel free to contact yz3889@columbia.edu if you have any questions about the code.

GraysHarborBC/comparison/README.md

@@ -0,0 +1,5 @@
+## Compare 3 methods, Randall (ocean forcing), Simple and Riemann.
+
+Move 6 files simple1102, simple1187, riemann1102, riemann1187, randall1102 and randall1187 (errors generated by KingTide.ipynb) to this directory
+
+Run comparison.py to generate error graphs and print 1 & 2 norm errors

GraysHarborBC/comparison/compare.py

@@ -0,0 +1,53 @@
+import numpy as np
+import matplotlib.pyplot as plt
+time = np.genfromtxt ('hours', delimiter=",")
+
+randall1102 = np.genfromtxt ('randall1102', delimiter=",")
+randall1187 = np.genfromtxt ('randall1187', delimiter=",")
+
+simple1102 =np.genfromtxt ('simple1102', delimiter=",")
+simple1187 = np.genfromtxt ('simple1187', delimiter=",")
+
+riemann1102 =np.genfromtxt ('riemann1102', delimiter=",")
+riemann1187 =np.genfromtxt ('riemann1187', delimiter=",")
+
+plt.figure(figsize=(13,5))
+plt.plot(time,randall1102, label='Randall1102')
+plt.plot(time,simple1102, label='simple1102')
+plt.plot(time,riemann1102, label='Riemann1102')
+
+
+plt.xlim(0,72)
+plt.ylim(-0.5,0.5)
+plt.legend(loc='upper left', fontsize=8)
+plt.xlabel('hours')
+plt.ylabel('m')
+plt.title('NOAA minus GeoClaw at 1102, Dec. 21-23, 2015')
+plt.grid(True)
+plt.savefig('1102error.png')
+
+plt.figure(figsize=(13,5))
+
+plt.plot(time,randall1187, label='Randall1187')
+plt.plot(time,simple1187, label='simple1187')
+plt.plot(time,riemann1187, label='Riemann1187')
+
+
+
+plt.xlim(0,72)
+plt.ylim(-0.5,0.5)
+plt.legend(loc='upper left', fontsize=8)
+plt.xlabel('hours')
+plt.ylabel('m')
+plt.title('NOAA minus GeoClaw at 1187, Dec. 21-23, 2015')
+plt.grid(True)
+
+plt.savefig('1187error.png')
+
+print("Randall, 1 norm error at 1102:", np.linalg.norm(randall1102,1))
+print("Simple, 1 norm error at 1102:",np.linalg.norm(simple1102,1))
+print("Riemann, 1 norm error at 1102:",np.linalg.norm(riemann1102,1))
+
+print("Randall, 2 norm error at 1102:", np.linalg.norm(randall1102,2))
+print("Simple, 2 norm error at 1102:",np.linalg.norm(simple1102,2))
+print("Riemann, 2 norm error at 1102:",np.linalg.norm(riemann1102,2))