diff --git a/.gitignore b/.gitignore deleted file mode 100644 index 259148f..0000000 --- a/.gitignore +++ /dev/null @@ -1,32 +0,0 @@ -# Prerequisites -*.d - -# Compiled Object files -*.slo -*.lo -*.o -*.obj - -# Precompiled Headers -*.gch -*.pch - -# Compiled Dynamic libraries -*.so -*.dylib -*.dll - -# Fortran module files -*.mod -*.smod - -# Compiled Static libraries -*.lai -*.la -*.a -*.lib - -# Executables -*.exe -*.out -*.app diff --git a/GraysHarborBC/3Methods.ipynb b/GraysHarborBC/3Methods.ipynb new file mode 100644 index 0000000..12e5fc5 --- /dev/null +++ b/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 +} diff --git a/GraysHarborBC/README.md b/GraysHarborBC/README.md new file mode 100644 index 0000000..9a0d9ba --- /dev/null +++ b/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. diff --git a/GraysHarborBC/comparison/README.md b/GraysHarborBC/comparison/README.md new file mode 100644 index 0000000..2f27e7b --- /dev/null +++ b/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 diff --git a/GraysHarborBC/comparison/compare.py b/GraysHarborBC/comparison/compare.py new file mode 100644 index 0000000..5229673 --- /dev/null +++ b/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)) diff --git a/GraysHarborBC/comparison/hours b/GraysHarborBC/comparison/hours new file mode 100644 index 0000000..f16d255 --- /dev/null +++ b/GraysHarborBC/comparison/hours @@ -0,0 +1,720 @@ +0.0 +0.1 +0.2 +0.3 +0.4 +0.5 +0.6 +0.7 +0.8 +0.9 +1.0 +1.1 +1.2 +1.3 +1.4 +1.5 +1.6 +1.7 +1.8 +1.9 +2.0 +2.1 +2.2 +2.3 +2.4 +2.5 +2.6 +2.7 +2.8 +2.9 +3.0 +3.1 +3.2 +3.3 +3.4 +3.5 +3.6 +3.7 +3.8 +3.9 +4.0 +4.1 +4.2 +4.3 +4.4 +4.5 +4.6 +4.7 +4.8 +4.9 +5.0 +5.1 +5.2 +5.3 +5.4 +5.5 +5.6 +5.7 +5.8 +5.9 +6.0 +6.1 +6.2 +6.3 +6.4 +6.5 +6.6 +6.7 +6.8 +6.9 +7.0 +7.1 +7.2 +7.3 +7.4 +7.5 +7.6 +7.7 +7.8 +7.9 +8.0 +8.1 +8.2 +8.3 +8.4 +8.5 +8.6 +8.7 +8.8 +8.9 +9.0 +9.1 +9.2 +9.3 +9.4 +9.5 +9.6 +9.7 +9.8 +9.9 +10.0 +10.1 +10.2 +10.3 +10.4 +10.5 +10.6 +10.7 +10.8 +10.9 +11.0 +11.1 +11.2 +11.3 +11.4 +11.5 +11.6 +11.7 +11.8 +11.9 +12.0 +12.1 +12.2 +12.3 +12.4 +12.5 +12.6 +12.7 +12.8 +12.9 +13.0 +13.1 +13.2 +13.3 +13.4 +13.5 +13.6 +13.7 +13.8 +13.9 +14.0 +14.1 +14.2 +14.3 +14.4 +14.5 +14.6 +14.7 +14.8 +14.9 +15.0 +15.1 +15.2 +15.3 +15.4 +15.5 +15.6 +15.7 +15.8 +15.9 +16.0 +16.1 +16.2 +16.3 +16.4 +16.5 +16.6 +16.7 +16.8 +16.9 +17.0 +17.1 +17.2 +17.3 +17.4 +17.5 +17.6 +17.7 +17.8 +17.9 +18.0 +18.1 +18.2 +18.3 +18.4 +18.5 +18.6 +18.7 +18.8 +18.9 +19.0 +19.1 +19.2 +19.3 +19.4 +19.5 +19.6 +19.7 +19.8 +19.9 +20.0 +20.1 +20.2 +20.3 +20.4 +20.5 +20.6 +20.7 +20.8 +20.9 +21.0 +21.1 +21.2 +21.3 +21.4 +21.5 +21.6 +21.7 +21.8 +21.9 +22.0 +22.1 +22.2 +22.3 +22.4 +22.5 +22.6 +22.7 +22.8 +22.9 +23.0 +23.1 +23.2 +23.3 +23.4 +23.5 +23.6 +23.7 +23.8 +23.9 +24.0 +24.1 +24.2 +24.3 +24.4 +24.5 +24.6 +24.7 +24.8 +24.9 +25.0 +25.1 +25.2 +25.3 +25.4 +25.5 +25.6 +25.7 +25.8 +25.9 +26.0 +26.1 +26.2 +26.3 +26.4 +26.5 +26.6 +26.7 +26.8 +26.9 +27.0 +27.1 +27.2 +27.3 +27.4 +27.5 +27.6 +27.7 +27.8 +27.9 +28.0 +28.1 +28.2 +28.3 +28.4 +28.5 +28.6 +28.7 +28.8 +28.9 +29.0 +29.1 +29.2 +29.3 +29.4 +29.5 +29.6 +29.7 +29.8 +29.9 +30.0 +30.1 +30.2 +30.3 +30.4 +30.5 +30.6 +30.7 +30.8 +30.9 +31.0 +31.1 +31.2 +31.3 +31.4 +31.5 +31.6 +31.7 +31.8 +31.9 +32.0 +32.1 +32.2 +32.3 +32.4 +32.5 +32.6 +32.7 +32.8 +32.9 +33.0 +33.1 +33.2 +33.3 +33.4 +33.5 +33.6 +33.7 +33.8 +33.9 +34.0 +34.1 +34.2 +34.3 +34.4 +34.5 +34.6 +34.7 +34.8 +34.9 +35.0 +35.1 +35.2 +35.3 +35.4 +35.5 +35.6 +35.7 +35.8 +35.9 +36.0 +36.1 +36.2 +36.3 +36.4 +36.5 +36.6 +36.7 +36.8 +36.9 +37.0 +37.1 +37.2 +37.3 +37.4 +37.5 +37.6 +37.7 +37.8 +37.9 +38.0 +38.1 +38.2 +38.3 +38.4 +38.5 +38.6 +38.7 +38.8 +38.9 +39.0 +39.1 +39.2 +39.3 +39.4 +39.5 +39.6 +39.7 +39.8 +39.9 +40.0 +40.1 +40.2 +40.3 +40.4 +40.5 +40.6 +40.7 +40.8 +40.9 +41.0 +41.1 +41.2 +41.3 +41.4 +41.5 +41.6 +41.7 +41.8 +41.9 +42.0 +42.1 +42.2 +42.3 +42.4 +42.5 +42.6 +42.7 +42.8 +42.9 +43.0 +43.1 +43.2 +43.3 +43.4 +43.5 +43.6 +43.7 +43.8 +43.9 +44.0 +44.1 +44.2 +44.3 +44.4 +44.5 +44.6 +44.7 +44.8 +44.9 +45.0 +45.1 +45.2 +45.3 +45.4 +45.5 +45.6 +45.7 +45.8 +45.9 +46.0 +46.1 +46.2 +46.3 +46.4 +46.5 +46.6 +46.7 +46.8 +46.9 +47.0 +47.1 +47.2 +47.3 +47.4 +47.5 +47.6 +47.7 +47.8 +47.9 +48.0 +48.1 +48.2 +48.3 +48.4 +48.5 +48.6 +48.7 +48.8 +48.9 +49.0 +49.1 +49.2 +49.3 +49.4 +49.5 +49.6 +49.7 +49.8 +49.9 +50.0 +50.1 +50.2 +50.3 +50.4 +50.5 +50.6 +50.7 +50.8 +50.9 +51.0 +51.1 +51.2 +51.3 +51.4 +51.5 +51.6 +51.7 +51.8 +51.9 +52.0 +52.1 +52.2 +52.3 +52.4 +52.5 +52.6 +52.7 +52.8 +52.9 +53.0 +53.1 +53.2 +53.3 +53.4 +53.5 +53.6 +53.7 +53.8 +53.9 +54.0 +54.1 +54.2 +54.3 +54.4 +54.5 +54.6 +54.7 +54.8 +54.9 +55.0 +55.1 +55.2 +55.3 +55.4 +55.5 +55.6 +55.7 +55.8 +55.9 +56.0 +56.1 +56.2 +56.3 +56.4 +56.5 +56.6 +56.7 +56.8 +56.9 +57.0 +57.1 +57.2 +57.3 +57.4 +57.5 +57.6 +57.7 +57.8 +57.9 +58.0 +58.1 +58.2 +58.3 +58.4 +58.5 +58.6 +58.7 +58.8 +58.9 +59.0 +59.1 +59.2 +59.3 +59.4 +59.5 +59.6 +59.7 +59.8 +59.9 +60.0 +60.1 +60.2 +60.3 +60.4 +60.5 +60.6 +60.7 +60.8 +60.9 +61.0 +61.1 +61.2 +61.3 +61.4 +61.5 +61.6 +61.7 +61.8 +61.9 +62.0 +62.1 +62.2 +62.3 +62.4 +62.5 +62.6 +62.7 +62.8 +62.9 +63.0 +63.1 +63.2 +63.3 +63.4 +63.5 +63.6 +63.7 +63.8 +63.9 +64.0 +64.1 +64.2 +64.3 +64.4 +64.5 +64.6 +64.7 +64.8 +64.9 +65.0 +65.1 +65.2 +65.3 +65.4 +65.5 +65.6 +65.7 +65.8 +65.9 +66.0 +66.1 +66.2 +66.3 +66.4 +66.5 +66.6 +66.7 +66.8 +66.9 +67.0 +67.1 +67.2 +67.3 +67.4 +67.5 +67.6 +67.7 +67.8 +67.9 +68.0 +68.1 +68.2 +68.3 +68.4 +68.5 +68.6 +68.7 +68.8 +68.9 +69.0 +69.1 +69.2 +69.3 +69.4 +69.5 +69.6 +69.7 +69.8 +69.9 +70.0 +70.1 +70.2 +70.3 +70.4 +70.5 +70.6 +70.7 +70.8 +70.9 +71.0 +71.1 +71.2 +71.3 +71.4 +71.5 +71.6 +71.7 +71.8 +71.9 diff --git a/GraysHarborBC/kingtide2015/KingTide.ipynb b/GraysHarborBC/kingtide2015/KingTide.ipynb new file mode 100644 index 0000000..ea5b66c --- /dev/null +++ b/GraysHarborBC/kingtide2015/KingTide.ipynb @@ -0,0 +1,352 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note: this notebook is run for the Simple case. The titles of the graphs and exported error file names need to change for Riemann and ocean forcing case.\n", + "\n", + "tshift and signul need to match what you get in ../sine" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Timeshift \n", + "tshift = 3220\n", + "#Signal should be multiplied by \n", + "sigmul = 1.04" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%pylab inline\n", + "import datetime\n", + "from scipy.interpolate import interp1d\n", + "from clawpack.geoclaw import util\n", + "import clawpack.pyclaw.gauges as gauges" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Download tide gauge data for Westport and Aberdeen" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "station_westport = '9441102'\n", + "MLW_westport = -1.068 # m relative to MTL\n", + "MHW_westport = 1.068 # m relative to MTL\n", + "\n", + "station_aberdeen = '9441187'\n", + "MLW_aberdeen = -1.21 # m relative to MTL\n", + "MHW_aberdeen = 1.21 # m relative to MTL\n", + "\n", + "begin_date = datetime.datetime(2015,12,21)\n", + "end_date = datetime.datetime(2015,12,26)\n", + "time_zone = 'GMT'\n", + "datum = 'MTL'\n", + "units = 'metric'\n", + "cache_dir = '.'\n", + "\n", + "time_westport, eta_westport, etap_westport = \\\n", + " util.fetch_noaa_tide_data(station_westport, begin_date, end_date,\n", + " time_zone, datum, units, cache_dir)\n", + "\n", + "time_aberdeen, eta_aberdeen, etap_aberdeen = \\\n", + " util.fetch_noaa_tide_data(station_aberdeen, begin_date, end_date,\n", + " time_zone, datum, units, cache_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "figure(figsize=(13,6))\n", + "plot(time_westport, etap_westport, 'b', label='Westport %s' % station_westport)\n", + "plot(time_westport, MHW_westport*ones(time_westport.shape),'b--', label=\"MLW/MHW Westport\")\n", + "plot(time_westport, MLW_westport*ones(time_westport.shape),'b--')\n", + "\n", + "plot(time_aberdeen, etap_aberdeen, 'r', label='Aberdeen %s' % station_aberdeen)\n", + "plot(time_aberdeen, MHW_aberdeen*ones(time_aberdeen.shape),'r--', label=\"MLW/MHW Aberdeen\")\n", + "plot(time_aberdeen, MLW_aberdeen*ones(time_aberdeen.shape),'r--')\n", + "\n", + "grid(True)\n", + "xticks(rotation=20)\n", + "legend()\n", + "title('NOAA predicted tides in Grays Harbor')\n", + "ylabel('meters relative to %s' % datum);" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Convert datetimes to elapsed seconds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# convert to elapsed seconds:\n", + "dt_westport = time_westport - time_westport[0]\n", + "t_westport = array([dt.item().total_seconds() for dt in dt_westport])\n", + "\n", + "dt_aberdeen = time_aberdeen - time_aberdeen[0]\n", + "t_aberdeen = array([dt.item().total_seconds() for dt in dt_aberdeen])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "figure(figsize=(13,6))\n", + "plot(t_westport, etap_westport, 'b', label='Westport %s' % station_westport)\n", + "plot(t_aberdeen, etap_aberdeen, 'r', label='Aberdeen %s' % station_aberdeen)\n", + "\n", + "\n", + "grid(True)\n", + "xticks(rotation=20)\n", + "xlabel('seconds')\n", + "legend()\n", + "title('NOAA predicted tides in Grays Harbor')\n", + "ylabel('meters relative to %s' % datum);" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#export data used for graphing/ocean forcing\n", + "t = t_westport - tshift\n", + "eta = etap_westport * sigmul\n", + "\n", + "fname = 'tidedata_dec2015.txt'\n", + "\n", + "etaprime = (eta[2:]-eta[:-2])/(t[2:]-t[:-2]) # central diff at interior t\n", + "etaprime = hstack(((eta[1]-eta[0])/(t[1]-t[0]), etaprime, \\\n", + " (eta[-1]-eta[-2])/(t[-1]-t[-2])))\n", + "\n", + "d = vstack((t, eta, etaprime)).T\n", + "\n", + "mt = d.shape[0]\n", + "header = '%i # mt' % mt\n", + "savetxt(fname, d, header=header, comments='')\n", + "print('Created ',fname)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Linear interpolate the tidal signal\n", + "For Riemann and Simple method, signal used on the left boundary is exported as tidal_signal.txt below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t_index_shift = int(tshift/360)\n", + "signal = eta[t_index_shift:]\n", + "numpy.savetxt(\"tidal_signal.txt\", signal, fmt=\"%s\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Since the data is discrete, we plot the linear interpolation of the data. \n", + "\n", + "In bc2amr.f90, the height jump at time t with jump $\\Delta t$ is implemented to be the derivative of the linear interpolation at t multiplied by $\\Delta t$. It's fairly accurate since $\\Delta t << 360s$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#linear interpolation of the tidal signals\n", + "from scipy.interpolate import interp1d\n", + "f = interp1d(t_westport -tshift, etap_westport)\n", + "x = np.linspace(0, 400000, num=100000, endpoint=True)\n", + "plt.plot(x, sigmul*f(x))\n", + "xlabel('seconds')\n", + "ylabel('meters relative to %s' % datum);\n", + "title('Continuous Tidal Signal')\n", + "plt.show()\n", + "print(\"Ocean level at time 0 =\",f(0))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Note : User needs to update ocean level at time 0 in setrun.py, it's -0.457 currently.\n", + "## Plot GeoClaw results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "import numpy\n", + "\n", + "outdir = '_output'\n", + "\n", + "figure(figsize=(13,5))\n", + "colors = ['b','r','m','g']\n", + "\n", + "for k,gaugeno in enumerate([1102,1187]):\n", + " gauge = gauges.GaugeSolution(gaugeno, outdir)\n", + " t = gauge.t / 3600. # convert to hours\n", + " q = gauge.q\n", + " eta = q[3,:]\n", + " if gaugeno==1102:\n", + " plot(t, eta, colors[k], label='GeoClaw 1102 - Westport')\n", + "\n", + " elif gaugeno==1187:\n", + " plot(t, eta, colors[k], label='GeoClaw 1187 - Aberdeen')\n", + " \n", + " else:\n", + " plot(t, eta, colors[k], label='GeoClaw Gauge %s' % gaugeno)\n", + " \n", + " \n", + "plot(t_westport/3600., etap_westport, 'c', label='NOAA 1102 - Westport')\n", + "plot(t_aberdeen/3600., etap_aberdeen, 'm', label='NOAA 1187 - Aberdeen')\n", + "\n", + "\n", + "xlim(0,72)\n", + "ylim(-2.5,2.5)\n", + "legend(loc='upper left', fontsize=8)\n", + "xlabel('hours')\n", + "ylabel('Surface relative to MTL (m)')\n", + "title('Grays Harbor tides (simple), Dec. 21-23, 2015')\n", + "grid(True)\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#Show error, export error\n", + "\n", + "def find_nearest(array, value):\n", + " array = numpy.asarray(array)\n", + " idx = (numpy.abs(array - value)).argmin()\n", + " return idx\n", + "\n", + "from operator import sub\n", + "outdir = '_output'\n", + "\n", + "figure(figsize=(13,5))\n", + "colors = ['b','r','m','g']\n", + "\n", + "for k,gaugeno in enumerate([1102,1187]):\n", + " gauge = gauges.GaugeSolution(gaugeno, outdir)\n", + " t = gauge.t / 3600. # convert to hours\n", + " q = gauge.q\n", + " eta = q[3,:]\n", + " \n", + " if gaugeno==1102:\n", + " i = 0\n", + " x_select = []\n", + " t_select = []\n", + "\n", + " while i < 720: \n", + " x_select.append(eta[find_nearest(t, i/10)]) \n", + " t_select.append(t[find_nearest(t, i/10)])\n", + " \n", + " i += 1 \n", + " plot(t_westport[0:720]/3600, list( map(sub, etap_westport[0:720], x_select) ), \n", + " colors[k], label='1102, NOAA minus GeoClaw') \n", + "\n", + " numpy.savetxt(\"simple1102\", \n", + " list( map(sub, etap_westport[0:720], x_select) ),\n", + " delimiter =\", \", \n", + " fmt ='% s') \n", + " \n", + " \n", + " \n", + " elif gaugeno==1187:\n", + " i = 0\n", + " x_select = []\n", + " t_select = []\n", + " while i < 720: \n", + " x_select.append(eta[find_nearest(t, i/10)]) \n", + " t_select.append(t[find_nearest(t, i/10)])\n", + " \n", + " i += 1 \n", + " plot(t_westport[0:720]/3600, list( map(sub, etap_aberdeen[0:720], x_select) ), \n", + " colors[k], label='1187, NOAA minus GeoClaw') \n", + " \n", + " numpy.savetxt(\"simple1187\", \n", + " list( map(sub, etap_aberdeen[0:720], x_select) ),\n", + " delimiter =\", \", \n", + " fmt ='% s') \n", + " \n", + " \n", + "xlim(0,72)\n", + "ylim(-2.5,2.5)\n", + "legend(loc='upper left', fontsize=8)\n", + "xlabel('hours')\n", + "ylabel('m')\n", + "title('Grays Harbor Difference, Dec. 21-23, 2015')\n", + "grid(True)\n" + ] + } + ], + "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 +} diff --git a/GraysHarborBC/kingtide2015/Makefile b/GraysHarborBC/kingtide2015/Makefile new file mode 100644 index 0000000..7e6fb5b --- /dev/null +++ b/GraysHarborBC/kingtide2015/Makefile @@ -0,0 +1,69 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + +MODULES = \ + setprob_module.f90 \ + +SOURCES = \ + bc2amr.f90 \ + setprob.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + python maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls diff --git a/GraysHarborBC/kingtide2015/README.md b/GraysHarborBC/kingtide2015/README.md new file mode 100644 index 0000000..9fdd18e --- /dev/null +++ b/GraysHarborBC/kingtide2015/README.md @@ -0,0 +1,14 @@ +The current directory is set up for the Simple method. + +To use Riemann method, replace ./bc2amr.f90 by ./src/Riemann/bc2amr.f90. + +To use ocean forcing method, go to ./src/ocean-forcing, copy every file to the current directory. + +For all methods, + +1. follow instructions in jupyter notebook, run until "Plot GeoClaw results" + +2. run the GeoClaw code + +3. use notebook to examine the results. + diff --git a/GraysHarborBC/kingtide2015/b4run.py b/GraysHarborBC/kingtide2015/b4run.py new file mode 100644 index 0000000..ecdee17 --- /dev/null +++ b/GraysHarborBC/kingtide2015/b4run.py @@ -0,0 +1,60 @@ +""" +Sample b4run.py file. If a file like this is in the rundir directory +used by runclaw.py (e.g. the application directory from which you +execute 'make .output') this will be executed after creating the +outdir and before running the Clawpack code. + +If no b4run.py file is found in the rundir then the environment variable +B4RUN, if set, is used to locate a file used more globally, e.g. you could +set this to point to $HOME/b4ruon.py for a default version to use if there +is no local version, or you can use this version by setting B4RUN to +$CLAW/clawutil/src/python/clawutil/b4run.py + +Use this, for example, to copy any files you want to the outdir, +either because they are needed for the run or to archive them along +with the output. + +This sample version copies the Makefile and any Python or Fortran codes. + +It also creates (or appends to) a runlog.txt file containing information +about this run: the rundir and the time/date the run started. + +Adapt this file to your own needs. +""" + + +def b4run(rundir, outdir): + + import os,sys,glob,shutil,time + + # --------------------------------------------------------------------- + # files to copy to outdir (in addition to *.data files always copied): + to_copy = ['Makefile', '*.py', '*.f*', 'tidedata_dec2015.txt','tidal_signal.txt'] + #to_copy = ['Makefile', '*.py', '*.f*'] + + if rundir != outdir: + for pattern in to_copy: + files = glob.glob(pattern) # files matching pattern + for file in files: + print('Copying %s to %s' % (file, outdir)) + shutil.copy(file, os.path.join(outdir,file)) + + # --------------------------------------------------------------------- + # also add to outdir/runlog.txt with info about when run was done + # and from what directory: + + runlog = os.path.join(outdir, 'runlog.txt') + + tm = time.localtime() + year = str(tm.tm_year).zfill(4) + month = str(tm.tm_mon).zfill(2) + day = str(tm.tm_mday).zfill(2) + hour = str(tm.tm_hour).zfill(2) + minute = str(tm.tm_min).zfill(2) + tz = tm.tm_zone + dt = '%s-%s-%s-%s:%s%s' % (year,month,day,hour,minute,tz) + + with open(runlog,'a') as f: + f.write('Run started at %s\n' % dt) + f.write('from RUNDIR =\n %s\n' % os.path.abspath(rundir)) + diff --git a/GraysHarborBC/kingtide2015/bc2amr.f90 b/GraysHarborBC/kingtide2015/bc2amr.f90 new file mode 100644 index 0000000..78c36ee --- /dev/null +++ b/GraysHarborBC/kingtide2015/bc2amr.f90 @@ -0,0 +1,324 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + use setprob_module, only: a + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt, t_num + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + ! Getting the jump in tidal signal using linear Intepolation on height.txt + t_num = floor(time/360) + delt = possk(level) + jump_h = (a(t_num+2) - a(t_num+1))*delt/360 + + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = 0 + val(1, i, j) = val(1, nxl + 1, j) + jump_h + val(2, i, j) = val(2, nxl + 1, j) + + end do + end do + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/GraysHarborBC/kingtide2015/setplot.py b/GraysHarborBC/kingtide2015/setplot.py new file mode 100644 index 0000000..b90fcb8 --- /dev/null +++ b/GraysHarborBC/kingtide2015/setplot.py @@ -0,0 +1,533 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d + +from clawpack.geoclaw import topotools +from clawpack.visclaw import gaugetools, plottools, colormaps + +cmin_water = -1. +cmax_water = 1. +cmin_speed = 0. +cmax_speed = 2. + +def speed(current_data): + """ + Return a masked array containing the speed only in wet cells. + Surface is eta = h+topo, assumed to be output as 4th column of fort.q + files. + """ + from clawpack.visclaw.geoplot import drytol_default + drytol = getattr(current_data.user, 'drytol', drytol_default) + q = current_data.q + h = q[0,:,:] + hu = q[1,:,:] + hv = q[2,:,:] + eta = q[3,:,:] + + hwet = np.ma.masked_where(h <= drytol, h) + speed = np.sqrt(hu**2 + hv**2) / hwet # will be masked where dry + + try: + # Use mask covering coarse regions if it's set: + m = current_data.mask_coarse + speed = numpy.ma.masked_where(m, speed) + except: + pass + + return speed + +sea_level = -0.468 +tide_data = np.loadtxt('tidedata_dec2015.txt', skiprows=1) +t_tide = tide_data[:,0] / 3600. +#offset = 2.561 # from MLLW to MHW +offset = 0. +eta_tide = interp1d(t_tide, tide_data[:,1]-offset, bounds_error=False, + fill_value = np.nan) + + + +def plot_tide(current_data): + from pylab import linspace,plot,grid,title,cos,pi,xlabel,ylabel,xticks + from pylab import xlim,ylim + tframe = current_data.t / 3600. + tt = linspace(t_tide[0],t_tide[-1],1000) + # forcing: + #tperiod = 16 + #eta_tide = lambda t: 1.5*(cos(2*pi*t/tperiod) - 1.) + plot(tt, eta_tide(tt), 'k') + plot([tframe],[eta_tide(tframe)], 'ro') + grid(True) + title('Tide stage in ocean') + xticks(np.arange(0,5.1*24,24)) + xlabel('hours') + ylabel('meters (0=MTL)') + xlim(0,120) + ylim(-2,2) + +def plot_src_boxes(): + + if 0: + # Johns River + x1rs = -123.951 + x2rs = -123.950 + y1rs = 46.8815 + y2rs = 46.8825 + + if 0: + # Chehalis River: + x1rs = -123.743 + x2rs = -123.740 + y1rs = 46.956 + y2rs = 46.958 + kwargs = {'color':'yellow', 'linewidth':1} + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + # Wishkah River + x1rs = -123.811 + x2rs = -123.810 + y1rs = 46.990 + y2rs = 46.991 + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + # Hoquiam River + x1rs = -123.886 + x2rs = -123.885 + y1rs = 46.995 + y2rs = 46.996 + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.format = 'binary' # 'ascii' or 'binary' to match setrun.py + + + # To plot gauge locations on pcolor or contour plot, use this as + # an afteraxis function: + + def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=[0,1102,1187], format_string='ko', add_labels=True) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Domain', figno=50) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + plotaxes.xlimits = [-124.65,-123.65] + plotaxes.ylimits = [46.8, 47.1] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface', figno=0) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.xlimits = [-124.23,-123.65] + plotaxes.xlimits = [-124.25,-123.65] + plotaxes.ylimits = [46.82, 47.06] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + + + #----------------------------------------- + # Figure for speed + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Speed', figno=1) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.xlimits = [-124.23,-123.65] + plotaxes.xlimits = [-124.25,-123.65] + plotaxes.ylimits = [46.82, 47.06] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Speed at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = speed + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = cmin_speed + plotitem.pcolor_cmax = cmax_speed + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'max'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Aberdeen', figno=10) + plotfigure.kwargs = {'figsize':(11,8)} + #plotfigure.show = False + xlimits = [-123.89,-123.73] + ylimits = [46.94, 46.996] + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.5])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + plotaxes.xlimits = xlimits + plotaxes.ylimits = ylimits + + + def fixup(current_data): + import pylab + #addgauges(current_data) + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=range(54,58), format_string='ko', add_labels=True) + t = current_data.t + t = t / 3600. # hours + pylab.title('Aberdeen at t = %.2f hours' % t, fontsize=15) + pylab.xticks(rotation=20,fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + pylab.ticklabel_format(useOffset=False) + plot_src_boxes() + + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.7 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Aberdeen speed', figno=11) + plotfigure.kwargs = {'figsize':(11,8)} + #plotfigure.show = False + + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.5])' + plotaxes.title = 'Speed' + plotaxes.scaled = True + plotaxes.xlimits = xlimits + plotaxes.ylimits = ylimits + + def fixup(current_data): + import pylab + #addgauges(current_data) + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=range(54,58), format_string='ko', add_labels=True) + t = current_data.t + t = t / 3600. # hours + pylab.title('Aberdeen speed at t = %.2f hours' % t, fontsize=15) + pylab.xticks(rotation=20,fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + pylab.ticklabel_format(useOffset=False) + plot_src_boxes() + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = speed + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = cmin_speed + plotitem.pcolor_cmax = cmax_speed + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.7 + plotitem.colorbar_kwargs = {'extend':'max'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface at gauges', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Surface' + + # Plot surface as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 3 + plotitem.plotstyle = 'b-' + + # Plot topo as green curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.show = False + + def gaugetopo(current_data): + q = current_data.q + h = q[0,:] + eta = q[3,:] + topo = eta - h + return topo + + plotitem.plot_var = gaugetopo + plotitem.plotstyle = 'g-' + + + + #----------------------------------------- + # Plots of timing (CPU and wall time): + + def make_timing_plots(plotdata): + from clawpack.visclaw import plot_timing_stats + import os,sys + try: + timing_plotdir = plotdata.plotdir + '/_timing_figures' + os.system('mkdir -p %s' % timing_plotdir) + # adjust units for plots based on problem: + units = {'comptime':'seconds', 'simtime':'hours', + 'cell':'millions'} + plot_timing_stats.make_plots(outdir=plotdata.outdir, + make_pngs=True, + plotdir=timing_plotdir, + units=units) + except: + print('*** Error making timing plots') + + otherfigure = plotdata.new_otherfigure(name='timing plots', + fname='_timing_figures/timing.html') + otherfigure.makefig = make_timing_plots + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata + diff --git a/GraysHarborBC/kingtide2015/setprob.f90 b/GraysHarborBC/kingtide2015/setprob.f90 new file mode 100644 index 0000000..80b9f95 --- /dev/null +++ b/GraysHarborBC/kingtide2015/setprob.f90 @@ -0,0 +1,27 @@ +subroutine setprob() + + use setprob_module, only: a + + implicit none + + + INTEGER :: col,rows,io + + rows = 0 + OPEN(UNIT=19, FILE="tidal_signal.txt") + !count # of rows + DO + READ(19,*,iostat=io) + IF (io/=0) EXIT + rows = rows + 1 + END DO + + rewind(19) + !export as a() + DO col = 1, rows + READ(19,*) a(col) + END DO + close(19) + + +end subroutine setprob diff --git a/GraysHarborBC/kingtide2015/setprob_module.f90 b/GraysHarborBC/kingtide2015/setprob_module.f90 new file mode 100644 index 0000000..0a9deaf --- /dev/null +++ b/GraysHarborBC/kingtide2015/setprob_module.f90 @@ -0,0 +1,8 @@ + +module setprob_module + + real, DIMENSION(1201) :: a + + save + +end module setprob_module diff --git a/GraysHarborBC/kingtide2015/setrun.py b/GraysHarborBC/kingtide2015/setrun.py new file mode 100644 index 0000000..56cae0b --- /dev/null +++ b/GraysHarborBC/kingtide2015/setrun.py @@ -0,0 +1,465 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os, sys +import numpy as np +from clawpack.amrclaw.data import FlagRegion + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + # x values should be integer multipes of 1/9" + # y values should be integer multipes of 1/9" + # Note: always satisfied if limits are multiples of 0.01 degree + + clawdata.lower[0] = -124.65 + clawdata.upper[0] = -123.65 + clawdata.lower[1] = 46.8 + clawdata.upper[1] = 47.1 + + clawdata.num_cells[0] = 60 + clawdata.num_cells[1] = 18 + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 3*24 + clawdata.tfinal = 3*24*3600. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 1 + clawdata.output_t0 = True + + + clawdata.output_format = 'binary' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.8 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 0 + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 1 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 4 + + # List of refinement ratios at each level (length at least mxnest-1) + # 1', 12", 4", 1", 1/3" + amrdata.refinement_ratios_x = [5,3,4,3] + amrdata.refinement_ratios_y = [5,3,4,3] + amrdata.refinement_ratios_t = [5,3,4,3] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + regions = rundata.regiondata.regions + + regions.append([1,1,0.,1e9,clawdata.lower[0]-0.1,clawdata.upper[0]+0.1, + clawdata.lower[1]-0.1,clawdata.upper[1]+0.1]) + regions.append([1,2,0.,1e9,-124.25, -123.65, 46.8, 47.1]) + regions.append([1,3,0.,1e9,-124.2, -123.65, 46.8, 47.05]) + regions.append([2,4,0.,1e9,-124.119,-124.097,46.903, 46.915]) # Wp marina + regions.append([2,4,0.,1e9,-123.89,-123.73,46.94, 46.996]) # Aberdeen + + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + gauges = rundata.gaugedata.gauges + rundata.gaugedata.min_time_increment = 1. # seconds between gauge output + rundata.gaugedata.display_format = "f15.7" # need enough digits for u,v + + gauges.append([0,-124.4,46.92,0,1e9]) # Offshore tide + gauges.append([1,-124.11,46.92,0,1e9]) # GH inlet + gauges.append([2,-123.92,46.95,0,1e9]) + gauges.append([3,-123.85,46.96,0,1e9]) + gauges.append([4,-124.00,46.90,0,1e9]) # Johns River + gauges.append([5,-123.808,46.972,0,1e9]) # Chehalis at Wishkah + + gauges.append([1187,-123.8571, 46.9657, 0,1e9]) # Aberdeen + gauges.append([1102,-124.10485, 46.90447, 0,1e9]) # Westport + + gauges.append([8,-124.6,46.96,0,1e9]) # at Left Boundary + + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = -0.457 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.05 + refinement_data.speed_tolerance = 6*[0.05] + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topofiles = topo_data.topofiles + + topodir = '../topo/topofiles/' + topofiles.append([3, topodir + 'etopo1_-125_-123_46_48_1min.asc']) + topofiles.append([3, topodir+'GH_13sec_mtl.asc']) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, minlevel,maxlevel,fname] + dtopo_data.dtopofiles = [] + dtopo_data.dt_max_dtopo = 1.0 + + + # == setqinit.data values == + + + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # variable_eta_init newly added to QinitData: + rundata.qinit_data.variable_eta_init = False + + + # == fgmax.data values == + # MODIFIED FORMAT: + # Now specify a list of objects of class fgmax_tools.FGmaxGrid + # specifying any fgmax grids. + + rundata.fgmax_data.fgmax_grids = [] + rundata.fgmax_data.num_fgmax_val = 5 # Save depth, speed, hs, hss, hmin + + return rundata + # end of function setgeo + # ----------------------: + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + from clawpack.geoclaw import kmltools + rundata = setrun(*sys.argv[1:]) + rundata.write() + + kmltools.make_input_data_kmls(rundata) diff --git a/GraysHarborBC/kingtide2015/src/Riemann/bc2amr.f90 b/GraysHarborBC/kingtide2015/src/Riemann/bc2amr.f90 new file mode 100644 index 0000000..976d602 --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/Riemann/bc2amr.f90 @@ -0,0 +1,329 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + use setprob_module, only: a + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt, t_num + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + ! Getting the jump in tidal signal using linear Intepolation on height.txt + t_num = floor(time/360) + delt = possk(level) + jump_h = (a(t_num+2) - a(t_num+1))*delt/360 + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = 0 + val(1, i, j) = val(1, nxl + 1, j) + jump_h + h_r = val(1, nxl + 1, j) + h_m = val(1, i, j) + u_r = val(2, nxl + 1, j) + if (h_r < h_m) then + 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 + else + val(2, i, j) = (h_m*u_r - 2*h_r*(sqrt(9.81*h_r)-sqrt(9.81*h_m)))/h_m + endif + end do + end do + + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/GraysHarborBC/kingtide2015/src/Simple/bc2amr.f90 b/GraysHarborBC/kingtide2015/src/Simple/bc2amr.f90 new file mode 100644 index 0000000..78c36ee --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/Simple/bc2amr.f90 @@ -0,0 +1,324 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + use setprob_module, only: a + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt, t_num + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + ! Getting the jump in tidal signal using linear Intepolation on height.txt + t_num = floor(time/360) + delt = possk(level) + jump_h = (a(t_num+2) - a(t_num+1))*delt/360 + + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = 0 + val(1, i, j) = val(1, nxl + 1, j) + jump_h + val(2, i, j) = val(2, nxl + 1, j) + + end do + end do + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/GraysHarborBC/kingtide2015/src/ocean-forcing/Makefile b/GraysHarborBC/kingtide2015/src/ocean-forcing/Makefile new file mode 100644 index 0000000..137985a --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/ocean-forcing/Makefile @@ -0,0 +1,71 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + + +MODULES = \ + setprob_module.f90 \ + +SOURCES = \ + src2.f90 \ + setprob.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + python maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls + diff --git a/GraysHarborBC/kingtide2015/src/ocean-forcing/setprob.f90 b/GraysHarborBC/kingtide2015/src/ocean-forcing/setprob.f90 new file mode 100644 index 0000000..64bef71 --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/ocean-forcing/setprob.f90 @@ -0,0 +1,29 @@ +subroutine setprob() + + use setprob_module, only: tide_t, tide_etaprime, tide_m + + implicit none + + real(kind=8) :: tide_eta ! read in but not used + + integer iunit, i + character(len=25) fname + + !if not using tide data... + !return + + iunit = 7 + fname = 'tidedata_dec2015.txt' + open(unit=iunit, file=fname, status='old') + + read(iunit,*) tide_m + allocate(tide_t(tide_m), tide_etaprime(tide_m)) + + do i=1,tide_m + read(iunit,*) tide_t(i), tide_eta, tide_etaprime(i) + enddo + + close(iunit) + + +end subroutine setprob diff --git a/GraysHarborBC/kingtide2015/src/ocean-forcing/setprob_module.f90 b/GraysHarborBC/kingtide2015/src/ocean-forcing/setprob_module.f90 new file mode 100644 index 0000000..36897a4 --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/ocean-forcing/setprob_module.f90 @@ -0,0 +1,9 @@ + +module setprob_module + + integer :: tide_m + real(kind=8), allocatable, dimension(:) :: tide_t, tide_etaprime + + save + +end module setprob_module diff --git a/GraysHarborBC/kingtide2015/src/ocean-forcing/setrun.py b/GraysHarborBC/kingtide2015/src/ocean-forcing/setrun.py new file mode 100644 index 0000000..226474a --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/ocean-forcing/setrun.py @@ -0,0 +1,464 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os, sys +import numpy as np +from clawpack.amrclaw.data import FlagRegion + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + # x values should be integer multipes of 1/9" + # y values should be integer multipes of 1/9" + # Note: always satisfied if limits are multiples of 0.01 degree + + clawdata.lower[0] = -124.65 + clawdata.upper[0] = -123.65 + clawdata.lower[1] = 46.8 + clawdata.upper[1] = 47.1 + + clawdata.num_cells[0] = 60 + clawdata.num_cells[1] = 18 + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 3*24 + clawdata.tfinal = 3*24*3600. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 1 + clawdata.output_t0 = True + + + clawdata.output_format = 'binary' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.8 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 1 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 4 + + # List of refinement ratios at each level (length at least mxnest-1) + # 1', 12", 4", 1", 1/3" + amrdata.refinement_ratios_x = [5,3,4,3] + amrdata.refinement_ratios_y = [5,3,4,3] + amrdata.refinement_ratios_t = [5,3,4,3] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + regions = rundata.regiondata.regions + + regions.append([1,1,0.,1e9,clawdata.lower[0]-0.1,clawdata.upper[0]+0.1, + clawdata.lower[1]-0.1,clawdata.upper[1]+0.1]) + regions.append([1,2,0.,1e9,-124.25, -123.65, 46.8, 47.1]) + regions.append([1,3,0.,1e9,-124.2, -123.65, 46.8, 47.05]) + regions.append([2,4,0.,1e9,-124.119,-124.097,46.903, 46.915]) # Wp marina + regions.append([2,4,0.,1e9,-123.89,-123.73,46.94, 46.996]) # Aberdeen + + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + gauges = rundata.gaugedata.gauges + rundata.gaugedata.min_time_increment = 1. # seconds between gauge output + rundata.gaugedata.display_format = "f15.7" # need enough digits for u,v + + gauges.append([0,-124.4,46.92,0,1e9]) # Offshore tide + gauges.append([1,-124.11,46.92,0,1e9]) # GH inlet + gauges.append([2,-123.92,46.95,0,1e9]) + gauges.append([3,-123.85,46.96,0,1e9]) + gauges.append([4,-124.00,46.90,0,1e9]) # Johns River + gauges.append([5,-123.808,46.972,0,1e9]) # Chehalis at Wishkah + + gauges.append([1187,-123.8571, 46.9657, 0,1e9]) # Aberdeen + gauges.append([1102,-124.10485, 46.90447, 0,1e9]) # Westport + + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = -0.468 + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.05 + refinement_data.speed_tolerance = 6*[0.05] + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topofiles = topo_data.topofiles + + topodir = '../topo/topofiles/' + topofiles.append([3, topodir + 'etopo1_-125_-123_46_48_1min.asc']) + topofiles.append([3, topodir+'GH_13sec_mtl.asc']) + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, minlevel,maxlevel,fname] + dtopo_data.dtopofiles = [] + dtopo_data.dt_max_dtopo = 1.0 + + + # == setqinit.data values == + + + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # variable_eta_init newly added to QinitData: + rundata.qinit_data.variable_eta_init = False + + + # == fgmax.data values == + # MODIFIED FORMAT: + # Now specify a list of objects of class fgmax_tools.FGmaxGrid + # specifying any fgmax grids. + + rundata.fgmax_data.fgmax_grids = [] + rundata.fgmax_data.num_fgmax_val = 5 # Save depth, speed, hs, hss, hmin + + return rundata + # end of function setgeo + # ----------------------: + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + from clawpack.geoclaw import kmltools + rundata = setrun(*sys.argv[1:]) + rundata.write() + + kmltools.make_input_data_kmls(rundata) + diff --git a/GraysHarborBC/kingtide2015/src/ocean-forcing/src2.f90 b/GraysHarborBC/kingtide2015/src/ocean-forcing/src2.f90 new file mode 100644 index 0000000..2f8038e --- /dev/null +++ b/GraysHarborBC/kingtide2015/src/ocean-forcing/src2.f90 @@ -0,0 +1,312 @@ +subroutine src2(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt) + + use geoclaw_module, only: g => grav, coriolis_forcing, coriolis + use geoclaw_module, only: friction_forcing, friction_depth + use geoclaw_module, only: manning_coefficient + use geoclaw_module, only: manning_break, num_manning + use geoclaw_module, only: spherical_distance, coordinate_system + use geoclaw_module, only: RAD2DEG, pi, dry_tolerance, DEG2RAD + use geoclaw_module, only: rho_air + + use storm_module, only: wind_forcing, pressure_forcing, wind_drag + use storm_module, only: wind_index, pressure_index + use storm_module, only: storm_direction, storm_location + + use friction_module, only: variable_friction, friction_index + + use setprob_module, only: tide_t, tide_etaprime, tide_m + + implicit none + + ! Input parameters + integer, intent(in) :: meqn,mbc,mx,my,maux + double precision, intent(in) :: xlower,ylower,dx,dy,t,dt + + ! Output + double precision, intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc) + double precision, intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + + ! Locals + integer :: i, j, nman, k + real(kind=8) :: h, hu, hv, gamma, dgamma, y, fdt, a(2,2), coeff + real(kind=8) :: xm, xc, xp, ym, yc, yp, dx_meters, dy_meters + real(kind=8) :: u, v, hu0, hv0 + real(kind=8) :: tau, wind_speed, theta, phi, psi, P_gradient(2), S(2) + real(kind=8) :: Ddt, sloc(2) + real(kind=8) :: eta,x,tperiod,deta_dt,ramp,alpha,ampl + real(kind=8) :: x1rs,x2rs,y1rs,y2rs,area,discharge_cfs,discharge,rs + + ! Algorithm parameters + ! Parameter controls when to zero out the momentum at a depth in the + ! friction source term + real(kind=8), parameter :: depth_tolerance = 1.0d-30 + + ! Physics + ! Nominal density of water + real(kind=8), parameter :: rho = 1025.d0 + + + ! Tidal forcing + ! This is applied as a source term based on desired eta'(t) offshore + + if (.false.) then + ! Sinusoidal forcing as a test: + tperiod = 12*3600. + ampl = 1.d0 + eta = ampl*sin(2.d0*pi*t/tperiod) + deta_dt = ampl*2.d0*pi/tperiod * cos(2.d0*pi*t/tperiod) + endif + + if (.true.) then + ! tidal forcing read in by setprob.f90: + do k=1,tide_m + if (t < tide_t(k)) then + exit + endif + enddo + + if ((k>1) .and. (k <= tide_m)) then + alpha = (t - tide_t(k-1)) / (tide_t(k) - tide_t(k-1)) + deta_dt = (1.d0-alpha)*tide_etaprime(k-1) + alpha*tide_etaprime(k) + else + deta_dt = 0.d0 + endif + endif + + do i=1-mbc,mx+mbc + x = xlower + (i-0.5d0)*dx + if (x < -124.3d0) then + do j=1-mbc,my+mbc + q(1,i,j) = q(1,i,j) + dt*deta_dt + enddo + else if (x < -124.19d0) then + ! ramp down to zero closer to shore: + ramp = 1.d0 - (x + 124.3d0)/(-124.19d0 + 124.3d0) + do j=1-mbc,my+mbc + q(1,i,j) = q(1,i,j) + dt*deta_dt*ramp + enddo + endif + enddo + + + ! River sources: + + if (.false.) then + + ! Johns River + x1rs = -123.951d0 + x2rs = -123.950d0 + y1rs = 46.8815d0 + y2rs = 46.8825d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 5000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + endif + + if (.false.) then + + ! Chehalis + + x1rs = -123.743d0 + x2rs = -123.740d0 + y1rs = 46.956d0 + y2rs = 46.958d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 20000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + + ! Wishkah River + x1rs = -123.811d0 + x2rs = -123.810d0 + y1rs = 46.990d0 + y2rs = 46.991d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 4000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + + + ! Hoquiam River + x1rs = -123.886d0 + x2rs = -123.885d0 + y1rs = 46.995d0 + y2rs = 46.996d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 9000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + + endif + + ! Friction source term + if (friction_forcing) then + do j=1,my + do i=1,mx + ! Extract appropriate momentum + if (q(1,i,j) < depth_tolerance) then + q(2:3,i,j) = 0.d0 + else + ! Apply friction source term only if in shallower water + if (q(1,i,j) <= friction_depth) then + if (.not.variable_friction) then + do nman = num_manning, 1, -1 + if (aux(1,i,j) .lt. manning_break(nman)) then + coeff = manning_coefficient(nman) + endif + enddo + else + coeff = aux(friction_index,i,j) + endif + + ! Calculate source term + gamma = sqrt(q(2,i,j)**2 + q(3,i,j)**2) * g & + * coeff**2 / (q(1,i,j)**(7.d0/3.d0)) + dgamma = 1.d0 + dt * gamma + q(2, i, j) = q(2, i, j) / dgamma + q(3, i, j) = q(3, i, j) / dgamma + endif + endif + enddo + enddo + endif + ! End of friction source term + + ! Coriolis source term + ! TODO: May want to remove the internal calls to coriolis as this could + ! lead to slow downs. + if (coriolis_forcing) then + do j=1,my + y = ylower + (j - 0.5d0) * dy + fdt = coriolis(y) * dt ! Calculate f dependent on coordinate system + + ! Calculate matrix components + a(1,1) = 1.d0 - 0.5d0 * fdt**2 + fdt**4 / 24.d0 + a(1,2) = fdt - fdt**3 / 6.d0 + a(2,1) = -fdt + fdt**3 / 6.d0 + a(2,2) = a(1,1) + + do i=1,mx + q(2,i,j) = q(2, i, j) * a(1,1) + q(3, i, j) * a(1,2) + q(3,i,j) = q(2, i, j) * a(2,1) + q(3, i, j) * a(2,2) + enddo + enddo + endif + ! End of coriolis source term + + ! wind ----------------------------------------------------------- + if (wind_forcing) then + ! Need storm location and direction for sector based wind drag + sloc = storm_location(t) + theta = storm_direction(t) + do j=1,my + yc = ylower + (j - 0.5d0) * dy + do i=1,mx + xc = xlower + (i - 0.5d0) * dx + if (q(1,i,j) > dry_tolerance) then + psi = atan2(yc - sloc(2), xc - sloc(1)) + if (theta > psi) then + phi = (2.d0 * pi - theta + psi) * RAD2DEG + else + phi = (psi - theta) * RAD2DEG + endif + wind_speed = sqrt(aux(wind_index,i,j)**2 & + + aux(wind_index+1,i,j)**2) + tau = wind_drag(wind_speed, phi) * rho_air * wind_speed / rho + q(2,i,j) = q(2,i,j) + dt * tau * aux(wind_index,i,j) + q(3,i,j) = q(3,i,j) + dt * tau * aux(wind_index+1,i,j) + endif + enddo + enddo + endif + ! ---------------------------------------------------------------- + + ! Atmosphere Pressure -------------------------------------------- + ! Handled in Riemann solver + ! if (pressure_forcing) then + ! do j=1,my + ! ym = ylower + (j - 1.d0) * dy + ! yc = ylower + (j - 0.5d0) * dy + ! yp = ylower + j * dy + ! do i=1,mx + ! xm = xlower + (i - 1.d0) * dx + ! xc = xlower + (i - 0.5d0) * dx + ! xp = xlower + i * dx + + ! if (coordinate_system == 2) then + ! ! Convert distance in lat-long to meters + ! dx_meters = spherical_distance(xp,yc,xm,yc) + ! dy_meters = spherical_distance(xc,yp,xc,ym) + ! else + ! dx_meters = dx + ! dy_meters = dy + ! endif + + ! ! Extract depths + ! h = q(1,i,j) + + ! ! Calculate gradient of Pressure + ! P_gradient(1) = (aux(pressure_index,i+1,j) & + ! - aux(pressure_index,i-1,j)) / (2.d0 * dx_meters) + ! P_gradient(2) = (aux(pressure_index,i,j+1) & + ! - aux(pressure_index,i,j-1)) / (2.d0 * dy_meters) + + ! ! Modify momentum in each layer + ! if (h > dry_tolerance) then + ! q(2, i, j) = q(2, i, j) - dt * h * P_gradient(1) / rho + ! q(3, i, j) = q(3, i, j) - dt * h * P_gradient(2) / rho + ! end if + ! enddo + ! enddo + ! endif + +end subroutine src2 diff --git a/GraysHarborBC/sine/GraysSine.ipynb b/GraysHarborBC/sine/GraysSine.ipynb new file mode 100644 index 0000000..532f606 --- /dev/null +++ b/GraysHarborBC/sine/GraysSine.ipynb @@ -0,0 +1,132 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The tidal signal $\\sin(2\\pi t / T)$ is implemented, a pure sine wave with amplitude 1 and period $T = 12\\times 3600 = $ 12 hours .\n", + "\n", + "For Riemann and Simple case, it's implemented on the left boundary, for the ocean forcing case, it is forced inside the domain.\n", + "\n", + "(This notebook is run for the Simple case. Only the title of the graph needs to be changed in the code for Riemann and ocean forcing case.)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pylab import *" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from IPython.display import Image\n", + "import clawpack.pyclaw.gauges as gauges" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "figure(400, figsize=(13,6))\n", + "clf()\n", + "colors = ['k','b','r']\n", + "\n", + "outdir = '_output'\n", + "\n", + "for k,gaugeno in enumerate([1102,1187]):\n", + " gauge = gauges.GaugeSolution(gaugeno, outdir)\n", + " t = gauge.t / 3600. # convert to hours\n", + " q = gauge.q\n", + " eta = q[3,:]\n", + " plot(t, eta, colors[k], label='Gauge %s' % gaugeno)\n", + "\n", + " # determine amplification and time shift:\n", + " m2 = int(floor(0.75*len(eta)))\n", + " eta2 = eta[m2:] # last part of eta signal\n", + " etamax2 = eta2.max()\n", + " etamin2 = eta2.min()\n", + " t2 = t[m2:]\n", + " jtmax = argmax(eta2) \n", + " tshift = (t2[jtmax] - 39.)*3600.\n", + " \n", + " print('At gauge %i, etamin2 = %.3f, etamax2 = %.3f at tshift = %.1f s' \\\n", + " % (gaugeno,etamin2,etamax2,tshift))\n", + "\n", + "tperiod = 12\n", + "eta = 1.*sin(2*pi*t/tperiod)\n", + "plot(t, eta, 'k--', label='Sine')\n", + "\n", + "legend(loc='upper right')\n", + "xlabel('hours')\n", + "ylabel('Surface relative to MTL (m)')\n", + "grid(True)\n", + "title('Simple sine condition and resulting GeoClaw gauge results');\n", + "\n", + "xticks(arange(0,t[-1]+0.1,12))\n", + "xlim(0,60)\n", + "\n", + "if 0:\n", + " fname = 'GaugeComparison.png'\n", + " savefig(fname, bbox_inches='tight')\n", + " print('Created %s' % fname) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Notes:\n", + "\n", + "- Given an observed tide at Westport that we want to match, we can shift it by about 3220 seconds and increase it by a factor 1.04 in order to obtain the signal to use at the left boundary. This is illustrated in the example `../kingtide2015`.\n" + ] + }, + { + "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 +} diff --git a/GraysHarborBC/sine/Makefile b/GraysHarborBC/sine/Makefile new file mode 100644 index 0000000..f8d638a --- /dev/null +++ b/GraysHarborBC/sine/Makefile @@ -0,0 +1,66 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + +SOURCES = \ + bc2amr.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + python maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls + diff --git a/GraysHarborBC/sine/README.md b/GraysHarborBC/sine/README.md new file mode 100644 index 0000000..b7bb4c2 --- /dev/null +++ b/GraysHarborBC/sine/README.md @@ -0,0 +1,7 @@ +The current directory is set up for the Simple method. + +First run the GeoClaw code and then use this notebook to examine the gauge results. + +To use Riemann method, replace ./bc2amr.f90 by ./src/Riemann/bc2amr.f90 . + +To use ocean forcing method, go to ./src/ocean-forcing, move every file to the current directory. diff --git a/GraysHarborBC/sine/bc2amr.f90 b/GraysHarborBC/sine/bc2amr.f90 new file mode 100644 index 0000000..6c729f2 --- /dev/null +++ b/GraysHarborBC/sine/bc2amr.f90 @@ -0,0 +1,323 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt, t_num + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + + ! Sine boundary condition + delt = possk(level) + jump_h = sin(2.d0*pi*(time+delt)/(12*3600)) - sin(2.d0*pi*(time)/(12*3600)) + + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = 0 + val(1, i, j) = val(1, nxl + 1, j) + jump_h + val(2, i, j) = val(2, nxl + 1, j) + end do + end do + + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/GraysHarborBC/sine/setplot.py b/GraysHarborBC/sine/setplot.py new file mode 100644 index 0000000..6f21eee --- /dev/null +++ b/GraysHarborBC/sine/setplot.py @@ -0,0 +1,538 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d + +from clawpack.geoclaw import topotools +from clawpack.visclaw import gaugetools, plottools, colormaps + +cmin_water = -1. +cmax_water = 1. +cmin_speed = 0. +cmax_speed = 2. + +def speed(current_data): + """ + Return a masked array containing the speed only in wet cells. + Surface is eta = h+topo, assumed to be output as 4th column of fort.q + files. + """ + from clawpack.visclaw.geoplot import drytol_default + drytol = getattr(current_data.user, 'drytol', drytol_default) + q = current_data.q + h = q[0,:,:] + hu = q[1,:,:] + hv = q[2,:,:] + eta = q[3,:,:] + + hwet = np.ma.masked_where(h <= drytol, h) + speed = np.sqrt(hu**2 + hv**2) / hwet # will be masked where dry + + try: + # Use mask covering coarse regions if it's set: + m = current_data.mask_coarse + speed = numpy.ma.masked_where(m, speed) + except: + pass + + return speed + +if 0: + sea_level = -0.468 + tide_data = np.loadtxt('tidedata_dec2015.txt', skiprows=1) + t_tide = tide_data[:,0] / 3600. + #offset = 2.561 # from MLLW to MHW + offset = 0. + eta_tide = interp1d(t_tide, tide_data[:,1]-offset, bounds_error=False, + fill_value = np.nan) + +if 1: + # forcing: + sea_level = 0. + tperiod = 12 # hours + eta_tide = lambda t: 1.*np.sin(2*np.pi*t/tperiod) + t_tide = np.linspace(0,48,1000) + + + +def plot_tide(current_data): + from pylab import linspace,plot,grid,title,cos,pi,xlabel,ylabel,xticks + from pylab import xlim,ylim + tframe = current_data.t / 3600. + tt = linspace(t_tide[0],t_tide[-1],1000) + plot(tt, eta_tide(tt), 'k') + plot([tframe],[eta_tide(tframe)], 'ro') + grid(True) + title('Tide stage in ocean') + xticks(np.arange(0,5.1*24,24)) + xlabel('hours') + ylabel('meters (0=MTL)') + xlim(0,48) + ylim(-1.2,1.2) + +def plot_src_boxes(): + + if 0: + # Johns River + x1rs = -123.951 + x2rs = -123.950 + y1rs = 46.8815 + y2rs = 46.8825 + + if 0: + # Chehalis River: + x1rs = -123.743 + x2rs = -123.740 + y1rs = 46.956 + y2rs = 46.958 + kwargs = {'color':'yellow', 'linewidth':1} + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + # Wishkah River + x1rs = -123.811 + x2rs = -123.810 + y1rs = 46.990 + y2rs = 46.991 + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + # Hoquiam River + x1rs = -123.886 + x2rs = -123.885 + y1rs = 46.995 + y2rs = 46.996 + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.format = 'binary' # 'ascii' or 'binary' to match setrun.py + + + # To plot gauge locations on pcolor or contour plot, use this as + # an afteraxis function: + + def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=[0,1102,1187], format_string='ko', add_labels=True) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Domain', figno=50) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + plotaxes.xlimits = [-124.65,-123.65] + plotaxes.ylimits = [46.8, 47.1] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface', figno=0) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.xlimits = [-124.23,-123.65] + plotaxes.xlimits = [-124.25,-123.65] + plotaxes.ylimits = [46.82, 47.06] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + + + #----------------------------------------- + # Figure for speed + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Speed', figno=1) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.xlimits = [-124.23,-123.65] + plotaxes.xlimits = [-124.25,-123.65] + plotaxes.ylimits = [46.82, 47.06] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Speed at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = speed + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = cmin_speed + plotitem.pcolor_cmax = cmax_speed + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'max'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Aberdeen', figno=10) + plotfigure.kwargs = {'figsize':(11,8)} + #plotfigure.show = False + xlimits = [-123.89,-123.73] + ylimits = [46.94, 46.996] + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.5])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + plotaxes.xlimits = xlimits + plotaxes.ylimits = ylimits + + + def fixup(current_data): + import pylab + #addgauges(current_data) + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=range(54,58), format_string='ko', add_labels=True) + t = current_data.t + t = t / 3600. # hours + pylab.title('Aberdeen at t = %.2f hours' % t, fontsize=15) + pylab.xticks(rotation=20,fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + pylab.ticklabel_format(useOffset=False) + plot_src_boxes() + + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.7 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Aberdeen speed', figno=11) + plotfigure.kwargs = {'figsize':(11,8)} + #plotfigure.show = False + + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.5])' + plotaxes.title = 'Speed' + plotaxes.scaled = True + plotaxes.xlimits = xlimits + plotaxes.ylimits = ylimits + + def fixup(current_data): + import pylab + #addgauges(current_data) + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=range(54,58), format_string='ko', add_labels=True) + t = current_data.t + t = t / 3600. # hours + pylab.title('Aberdeen speed at t = %.2f hours' % t, fontsize=15) + pylab.xticks(rotation=20,fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + pylab.ticklabel_format(useOffset=False) + plot_src_boxes() + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = speed + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = cmin_speed + plotitem.pcolor_cmax = cmax_speed + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.7 + plotitem.colorbar_kwargs = {'extend':'max'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface at gauges', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Surface' + + # Plot surface as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 3 + plotitem.plotstyle = 'b-' + + # Plot topo as green curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.show = False + + def gaugetopo(current_data): + q = current_data.q + h = q[0,:] + eta = q[3,:] + topo = eta - h + return topo + + plotitem.plot_var = gaugetopo + plotitem.plotstyle = 'g-' + + + + #----------------------------------------- + # Plots of timing (CPU and wall time): + + def make_timing_plots(plotdata): + from clawpack.visclaw import plot_timing_stats + import os,sys + try: + timing_plotdir = plotdata.plotdir + '/_timing_figures' + os.system('mkdir -p %s' % timing_plotdir) + # adjust units for plots based on problem: + units = {'comptime':'seconds', 'simtime':'hours', + 'cell':'millions'} + plot_timing_stats.make_plots(outdir=plotdata.outdir, + make_pngs=True, + plotdir=timing_plotdir, + units=units) + except: + print('*** Error making timing plots') + + otherfigure = plotdata.new_otherfigure(name='timing plots', + fname='_timing_figures/timing.html') + otherfigure.makefig = make_timing_plots + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata + diff --git a/GraysHarborBC/sine/setrun.py b/GraysHarborBC/sine/setrun.py new file mode 100644 index 0000000..0f3046a --- /dev/null +++ b/GraysHarborBC/sine/setrun.py @@ -0,0 +1,467 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os, sys +import numpy as np +from clawpack.amrclaw.data import FlagRegion + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + # x values should be integer multipes of 1/9" + # y values should be integer multipes of 1/9" + # Note: always satisfied if limits are multiples of 0.01 degree + + clawdata.lower[0] = -124.65 + clawdata.upper[0] = -123.65 + clawdata.lower[1] = 46.8 + clawdata.upper[1] = 47.1 + + clawdata.num_cells[0] = 60 + clawdata.num_cells[1] = 18 + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 2*24 + clawdata.tfinal = 2*24*3600. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 1 + clawdata.output_t0 = True + + + clawdata.output_format = 'binary' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = 1 + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.8 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 0 + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 1 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least mxnest-1) + # 1', 12", 4", 1", 1/3" + amrdata.refinement_ratios_x = [5,3,4,3] + amrdata.refinement_ratios_y = [5,3,4,3] + amrdata.refinement_ratios_t = [5,3,4,3] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + regions = rundata.regiondata.regions + + regions.append([1,1,0.,1e9,clawdata.lower[0]-0.1,clawdata.upper[0]+0.1, + clawdata.lower[1]-0.1,clawdata.upper[1]+0.1]) + regions.append([1,2,0.,1e9,-124.25, -123.65, 46.8, 47.1]) + regions.append([1,3,0.,1e9,-124.2, -123.65, 46.8, 47.05]) + regions.append([2,4,0.,1e9,-124.119,-124.097,46.903, 46.915]) # Wp marina + regions.append([2,4,0.,1e9,-123.89,-123.73,46.94, 46.996]) # Aberdeen + + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + gauges = rundata.gaugedata.gauges + rundata.gaugedata.min_time_increment = 1. # seconds between gauge output + rundata.gaugedata.display_format = "f15.7" # need enough digits for u,v + + gauges.append([0,-124.4,46.92,0,1e9]) # Offshore tide + gauges.append([1,-124.11,46.92,0,1e9]) # GH inlet + gauges.append([2,-123.92,46.95,0,1e9]) + gauges.append([3,-123.85,46.96,0,1e9]) + gauges.append([4,-124.00,46.90,0,1e9]) # Johns River + gauges.append([5,-123.808,46.972,0,1e9]) # Chehalis at Wishkah + + gauges.append([8,-124.65,46.96,0,1e9]) # at Left Boundary + gauges.append([9,-124.65,46.97,0,1e9]) + + gauges.append([1187,-123.8571, 46.9657, 0,1e9]) # Aberdeen + gauges.append([1102,-124.10485, 46.90447, 0,1e9]) # Westport + + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0. + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.05 + refinement_data.speed_tolerance = 6*[0.05] + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topofiles = topo_data.topofiles + + topodir = '../topo/topofiles/' + topofiles.append([3, topodir + 'etopo1_-125_-123_46_48_1min.asc']) + topofiles.append([3, topodir+'GH_13sec_mtl.asc']) + + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, minlevel,maxlevel,fname] + dtopo_data.dtopofiles = [] + dtopo_data.dt_max_dtopo = 1.0 + + + # == setqinit.data values == + + + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # variable_eta_init newly added to QinitData: + rundata.qinit_data.variable_eta_init = False + + + # == fgmax.data values == + # MODIFIED FORMAT: + # Now specify a list of objects of class fgmax_tools.FGmaxGrid + # specifying any fgmax grids. + + rundata.fgmax_data.fgmax_grids = [] + rundata.fgmax_data.num_fgmax_val = 5 # Save depth, speed, hs, hss, hmin + + return rundata + # end of function setgeo + # ----------------------: + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + from clawpack.geoclaw import kmltools + rundata = setrun(*sys.argv[1:]) + rundata.write() + + kmltools.make_input_data_kmls(rundata) diff --git a/GraysHarborBC/sine/src/Riemann/bc2amr.f90 b/GraysHarborBC/sine/src/Riemann/bc2amr.f90 new file mode 100644 index 0000000..82174ca --- /dev/null +++ b/GraysHarborBC/sine/src/Riemann/bc2amr.f90 @@ -0,0 +1,329 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + !jump_h and jump in time defined + delt = possk(level) + jump_h = sin(2.d0*pi*(time+delt)/(12*3600)) - sin(2.d0*pi*(time)/(12*3600)) + + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = 0 + val(1, i, j) = val(1, nxl + 1, j) + jump_h + h_r = val(1, nxl + 1, j) + h_m = val(1, i, j) + u_r = val(2, nxl + 1, j) + if (h_r < h_m) then + 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 + else + val(2, i, j) = (h_m*u_r - 2*h_r*(sqrt(9.81*h_r)-sqrt(9.81*h_m)))/h_m + endif + + end do + end do + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/GraysHarborBC/sine/src/Simple/bc2amr.f90 b/GraysHarborBC/sine/src/Simple/bc2amr.f90 new file mode 100644 index 0000000..6c729f2 --- /dev/null +++ b/GraysHarborBC/sine/src/Simple/bc2amr.f90 @@ -0,0 +1,323 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt, t_num + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + + ! Sine boundary condition + delt = possk(level) + jump_h = sin(2.d0*pi*(time+delt)/(12*3600)) - sin(2.d0*pi*(time)/(12*3600)) + + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = 0 + val(1, i, j) = val(1, nxl + 1, j) + jump_h + val(2, i, j) = val(2, nxl + 1, j) + end do + end do + + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + stop "A user defined boundary condition was not provided." + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/GraysHarborBC/sine/src/ocean-forcing/Makefile b/GraysHarborBC/sine/src/ocean-forcing/Makefile new file mode 100644 index 0000000..137985a --- /dev/null +++ b/GraysHarborBC/sine/src/ocean-forcing/Makefile @@ -0,0 +1,71 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + + +MODULES = \ + setprob_module.f90 \ + +SOURCES = \ + src2.f90 \ + setprob.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + python maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls + diff --git a/GraysHarborBC/sine/src/ocean-forcing/setprob.f90 b/GraysHarborBC/sine/src/ocean-forcing/setprob.f90 new file mode 100644 index 0000000..2ecd212 --- /dev/null +++ b/GraysHarborBC/sine/src/ocean-forcing/setprob.f90 @@ -0,0 +1,30 @@ +subroutine setprob() + + use setprob_module, only: tide_t, tide_etaprime, tide_m + + implicit none + + real(kind=8) :: tide_eta ! read in but not used + + integer iunit, i + character(len=25) fname + + ! if not using tide data... + return + + + iunit = 7 + fname = 'XXX.txt' + open(unit=iunit, file=fname, status='old') + + read(iunit,*) tide_m + allocate(tide_t(tide_m), tide_etaprime(tide_m)) + + do i=1,tide_m + read(iunit,*) tide_t(i), tide_eta, tide_etaprime(i) + enddo + + close(iunit) + + +end subroutine setprob diff --git a/GraysHarborBC/sine/src/ocean-forcing/setprob_module.f90 b/GraysHarborBC/sine/src/ocean-forcing/setprob_module.f90 new file mode 100644 index 0000000..36897a4 --- /dev/null +++ b/GraysHarborBC/sine/src/ocean-forcing/setprob_module.f90 @@ -0,0 +1,9 @@ + +module setprob_module + + integer :: tide_m + real(kind=8), allocatable, dimension(:) :: tide_t, tide_etaprime + + save + +end module setprob_module diff --git a/GraysHarborBC/sine/src/ocean-forcing/setrun.py b/GraysHarborBC/sine/src/ocean-forcing/setrun.py new file mode 100644 index 0000000..3f8cc7e --- /dev/null +++ b/GraysHarborBC/sine/src/ocean-forcing/setrun.py @@ -0,0 +1,465 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os, sys +import numpy as np +from clawpack.amrclaw.data import FlagRegion + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + # x values should be integer multipes of 1/9" + # y values should be integer multipes of 1/9" + # Note: always satisfied if limits are multiples of 0.01 degree + + clawdata.lower[0] = -124.65 + clawdata.upper[0] = -123.65 + clawdata.lower[1] = 46.8 + clawdata.upper[1] = 47.1 + + clawdata.num_cells[0] = 60 + clawdata.num_cells[1] = 18 + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 2*24 + clawdata.tfinal = 2*24*3600. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 1 + clawdata.output_t0 = True + + + clawdata.output_format = 'binary' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.8 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 'extrap' + clawdata.bc_upper[0] = 'extrap' + + clawdata.bc_lower[1] = 'extrap' + clawdata.bc_upper[1] = 'extrap' + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 1 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least mxnest-1) + # 1', 12", 4", 1", 1/3" + amrdata.refinement_ratios_x = [5,3,4,3] + amrdata.refinement_ratios_y = [5,3,4,3] + amrdata.refinement_ratios_t = [5,3,4,3] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + regions = rundata.regiondata.regions + + regions.append([1,1,0.,1e9,clawdata.lower[0]-0.1,clawdata.upper[0]+0.1, + clawdata.lower[1]-0.1,clawdata.upper[1]+0.1]) + regions.append([1,2,0.,1e9,-124.25, -123.65, 46.8, 47.1]) + regions.append([1,3,0.,1e9,-124.2, -123.65, 46.8, 47.05]) + regions.append([2,4,0.,1e9,-124.119,-124.097,46.903, 46.915]) # Wp marina + regions.append([2,4,0.,1e9,-123.89,-123.73,46.94, 46.996]) # Aberdeen + + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + gauges = rundata.gaugedata.gauges + rundata.gaugedata.min_time_increment = 1. # seconds between gauge output + rundata.gaugedata.display_format = "f15.7" # need enough digits for u,v + + gauges.append([0,-124.4,46.92,0,1e9]) # Offshore tide + gauges.append([1,-124.11,46.92,0,1e9]) # GH inlet + gauges.append([2,-123.92,46.95,0,1e9]) + gauges.append([3,-123.85,46.96,0,1e9]) + gauges.append([4,-124.00,46.90,0,1e9]) # Johns River + gauges.append([5,-123.808,46.972,0,1e9]) # Chehalis at Wishkah + + gauges.append([1187,-123.8571, 46.9657, 0,1e9]) # Aberdeen + gauges.append([1102,-124.10485, 46.90447, 0,1e9]) # Westport + + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0. + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.05 + refinement_data.speed_tolerance = 6*[0.05] + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topofiles = topo_data.topofiles + + topodir = '../topo/topofiles/' + topofiles.append([3, topodir + 'etopo1_-125_-123_46_48_1min.asc']) + topofiles.append([3, topodir+'GH_13sec_mtl.asc']) + + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, minlevel,maxlevel,fname] + dtopo_data.dtopofiles = [] + dtopo_data.dt_max_dtopo = 1.0 + + + # == setqinit.data values == + + + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # variable_eta_init newly added to QinitData: + rundata.qinit_data.variable_eta_init = False + + + # == fgmax.data values == + # MODIFIED FORMAT: + # Now specify a list of objects of class fgmax_tools.FGmaxGrid + # specifying any fgmax grids. + + rundata.fgmax_data.fgmax_grids = [] + rundata.fgmax_data.num_fgmax_val = 5 # Save depth, speed, hs, hss, hmin + + return rundata + # end of function setgeo + # ----------------------: + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + from clawpack.geoclaw import kmltools + rundata = setrun(*sys.argv[1:]) + rundata.write() + + kmltools.make_input_data_kmls(rundata) + diff --git a/GraysHarborBC/sine/src/ocean-forcing/src2.f90 b/GraysHarborBC/sine/src/ocean-forcing/src2.f90 new file mode 100644 index 0000000..df3467e --- /dev/null +++ b/GraysHarborBC/sine/src/ocean-forcing/src2.f90 @@ -0,0 +1,312 @@ +subroutine src2(meqn,mbc,mx,my,xlower,ylower,dx,dy,q,maux,aux,t,dt) + + use geoclaw_module, only: g => grav, coriolis_forcing, coriolis + use geoclaw_module, only: friction_forcing, friction_depth + use geoclaw_module, only: manning_coefficient + use geoclaw_module, only: manning_break, num_manning + use geoclaw_module, only: spherical_distance, coordinate_system + use geoclaw_module, only: RAD2DEG, pi, dry_tolerance, DEG2RAD + use geoclaw_module, only: rho_air + + use storm_module, only: wind_forcing, pressure_forcing, wind_drag + use storm_module, only: wind_index, pressure_index + use storm_module, only: storm_direction, storm_location + + use friction_module, only: variable_friction, friction_index + + use setprob_module, only: tide_t, tide_etaprime, tide_m + + implicit none + + ! Input parameters + integer, intent(in) :: meqn,mbc,mx,my,maux + double precision, intent(in) :: xlower,ylower,dx,dy,t,dt + + ! Output + double precision, intent(inout) :: q(meqn,1-mbc:mx+mbc,1-mbc:my+mbc) + double precision, intent(inout) :: aux(maux,1-mbc:mx+mbc,1-mbc:my+mbc) + + ! Locals + integer :: i, j, nman, k + real(kind=8) :: h, hu, hv, gamma, dgamma, y, fdt, a(2,2), coeff + real(kind=8) :: xm, xc, xp, ym, yc, yp, dx_meters, dy_meters + real(kind=8) :: u, v, hu0, hv0 + real(kind=8) :: tau, wind_speed, theta, phi, psi, P_gradient(2), S(2) + real(kind=8) :: Ddt, sloc(2) + real(kind=8) :: eta,x,tperiod,deta_dt,ramp,alpha,ampl + real(kind=8) :: x1rs,x2rs,y1rs,y2rs,area,discharge_cfs,discharge,rs + + ! Algorithm parameters + ! Parameter controls when to zero out the momentum at a depth in the + ! friction source term + real(kind=8), parameter :: depth_tolerance = 1.0d-30 + + ! Physics + ! Nominal density of water + real(kind=8), parameter :: rho = 1025.d0 + + + ! Tidal forcing + ! This is applied as a source term based on desired eta'(t) offshore + + if (.true.) then + ! Sinusoidal forcing as a test: + tperiod = 12*3600. + ampl = 1.d0 + eta = ampl*sin(2.d0*pi*t/tperiod) + deta_dt = ampl*2.d0*pi/tperiod * cos(2.d0*pi*t/tperiod) + endif + + if (.false.) then + ! tidal forcing read in by setprob.f90: + do k=1,tide_m + if (t < tide_t(k)) then + exit + endif + enddo + + if ((k>1) .and. (k <= tide_m)) then + alpha = (t - tide_t(k-1)) / (tide_t(k) - tide_t(k-1)) + deta_dt = (1.d0-alpha)*tide_etaprime(k-1) + alpha*tide_etaprime(k) + else + deta_dt = 0.d0 + endif + endif + + do i=1-mbc,mx+mbc + x = xlower + (i-0.5d0)*dx + if (x < -124.3d0) then + do j=1-mbc,my+mbc + q(1,i,j) = q(1,i,j) + dt*deta_dt + enddo + else if (x < -124.19d0) then + ! ramp down to zero closer to shore: + ramp = 1.d0 - (x + 124.3d0)/(-124.19d0 + 124.3d0) + do j=1-mbc,my+mbc + q(1,i,j) = q(1,i,j) + dt*deta_dt*ramp + enddo + endif + enddo + + + ! River sources: + + if (.false.) then + + ! Johns River + x1rs = -123.951d0 + x2rs = -123.950d0 + y1rs = 46.8815d0 + y2rs = 46.8825d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 5000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + endif + + if (.false.) then + + ! Chehalis + + x1rs = -123.743d0 + x2rs = -123.740d0 + y1rs = 46.956d0 + y2rs = 46.958d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 20000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + + ! Wishkah River + x1rs = -123.811d0 + x2rs = -123.810d0 + y1rs = 46.990d0 + y2rs = 46.991d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 4000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + + + ! Hoquiam River + x1rs = -123.886d0 + x2rs = -123.885d0 + y1rs = 46.995d0 + y2rs = 46.996d0 + + area = (x2rs-x1rs)*(y2rs-y1rs)*(111d3)**2 * cos(y1rs*DEG2RAD) ! m**2 + discharge_cfs = 9000d0 ! cubic feet / second + discharge = discharge_cfs *(0.3048)**3 ! m**3 / second + rs = discharge/area ! m/sec needed in source area + + do i=1,mx + x = xlower + (i-0.5d0)*dx + do j=1,my + y = ylower + (j-0.5d0)*dy + if ((x >= x1rs) .and. (x <= x2rs) .and. & + (y >= y1rs) .and. (y <= y2rs)) then + q(1,i,j) = q(1,i,j) + dt*rs + endif + enddo + enddo + + endif + + ! Friction source term + if (friction_forcing) then + do j=1,my + do i=1,mx + ! Extract appropriate momentum + if (q(1,i,j) < depth_tolerance) then + q(2:3,i,j) = 0.d0 + else + ! Apply friction source term only if in shallower water + if (q(1,i,j) <= friction_depth) then + if (.not.variable_friction) then + do nman = num_manning, 1, -1 + if (aux(1,i,j) .lt. manning_break(nman)) then + coeff = manning_coefficient(nman) + endif + enddo + else + coeff = aux(friction_index,i,j) + endif + + ! Calculate source term + gamma = sqrt(q(2,i,j)**2 + q(3,i,j)**2) * g & + * coeff**2 / (q(1,i,j)**(7.d0/3.d0)) + dgamma = 1.d0 + dt * gamma + q(2, i, j) = q(2, i, j) / dgamma + q(3, i, j) = q(3, i, j) / dgamma + endif + endif + enddo + enddo + endif + ! End of friction source term + + ! Coriolis source term + ! TODO: May want to remove the internal calls to coriolis as this could + ! lead to slow downs. + if (coriolis_forcing) then + do j=1,my + y = ylower + (j - 0.5d0) * dy + fdt = coriolis(y) * dt ! Calculate f dependent on coordinate system + + ! Calculate matrix components + a(1,1) = 1.d0 - 0.5d0 * fdt**2 + fdt**4 / 24.d0 + a(1,2) = fdt - fdt**3 / 6.d0 + a(2,1) = -fdt + fdt**3 / 6.d0 + a(2,2) = a(1,1) + + do i=1,mx + q(2,i,j) = q(2, i, j) * a(1,1) + q(3, i, j) * a(1,2) + q(3,i,j) = q(2, i, j) * a(2,1) + q(3, i, j) * a(2,2) + enddo + enddo + endif + ! End of coriolis source term + + ! wind ----------------------------------------------------------- + if (wind_forcing) then + ! Need storm location and direction for sector based wind drag + sloc = storm_location(t) + theta = storm_direction(t) + do j=1,my + yc = ylower + (j - 0.5d0) * dy + do i=1,mx + xc = xlower + (i - 0.5d0) * dx + if (q(1,i,j) > dry_tolerance) then + psi = atan2(yc - sloc(2), xc - sloc(1)) + if (theta > psi) then + phi = (2.d0 * pi - theta + psi) * RAD2DEG + else + phi = (psi - theta) * RAD2DEG + endif + wind_speed = sqrt(aux(wind_index,i,j)**2 & + + aux(wind_index+1,i,j)**2) + tau = wind_drag(wind_speed, phi) * rho_air * wind_speed / rho + q(2,i,j) = q(2,i,j) + dt * tau * aux(wind_index,i,j) + q(3,i,j) = q(3,i,j) + dt * tau * aux(wind_index+1,i,j) + endif + enddo + enddo + endif + ! ---------------------------------------------------------------- + + ! Atmosphere Pressure -------------------------------------------- + ! Handled in Riemann solver + ! if (pressure_forcing) then + ! do j=1,my + ! ym = ylower + (j - 1.d0) * dy + ! yc = ylower + (j - 0.5d0) * dy + ! yp = ylower + j * dy + ! do i=1,mx + ! xm = xlower + (i - 1.d0) * dx + ! xc = xlower + (i - 0.5d0) * dx + ! xp = xlower + i * dx + + ! if (coordinate_system == 2) then + ! ! Convert distance in lat-long to meters + ! dx_meters = spherical_distance(xp,yc,xm,yc) + ! dy_meters = spherical_distance(xc,yp,xc,ym) + ! else + ! dx_meters = dx + ! dy_meters = dy + ! endif + + ! ! Extract depths + ! h = q(1,i,j) + + ! ! Calculate gradient of Pressure + ! P_gradient(1) = (aux(pressure_index,i+1,j) & + ! - aux(pressure_index,i-1,j)) / (2.d0 * dx_meters) + ! P_gradient(2) = (aux(pressure_index,i,j+1) & + ! - aux(pressure_index,i,j-1)) / (2.d0 * dy_meters) + + ! ! Modify momentum in each layer + ! if (h > dry_tolerance) then + ! q(2, i, j) = q(2, i, j) - dt * h * P_gradient(1) / rho + ! q(3, i, j) = q(3, i, j) - dt * h * P_gradient(2) / rho + ! end if + ! enddo + ! enddo + ! endif + +end subroutine src2 diff --git a/GraysHarborBC/topo/GraysHarborMTL.ipynb b/GraysHarborBC/topo/GraysHarborMTL.ipynb new file mode 100644 index 0000000..fb3c634 --- /dev/null +++ b/GraysHarborBC/topo/GraysHarborMTL.ipynb @@ -0,0 +1,170 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Grays Harbor Topography\n", + "\n", + "Modify topography from the Astoria 1/3\" DEM, which was referenced to MHW, by an approximate adjustment to Mean Tide Level (MTL). Around Grays Harbor, this adjustment is approximately linear in longitude, as verified separately using the [VDatum software](https://vdatum.noaa.gov/vdatumweb/). We use the datums at Westport and Aberdeen to estimate the linear function." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-py" + ] + }, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pylab import *" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from clawpack.geoclaw import topotools" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "server = 'https://www.ngdc.noaa.gov/thredds/dodsC/regional/'\n", + "url_astoria = server + 'astoria_13_mhw_2012.nc'\n", + "extent = [-124.2, -123.65, 46.8, 47.15]\n", + "GH_13sec_mhw = topotools.read_netcdf(url_astoria, extent=extent)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "figure(figsize=(13,7))\n", + "ax = axes()\n", + "GH_13sec_mhw.crop(coarsen=3).plot(axes=ax)\n", + "title('Astoria MHW');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Approximate conversion from MHW to MTL\n", + "\n", + "Datums at Westport and Aberdeen tide gauges:\n", + "\n", + "Datum | Westport | Aberdeen\n", + "------|----------|------------\n", + " MHW | 2.561 | 2.869\n", + " diff | 1.068 | 1.21\n", + " MTL | 1.493 | 1.659\n", + " diff | 1.068 | 1.21\n", + " MLW | 0.425 | 0.449\n", + " \n", + "Topo referenced to MHW should be increased by 1.068m at Westport and by 1.21m at Aberdeen to reference to MTL. We use a linear function based on longitude:\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Zmtl = GH_13sec_mhw.Z + (1.068 + (1.21-1.068)*(GH_13sec_mhw.X+124.105)/(-123.85+124.105))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GH_13sec_mtl = topotools.Topography()\n", + "GH_13sec_mtl.set_xyZ(GH_13sec_mhw.X, GH_13sec_mhw.Y, Zmtl)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "extent = [-123.9,-123.7,46.93,46.99]\n", + "figure(figsize=(13,13))\n", + "ax = subplot(211)\n", + "GH_13sec_mhw.crop(extent).plot(axes=ax,limits=(-10,10),\n", + " add_colorbar=True, cb_kwargs={'shrink':0.7})\n", + "title('Topo relative to MHW near Aberdeen')\n", + "ax = subplot(212)\n", + "GH_13sec_mtl.crop(extent).plot(axes=ax,limits=(-10,10),\n", + " add_colorbar=True, cb_kwargs={'shrink':0.7})\n", + "title('Topo relative to MTL near Aberdeen');" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create topo file:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GH_13sec_mtl.write('GH_13sec_mtl.asc', topo_type=3, \n", + " header_style='asc', Z_format='%10.3f')" + ] + }, + { + "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 +} diff --git a/GraysHarborBC/topo/README.md b/GraysHarborBC/topo/README.md new file mode 100644 index 0000000..917c35b --- /dev/null +++ b/GraysHarborBC/topo/README.md @@ -0,0 +1,11 @@ +# GraysHarbor/topo + +Use [GraysHarborMTL.ipynb](topo/GraysHarborMTL.html) to download 1/3 arcsecond +topography data that is referenced to MHW and adjust it so that it is +more closely aligned with Mean Tide Level (MTL). Since we are modeling +the tide, we want to use topography referenced to MTL if possible. + +Use `fetch_etopo.py` to download 1 arcminute topography of the coastal +region from the etopo1 database. + +Creat`topofiles` subdirectory and move the resulting .asc files into it. diff --git a/GraysHarborBC/topo/fetch_etopo.py b/GraysHarborBC/topo/fetch_etopo.py new file mode 100644 index 0000000..c3e0a93 --- /dev/null +++ b/GraysHarborBC/topo/fetch_etopo.py @@ -0,0 +1,53 @@ +""" +Simple script to download etopo1 topography/bathymetry data from + http://www.ngdc.noaa.gov/mgg/global/global.html + +The etopo1 data has 1-arcminute resolution, but you can request coarsening. +E.g. set resolution = 4./60. for 4-arcminute resolution. + +""" + +import os +from clawpack.geoclaw import topotools + +plot_topo = True + +# Set the limits of the domain and resolution: + +resolution = 1./60. # degrees + +# If user environment variable ETOPO is set to a valid path to a directory, +# downloaded data will be stored there. Allows sharing same data among +# various projects. + +try: + etopo_dir = os.environ['ETOPO'] + os.chdir(etopo_dir) # make sure it's a valid directory +except: + print("ETOPO not set or invalid directory, setting etopo_dir='.'") + etopo_dir = '.' + +extent = [-125,-123,46,48] +coarsen = 1 +topo = topotools.read_netcdf('etopo1', extent=extent, + coarsen=coarsen, verbose=True) + +name = 'etopo1_%s_%s_%s_%s_%smin' % tuple(extent + [coarsen]) +print('name = ',name) +fname = os.path.join(etopo_dir, name + '.asc') +topo.write(fname, topo_type=3, header_style='asc', + grid_registration='llcorner', Z_format='%.0f') + + + +if plot_topo: + # plot the topo and save as a png file... + import matplotlib.pyplot as plt + topo.plot() + plt.title('Topo file %s' % name) + fname = name + '.png' + fname = os.path.join(etopo_dir, fname) + plt.savefig(fname) + print('Created %s' % fname) + + diff --git a/NewYorkBC/README.md b/NewYorkBC/README.md new file mode 100644 index 0000000..a71ccee --- /dev/null +++ b/NewYorkBC/README.md @@ -0,0 +1,13 @@ +## NewYorkBC +Sets up basic New York tidal example using boundary conditions + +## `topography` directory + +Use `fetch_etopo.py` to download 1 arcminute topography of the coastal +region from the etopo1 database. + +## `sine` directory + +In this directory, the right boundary condition is implemented using +pure sine wave oscillation with amplitude 1. The notebook shows the simulation. + diff --git a/NewYorkBC/sine/Makefile b/NewYorkBC/sine/Makefile new file mode 100644 index 0000000..de2fa0a --- /dev/null +++ b/NewYorkBC/sine/Makefile @@ -0,0 +1,65 @@ +# Makefile for Clawpack code in this directory. +# This version only sets the local files and frequently changed +# options, and then includes the standard makefile pointed to by CLAWMAKE. +CLAWMAKE = $(CLAW)/clawutil/src/Makefile.common + +# See the above file for details and a list of make options, or type +# make .help +# at the unix prompt. + + +# Adjust these variables if desired: +# ---------------------------------- + +CLAW_PKG = geoclaw # Clawpack package to use +EXE = xgeoclaw # Executable to create +SETRUN_FILE = setrun.py # File containing function to make data +OUTDIR = _output # Directory for output +SETPLOT_FILE = setplot.py # File containing function to set plots +PLOTDIR = _plots # Directory for plots + +# Environment variable FC should be set to fortran compiler, e.g. gfortran + +# Compiler flags can be specified here or set as an environment variable +FFLAGS ?= + +# --------------------------------- +# package sources for this program: +# --------------------------------- + +GEOLIB = $(CLAW)/geoclaw/src/2d/shallow +include $(GEOLIB)/Makefile.geoclaw + +# --------------------------------------- +# package sources specifically to exclude +# (i.e. if a custom replacement source +# under a different name is provided) +# --------------------------------------- + +EXCLUDE_MODULES = \ + +EXCLUDE_SOURCES = \ + +# ---------------------------------------- +# List of custom sources for this program: +# ---------------------------------------- + +SOURCES = \ + bc2amr.f90 \ + $(CLAW)/riemann/src/rpn2_geoclaw.f \ + $(CLAW)/riemann/src/rpt2_geoclaw.f \ + $(CLAW)/riemann/src/geoclaw_riemann_utils.f \ + +#------------------------------------------------------------------- +# Include Makefile containing standard definitions and make options: +include $(CLAWMAKE) + +# Construct the topography data +.PHONY: topo all +topo: + python maketopo.py + +all: + $(MAKE) topo + $(MAKE) .plots + $(MAKE) .htmls diff --git a/NewYorkBC/sine/NYSine.ipynb b/NewYorkBC/sine/NYSine.ipynb new file mode 100644 index 0000000..70a072f --- /dev/null +++ b/NewYorkBC/sine/NYSine.ipynb @@ -0,0 +1,116 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# New York Tides with Sine Boundary Condition\n", + "\n", + "First run the GeoClaw code and then use this notebook to examine the gauge results.\n", + "\n", + "In this example, we set tidal signal to be $\\sin(2\\pi t / T)$ on the right boundary, a pure sine wave with amplitude 1 and period $T = 12\\times 3600 = $ 12 hours, as shown by the dashed line in the gauge time series plots below. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pylab import *" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from IPython.display import Image\n", + "import clawpack.pyclaw.gauges as gauges" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "figure(400, figsize=(13,6))\n", + "clf()\n", + "colors = ['k','b','r','g']\n", + "\n", + "outdir = '_output'\n", + "\n", + "for k,gaugeno in enumerate([0,1,5186,9050]):\n", + " gauge = gauges.GaugeSolution(gaugeno, outdir)\n", + " t = gauge.t / 3600. # convert to hours\n", + " q = gauge.q\n", + " eta = q[3,:]\n", + " plot(t, eta, colors[k], label='Gauge %s' % gaugeno)\n", + " \n", + " \n", + " # determine amplification and time shift:\n", + " m2 = int(floor(0.75*len(eta)))\n", + " eta2 = eta[m2:] # last part of eta signal\n", + " etamax2 = eta2.max()\n", + " etamin2 = eta2.min()\n", + " t2 = t[m2:]\n", + " jtmax = argmax(eta2) \n", + " tshift = t2[jtmax] *3600.\n", + " \n", + " print('At gauge %i, etamin2 = %.3f, etamax2 = %.3f at tshift = %.1f s' \\\n", + " % (gaugeno,etamin2,etamax2,tshift))\n", + " \n", + "#Sine Wave \n", + "tperiod = 12\n", + "eta = 1.*sin(2*pi*t/tperiod)\n", + "plot(t, eta, 'k--', label='Sine')\n", + "\n", + "legend(loc='upper right')\n", + "xlabel('hours')\n", + "ylabel('Surface relative to MTL (m)')\n", + "grid(True)\n", + "title('Sine wave BC and resulting GeoClaw gauge results');\n", + "\n", + "xticks(arange(0,t[-1]+0.1,12))\n", + "xlim(0,60)\n", + "\n", + "if 1:\n", + " fname = 'GaugeComparison.png'\n", + " savefig(fname, bbox_inches='tight')\n", + " print('Created %s' % fname)" + ] + } + ], + "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 +} diff --git a/NewYorkBC/sine/bc2amr.f90 b/NewYorkBC/sine/bc2amr.f90 new file mode 100644 index 0000000..f0bde62 --- /dev/null +++ b/NewYorkBC/sine/bc2amr.f90 @@ -0,0 +1,342 @@ +! :::::::::: bc2amr ::::::::::::::::::::::::::::::::::::::::::::::; +!> \callgraph +!! \callergraph +!! Take a grid patch with mesh widths **hx**,**hy**, of dimensions **nrow** by +!! **ncol**, and set the values of any piece of +!! of the patch which extends outside the physical domain +!! using the boundary conditions. +!! +!! +!! Specific to geoclaw: extrapolates aux(i,j,1) at boundaries +!! to constant. +!! +!! ### Standard boundary condition choices for amr2ez in clawpack +!! +!! At each boundary k = 1 (left), 2 (right), 3 (bottom), 4 (top): +!! +!! mthbc(k) = +!! * 0 for user-supplied BC's (must be inserted!) +!! * 1 for zero-order extrapolation +!! * 2 for periodic boundary conditions +!! * 3 for solid walls, assuming this can be implemented +!! by reflecting the data about the boundary and then +!! negating the 2'nd (for k=1,2) or 3'rd (for k=3,4) +!! component of q. +!! * 4 for sphere bcs (left half maps to right half of same side, and vice versa), as if domain folded in half +!! +!! The corners of the grid patch are at +!! (xlo_patch,ylo_patch) -- lower left corner +!! (xhi_patch,yhi_patch) -- upper right corner +!! +!! The physical domain itself is a rectangle bounded by +!! (xlower,ylower) -- lower left corner +!! (xupper,yupper) -- upper right corner +!! +! This figure below does not work with doxygen +! the picture is the following: +! ____________________________________________________ +! +! _____________________ (xupper,yupper) +! | | +! ____|____ (xhi_patch,yhi_patch) +! | | | | +! | | | | +! | | | | +! |___|____| | +! (xlo_patch,ylo_patch) | | +! | | +! |_____________________| +! (xlower,ylower) +! ____________________________________________________ +!! +!! +!> Any cells that lie outside the physical domain are ghost cells whose +!! values should be set in this routine. This is tested for by comparing +!! xlo_patch with xlower to see if values need to be set at the left +! as in the figure above, +! +!> and similarly at the other boundaries. +!! Patches are guaranteed to have at least 1 row of cells filled +!! with interior values so it is possible to extrapolate. +!! Fix [trimbd()](@ref trimbd) if you want more than 1 row pre-set. +!! +!! Make sure the order the boundaries are specified is correct +!! so that diagonal corner cells are also properly taken care of. +!! +!! Periodic boundaries are set before calling this routine, so if the +!! domain is periodic in one direction only you +!! can safely extrapolate in the other direction. +!! +!! Don't overwrite ghost cells in periodic directions! +!! +!! \param val data array for solution \f$q \f$ (cover the whole grid **msrc**) +!! \param aux data array for auxiliary variables +!! \param nrow number of cells in *i* direction on this grid +!! \param ncol number of cells in *j* direction on this grid +!! \param meqn number of equations for the system +!! \param naux number of auxiliary variables +!! \param hx spacing (mesh size) in *i* direction +!! \param hy spacing (mesh size) in *j* direction +!! \param level AMR level of this grid +!! \param time setting ghost cell values at time **time** +!! \param xlo_patch left bound of the input grid +!! \param xhi_patch right bound of the input grid +!! \param ylo_patch lower bound of the input grid +!! \param yhi_patch upper bound of the input grid +! ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::; + +subroutine bc2amr(val,aux,nrow,ncol,meqn,naux, hx, hy, level, time, & + xlo_patch, xhi_patch, ylo_patch, yhi_patch) + use amr_module + use geoclaw_module, only: pi + !use amr_module, only: mthbc, xlower, ylower, xupper, yupper + !use amr_module, only: xperdom,yperdom,spheredom + + implicit none + + ! Input/Output + integer, intent(in) :: nrow, ncol, meqn, naux, level + real(kind=8), intent(in) :: hx, hy, time + real(kind=8), intent(in) :: xlo_patch, xhi_patch + real(kind=8), intent(in) :: ylo_patch, yhi_patch + real(kind=8), intent(in out) :: val(meqn, nrow, ncol) + real(kind=8), intent(in out) :: aux(naux, nrow, ncol) + + ! Local storage + integer :: i, j, ibeg, jbeg, nxl, nxr, nyb, nyt, t_num + real(kind=8) :: hxmarg, hymarg, jump_h, h_r, h_m, u_r, delt + real, DIMENSION(1201) :: a + INTEGER :: col,max_cols + delt = possk(level) + jump_h = sin(2*pi*(time+delt)/(12*3600))- sin(2*pi*(time)/(12*3600)) + hxmarg = hx * .01d0 + hymarg = hy * .01d0 + + ! Use periodic boundary condition specialized code only, if only one + ! boundary is periodic we still proceed below + if (xperdom .and. (yperdom .or. spheredom)) then + return + end if + + ! Each check has an initial check to ensure that the boundary is a real + ! boundary condition and otherwise skips the code. Otherwise + !------------------------------------------------------- + ! Left boundary: + !------------------------------------------------------- + if (xlo_patch < xlower-hxmarg) then + ! number of grid cells from this patch lying outside physical domain: + nxl = int((xlower + hxmarg - xlo_patch) / hx) + + select case(mthbc(1)) + case(0) ! User defined boundary condition + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(3, i, j) = val(3, nxl + 1, j) + val(1, i, j) = val(1, nxl + 1, j) + jump_h + val(2, i, j) = val(2, nxl + 1, j) + + end do + end do + case(1) ! Zero-order extrapolation + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, nxl + 1, j) + val(:, i, j) = val(:, nxl + 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, ncol + do i=1, nxl + aux(:, i, j) = aux(:, 2 * nxl + 1 - i, j) + val(:, i, j) = val(:, 2 * nxl + 1 - i, j) + end do + end do + ! negate the normal velocity: + do j = 1, ncol + do i=1, nxl + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + end select + end if + + !------------------------------------------------------- + ! Right boundary: + !------------------------------------------------------- + if (xhi_patch > xupper+hxmarg) then + + ! number of grid cells lying outside physical domain: + nxr = int((xhi_patch - xupper + hxmarg) / hx) + ibeg = max(nrow - nxr + 1, 1) + + select case(mthbc(2)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(3, i, j) = val(3, ibeg - 1, j) + val(1, i, j) = val(1, ibeg - 1, j) + jump_h + val(2, i, j) = val(2, ibeg - 1, j) + end do + end do + + case(1) ! Zero-order extrapolation + do i = ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, ibeg - 1, j) + val(:, i, j) = val(:, ibeg - 1, j) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do i=ibeg, nrow + do j = 1, ncol + aux(:, i, j) = aux(:, 2 * ibeg - 1 - i, j) + val(:, i, j) = val(:, 2 * ibeg - 1 - i, j) + end do + end do + ! negate the normal velocity: + do i = ibeg, nrow + do j = 1, ncol + val(2, i, j) = -val(2, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Bottom boundary: + !------------------------------------------------------- + if (ylo_patch < ylower - hymarg) then + + ! number of grid cells lying outside physical domain: + nyb = int((ylower + hymarg - ylo_patch) / hy) + + select case(mthbc(3)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + do j = 1, nyb + do i = 1, nrow + aux(:, i, j) = aux(:, i, nyb + 1) + val(3, i, j) = val(3, i, nyb + 1) + val(1, i, j) = val(1, i, nyb + 1) + jump_h + val(2, i, j) = val(2, i, nyb + 1) + end do + end do + + case(1) ! Zero-order extrapolation + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, nyb + 1) + val(:,i,j) = val(:, i, nyb + 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = 1, nyb + do i = 1, nrow + aux(:,i,j) = aux(:, i, 2 * nyb + 1 - j) + val(:,i,j) = val(:, i, 2 * nyb + 1 - j) + end do + end do + ! negate the normal velocity: + do j = 1, nyb + do i = 1, nrow + val(3,i,j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + + !------------------------------------------------------- + ! Top boundary: + !------------------------------------------------------- + if (yhi_patch > yupper + hymarg) then + + ! number of grid cells lying outside physical domain: + nyt = int((yhi_patch - yupper + hymarg) / hy) + jbeg = max(ncol - nyt + 1, 1) + + select case(mthbc(4)) + case(0) ! User defined boundary condition + ! Replace this code with a user defined boundary condition + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(3, i, j) = val(3, i, jbeg - 1) + val(1, i, j) = val(1, i, jbeg - 1) + jump_h + val(2, i, j) = val(2, i, jbeg - 1) + end do + end do + + case(1) ! Zero-order extrapolation + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, jbeg - 1) + val(:, i, j) = val(:, i, jbeg - 1) + end do + end do + + case(2) ! Periodic boundary condition + continue + + case(3) ! Wall boundary conditions + do j = jbeg, ncol + do i = 1, nrow + aux(:, i, j) = aux(:, i, 2 * jbeg - 1 - j) + val(:, i, j) = val(:, i, 2 * jbeg - 1 - j) + end do + end do + ! negate the normal velocity: + do j = jbeg, ncol + do i = 1, nrow + val(3, i, j) = -val(3, i, j) + end do + end do + + case(4) ! Spherical domain + continue + + case default + print *, "Invalid boundary condition requested." + stop + + end select + end if + +end subroutine bc2amr diff --git a/NewYorkBC/sine/setplot.py b/NewYorkBC/sine/setplot.py new file mode 100644 index 0000000..5dfdb8b --- /dev/null +++ b/NewYorkBC/sine/setplot.py @@ -0,0 +1,528 @@ + +""" +Set up the plot figures, axes, and items to be done for each frame. + +This module is imported by the plotting routines and then the +function setplot is called to set the plot parameters. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d + +from clawpack.geoclaw import topotools +from clawpack.visclaw import gaugetools, plottools, colormaps + +cmin_water = -1. +cmax_water = 1. +cmin_speed = 0. +cmax_speed = 2. + +def speed(current_data): + """ + Return a masked array containing the speed only in wet cells. + Surface is eta = h+topo, assumed to be output as 4th column of fort.q + files. + """ + from clawpack.visclaw.geoplot import drytol_default + drytol = getattr(current_data.user, 'drytol', drytol_default) + q = current_data.q + h = q[0,:,:] + hu = q[1,:,:] + hv = q[2,:,:] + eta = q[3,:,:] + + hwet = np.ma.masked_where(h <= drytol, h) + speed = np.sqrt(hu**2 + hv**2) / hwet # will be masked where dry + + try: + # Use mask covering coarse regions if it's set: + m = current_data.mask_coarse + speed = numpy.ma.masked_where(m, speed) + except: + pass + + return speed + +if 1: + # boundary: + sea_level = 0. + tperiod = 12 # hours + eta_tide = lambda t: 1.*np.sin(2*np.pi*t/tperiod) + t_tide = np.linspace(0,48,1000) + + + +def plot_tide(current_data): + from pylab import linspace,plot,grid,title,cos,pi,xlabel,ylabel,xticks + from pylab import xlim,ylim + tframe = current_data.t / 3600. + tt = linspace(t_tide[0],t_tide[-1],1000) + plot(tt, eta_tide(tt), 'k') + plot([tframe],[eta_tide(tframe)], 'ro') + grid(True) + title('Tide stage in ocean') + xticks(np.arange(0,5.1*24,24)) + xlabel('hours') + ylabel('meters (0=MTL)') + xlim(0,48) + ylim(-1.2,1.2) + +def plot_src_boxes(): + # Showing Rivers, GraysHabor Example + if 0: + # Johns River + x1rs = -123.951 + x2rs = -123.950 + y1rs = 46.8815 + y2rs = 46.8825 + + if 0: + # Chehalis River: + x1rs = -123.743 + x2rs = -123.740 + y1rs = 46.956 + y2rs = 46.958 + kwargs = {'color':'yellow', 'linewidth':1} + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + # Wishkah River + x1rs = -123.811 + x2rs = -123.810 + y1rs = 46.990 + y2rs = 46.991 + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + # Hoquiam River + x1rs = -123.886 + x2rs = -123.885 + y1rs = 46.995 + y2rs = 46.996 + plottools.plotbox([x1rs,x2rs,y1rs,y2rs],kwargs) + + +#-------------------------- +def setplot(plotdata=None): +#-------------------------- + + """ + Specify what is to be plotted at each frame. + Input: plotdata, an instance of pyclaw.plotters.data.ClawPlotData. + Output: a modified version of plotdata. + + """ + + + from clawpack.visclaw import colormaps, geoplot + from numpy import linspace + + if plotdata is None: + from clawpack.visclaw.data import ClawPlotData + plotdata = ClawPlotData() + + + plotdata.clearfigures() # clear any old figures,axes,items data + plotdata.format = 'binary' # 'ascii' or 'binary' to match setrun.py + + + # To plot gauge locations on pcolor or contour plot, use this as + # an afteraxis function: + + def addgauges(current_data): + from clawpack.visclaw import gaugetools + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=[0,1,5186,9050], format_string='ko', add_labels=True) + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Domain', figno=50) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + plotaxes.xlimits = [-74.4,-71.1] + plotaxes.ylimits = [40.1, 41.3] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + + + #----------------------------------------- + # Figure for surface + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface', figno=0) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.xlimits = [-124.23,-123.65] + plotaxes.xlimits = [-74.4,-71.1] + plotaxes.ylimits = [40.1, 41.3] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Surface at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + + + #----------------------------------------- + # Figure for speed + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Speed', figno=1) + plotfigure.kwargs = {'figsize':(13,8)} + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.6])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + #plotaxes.xlimits = [-124.23,-123.65] + plotaxes.xlimits = [-74.4,-71.1] + plotaxes.ylimits = [40.1, 41.3] + + def fixup(current_data): + import pylab + addgauges(current_data) + t = current_data.t + t = t / 3600. # hours + pylab.title('Speed at %4.2f hours' % t, fontsize=15) + pylab.xticks(fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = speed + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = cmin_speed + plotitem.pcolor_cmax = cmax_speed + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.6 + plotitem.colorbar_kwargs = {'extend':'max'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # add contour lines of bathy if desired: + plotitem = plotaxes.new_plotitem(plot_type='2d_contour') + plotitem.show = False + plotitem.plot_var = geoplot.topo + plotitem.contour_levels = linspace(-3000,-3000,1) + plotitem.amr_contour_colors = ['y'] # color on each level + plotitem.kwargs = {'linestyles':'solid','linewidths':2} + plotitem.amr_contour_show = [1,0,0] + plotitem.celledges_show = 0 + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Aberdeen', figno=10) + plotfigure.kwargs = {'figsize':(11,8)} + #plotfigure.show = False + xlimits = [-74.4,-71.1] + ylimits = [40.1, 41.3] + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.5])' + plotaxes.title = 'Surface' + plotaxes.scaled = True + plotaxes.xlimits = xlimits + plotaxes.ylimits = ylimits + + + def fixup(current_data): + import pylab + #addgauges(current_data) + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=range(54,58), format_string='ko', add_labels=True) + t = current_data.t + t = t / 3600. # hours + pylab.title('Aberdeen at t = %.2f hours' % t, fontsize=15) + pylab.xticks(rotation=20,fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + pylab.ticklabel_format(useOffset=False) + plot_src_boxes() + + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + #plotitem.plot_var = geoplot.surface + plotitem.plot_var = geoplot.surface_or_depth + plotitem.pcolor_cmap = geoplot.tsunami_colormap + plotitem.pcolor_cmin = cmin_water + plotitem.pcolor_cmax = cmax_water + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.7 + plotitem.colorbar_kwargs = {'extend':'both'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figure for zoom + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Aberdeen speed', figno=11) + plotfigure.kwargs = {'figsize':(11,8)} + #plotfigure.show = False + + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes('pcolor') + plotaxes.axescmd = 'axes([.1,.1,.85,.5])' + plotaxes.title = 'Speed' + plotaxes.scaled = True + plotaxes.xlimits = xlimits + plotaxes.ylimits = ylimits + + def fixup(current_data): + import pylab + #addgauges(current_data) + gaugetools.plot_gauge_locations(current_data.plotdata, \ + gaugenos=range(54,58), format_string='ko', add_labels=True) + t = current_data.t + t = t / 3600. # hours + pylab.title('Aberdeen speed at t = %.2f hours' % t, fontsize=15) + pylab.xticks(rotation=20,fontsize=10) + pylab.yticks(fontsize=10) + pylab.gca().set_aspect(1/pylab.cos(47*pylab.pi/180)) + pylab.ticklabel_format(useOffset=False) + plot_src_boxes() + plotaxes.afteraxes = fixup + + # Water + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = speed + plotitem.pcolor_cmap = colormaps.blue_yellow_red + plotitem.pcolor_cmin = cmin_speed + plotitem.pcolor_cmax = cmax_speed + plotitem.add_colorbar = True + plotitem.colorbar_shrink = 0.7 + plotitem.colorbar_kwargs = {'extend':'max'} + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + # Land + plotitem = plotaxes.new_plotitem(plot_type='2d_pcolor') + plotitem.plot_var = geoplot.land + plotitem.pcolor_cmap = geoplot.land_colors + plotitem.pcolor_cmin = 0.0 + plotitem.pcolor_cmax = 100.0 + plotitem.add_colorbar = False + plotitem.amr_celledges_show = [0,0,0] + plotitem.patchedges_show = 0 + + plotaxes = plotfigure.new_plotaxes('1d_plot') + plotaxes.axescmd = 'axes([.1,.75,.7,.15])' + plotaxes.afteraxes = plot_tide + + #----------------------------------------- + # Figures for gauges + #----------------------------------------- + plotfigure = plotdata.new_plotfigure(name='Surface at gauges', figno=300, \ + type='each_gauge') + plotfigure.clf_each_gauge = True + + # Set up for axes in this figure: + plotaxes = plotfigure.new_plotaxes() + plotaxes.xlimits = 'auto' + plotaxes.ylimits = 'auto' + plotaxes.title = 'Surface' + + # Plot surface as blue curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.plot_var = 3 + plotitem.plotstyle = 'b-' + + # Plot topo as green curve: + plotitem = plotaxes.new_plotitem(plot_type='1d_plot') + plotitem.show = False + + def gaugetopo(current_data): + q = current_data.q + h = q[0,:] + eta = q[3,:] + topo = eta - h + return topo + + plotitem.plot_var = gaugetopo + plotitem.plotstyle = 'g-' + + + + #----------------------------------------- + # Plots of timing (CPU and wall time): + + def make_timing_plots(plotdata): + from clawpack.visclaw import plot_timing_stats + import os,sys + try: + timing_plotdir = plotdata.plotdir + '/_timing_figures' + os.system('mkdir -p %s' % timing_plotdir) + # adjust units for plots based on problem: + units = {'comptime':'seconds', 'simtime':'hours', + 'cell':'millions'} + plot_timing_stats.make_plots(outdir=plotdata.outdir, + make_pngs=True, + plotdir=timing_plotdir, + units=units) + except: + print('*** Error making timing plots') + + otherfigure = plotdata.new_otherfigure(name='timing plots', + fname='_timing_figures/timing.html') + otherfigure.makefig = make_timing_plots + + + #----------------------------------------- + + # Parameters used only when creating html and/or latex hardcopy + # e.g., via pyclaw.plotters.frametools.printframes: + + plotdata.printfigs = True # print figures + plotdata.print_format = 'png' # file format + plotdata.print_framenos = 'all' # list of frames to print + plotdata.print_gaugenos = 'all' # list of gauges to print + plotdata.print_fignos = 'all' # list of figures to print + plotdata.html = True # create html files of plots? + plotdata.html_homelink = '../README.html' # pointer for top of index + plotdata.latex = True # create latex file of plots? + plotdata.latex_figsperline = 2 # layout of plots + plotdata.latex_framesperline = 1 # layout of plots + plotdata.latex_makepdf = False # also run pdflatex? + plotdata.parallel = True # make multiple frame png's at once + + return plotdata diff --git a/NewYorkBC/sine/setrun.py b/NewYorkBC/sine/setrun.py new file mode 100644 index 0000000..7866250 --- /dev/null +++ b/NewYorkBC/sine/setrun.py @@ -0,0 +1,451 @@ +""" +Module to set up run time parameters for Clawpack. + +The values set in the function setrun are then written out to data files +that will be read in by the Fortran code. + +""" + +from __future__ import absolute_import +from __future__ import print_function +import os, sys +import numpy as np +from clawpack.amrclaw.data import FlagRegion + + +try: + CLAW = os.environ['CLAW'] +except: + raise Exception("*** Must first set CLAW enviornment variable") + + + +#------------------------------ +def setrun(claw_pkg='geoclaw'): +#------------------------------ + + """ + Define the parameters used for running Clawpack. + + INPUT: + claw_pkg expected to be "geoclaw" for this setrun. + + OUTPUT: + rundata - object of class ClawRunData + + """ + + from clawpack.clawutil import data + + assert claw_pkg.lower() == 'geoclaw', "Expected claw_pkg = 'geoclaw'" + + num_dim = 2 + rundata = data.ClawRunData(claw_pkg, num_dim) + + + #------------------------------------------------------------------ + # Problem-specific parameters to be written to setprob.data: + #------------------------------------------------------------------ + + #probdata = rundata.new_UserData(name='probdata',fname='setprob.data') + + #------------------------------------------------------------------ + # GeoClaw specific parameters: + #------------------------------------------------------------------ + rundata = setgeo(rundata) + + #------------------------------------------------------------------ + # Standard Clawpack parameters to be written to claw.data: + # (or to amr2ez.data for AMR) + #------------------------------------------------------------------ + clawdata = rundata.clawdata # initialized when rundata instantiated + + + # Set single grid parameters first. + # See below for AMR parameters. + + + # --------------- + # Spatial domain: + # --------------- + + # Number of space dimensions: + clawdata.num_dim = num_dim + + # Lower and upper edge of computational domain: + # x values should be integer multipes of 1/9" + # y values should be integer multipes of 1/9" + # Note: always satisfied if limits are multiples of 0.01 degree + + clawdata.lower[0] = -74.4 + clawdata.upper[0] = -71.1 + clawdata.lower[1] = 40.1 + clawdata.upper[1] = 41.3 + + clawdata.num_cells[0] = 40 + clawdata.num_cells[1] = 20 + + + # --------------- + # Size of system: + # --------------- + + # Number of equations in the system: + clawdata.num_eqn = 3 + + # Number of auxiliary variables in the aux array (initialized in setaux) + clawdata.num_aux = 3 + + # Index of aux array corresponding to capacity function, if there is one: + clawdata.capa_index = 2 + + + + # ------------- + # Initial time: + # ------------- + + clawdata.t0 = 0.0 + + + # Restart from checkpoint file of a previous run? + # If restarting, t0 above should be from original run, and the + # restart_file 'fort.chkNNNNN' specified below should be in + # the OUTDIR indicated in Makefile. + + clawdata.restart = False + clawdata.restart_file = '' + + # ------------- + # Output times: + #-------------- + + # Specify at what times the results should be written to fort.q files. + # Note that the time integration stops after the final output time. + # The solution at initial time t0 is always written in addition. + + clawdata.output_style = 1 + + if clawdata.output_style==1: + # Output nout frames at equally spaced times up to tfinal: + clawdata.num_output_times = 2*24 + clawdata.tfinal = 2*24*3600. + clawdata.output_t0 = True # output at initial (or restart) time? + + elif clawdata.output_style == 2: + # Specify a list of output times. + clawdata.output_times = [] + + elif clawdata.output_style == 3: + # Output every iout timesteps with a total of ntot time steps: + clawdata.output_step_interval = 1 + clawdata.total_steps = 1 + clawdata.output_t0 = True + + + clawdata.output_format = 'binary' + + clawdata.output_q_components = 'all' # need all + clawdata.output_aux_components = 'none' # eta=h+B is in q + clawdata.output_aux_onlyonce = False # output aux arrays each frame + + + + # --------------------------------------------------- + # Verbosity of messages to screen during integration: + # --------------------------------------------------- + + # The current t, dt, and cfl will be printed every time step + # at AMR levels <= verbosity. Set verbosity = 0 for no printing. + # (E.g. verbosity == 2 means print only on levels 1 and 2.) + clawdata.verbosity = 1 + + + + # -------------- + # Time stepping: + # -------------- + + # if dt_variable==1: variable time steps used based on cfl_desired, + # if dt_variable==0: fixed time steps dt = dt_initial will always be used. + clawdata.dt_variable = True + + # Initial time step for variable dt. + # If dt_variable==0 then dt=dt_initial for all steps: + clawdata.dt_initial = 0.2 + + # Max time step to be allowed if variable dt used: + clawdata.dt_max = 1e+99 + + # Desired Courant number if variable dt used, and max to allow without + # retaking step with a smaller dt: + clawdata.cfl_desired = 0.8 + clawdata.cfl_max = 1.0 + + # Maximum number of time steps to allow between output times: + clawdata.steps_max = 50000 + + + + + # ------------------ + # Method to be used: + # ------------------ + + # Order of accuracy: 1 => Godunov, 2 => Lax-Wendroff plus limiters + clawdata.order = 2 + + # Use dimensional splitting? (not yet available for AMR) + clawdata.dimensional_split = 'unsplit' + + # For unsplit method, transverse_waves can be + # 0 or 'none' ==> donor cell (only normal solver used) + # 1 or 'increment' ==> corner transport of waves + # 2 or 'all' ==> corner transport of 2nd order corrections too + clawdata.transverse_waves = 2 + + # Number of waves in the Riemann solution: + clawdata.num_waves = 3 + + # List of limiters to use for each wave family: + # Required: len(limiter) == num_waves + # Some options: + # 0 or 'none' ==> no limiter (Lax-Wendroff) + # 1 or 'minmod' ==> minmod + # 2 or 'superbee' ==> superbee + # 3 or 'mc' ==> MC limiter + # 4 or 'vanleer' ==> van Leer + clawdata.limiter = ['mc', 'mc', 'mc'] + + clawdata.use_fwaves = True # True ==> use f-wave version of algorithms + + # Source terms splitting: + # src_split == 0 or 'none' ==> no source term (src routine never called) + # src_split == 1 or 'godunov' ==> Godunov (1st order) splitting used, + # src_split == 2 or 'strang' ==> Strang (2nd order) splitting used, not recommended. + clawdata.source_split = 'godunov' + + + # -------------------- + # Boundary conditions: + # -------------------- + + # Number of ghost cells (usually 2) + clawdata.num_ghost = 2 + + # Choice of BCs at xlower and xupper: + # 0 => user specified (must modify bcN.f to use this option) + # 1 => extrapolation (non-reflecting outflow) + # 2 => periodic (must specify this at both boundaries) + # 3 => solid wall for systems where q(2) is normal velocity + + clawdata.bc_lower[0] = 1 + clawdata.bc_upper[0] = 0 + + clawdata.bc_lower[1] = 1 + clawdata.bc_upper[1] = 1 + + + + # -------------- + # Checkpointing: + # -------------- + + # Specify when checkpoint files should be created that can be + # used to restart a computation. + + # negative checkpoint_style means alternate between aaaaa and bbbbb files + # so that at most 2 checkpoint files exist at any time, useful when + # doing frequent checkpoints of large problems. + + clawdata.checkpt_style = 1 + + if clawdata.checkpt_style == 0: + # Do not checkpoint at all + pass + + elif clawdata.checkpt_style == 1: + # Checkpoint only at tfinal. + pass + + elif abs(clawdata.checkpt_style) == 2: + # Specify a list of checkpoint times. + clawdata.checkpt_times = 3600.*np.arange(1,16,1) + + elif abs(clawdata.checkpt_style) == 3: + # Checkpoint every checkpt_interval timesteps (on Level 1) + # and at the final time. + clawdata.checkpt_interval = 5 + + + # --------------- + # AMR parameters: + # --------------- + amrdata = rundata.amrdata + + # max number of refinement levels: + amrdata.amr_levels_max = 3 + + # List of refinement ratios at each level (length at least mxnest-1) + # 1', 12", 4", 1", 1/3" + amrdata.refinement_ratios_x = [5,3,4,3] + amrdata.refinement_ratios_y = [5,3,4,3] + amrdata.refinement_ratios_t = [5,3,4,3] + + + + # Specify type of each aux variable in amrdata.auxtype. + # This must be a list of length maux, each element of which is one of: + # 'center', 'capacity', 'xleft', or 'yleft' (see documentation). + + amrdata.aux_type = ['center','capacity','yleft'] + + + # Flag using refinement routine flag2refine rather than richardson error + amrdata.flag_richardson = False # use Richardson? + amrdata.flag2refine = True + + # steps to take on each level L between regriddings of level L+1: + amrdata.regrid_interval = 3 + + # width of buffer zone around flagged points: + # (typically the same as regrid_interval so waves don't escape): + amrdata.regrid_buffer_width = 2 + + # clustering alg. cutoff for (# flagged pts) / (total # of cells refined) + # (closer to 1.0 => more small grids may be needed to cover flagged cells) + amrdata.clustering_cutoff = 0.700000 + + # print info about each regridding up to this level: + amrdata.verbosity_regrid = 1 + + # ----- For developers ----- + # Toggle debugging print statements: + amrdata.dprint = False # print domain flags + amrdata.eprint = False # print err est flags + amrdata.edebug = False # even more err est flags + amrdata.gprint = False # grid bisection/clustering + amrdata.nprint = False # proper nesting output + amrdata.pprint = False # proj. of tagged points + amrdata.rprint = False # print regridding summary + amrdata.sprint = False # space/memory output + amrdata.tprint = False # time step reporting each level + amrdata.uprint = False # update/upbnd reporting + + # More AMR parameters can be set -- see the defaults in pyclaw/data.py + + # --------------- + # Regions: + # --------------- + #rundata.regiondata.regions = [] + # to specify regions of refinement append lines of the form + # [minlevel,maxlevel,t1,t2,x1,x2,y1,y2] + + # --------------- + # Gauges: + # --------------- + # for gauges append lines of the form [gaugeno, x, y, t1, t2] + gauges = rundata.gaugedata.gauges + rundata.gaugedata.min_time_increment = 1. # seconds between gauge output + rundata.gaugedata.display_format = "f15.7" # need enough digits for u,v + + gauges.append([0,-71.1,40.6,0,1e9]) + gauges.append([1,-72.5,40.4,0,1e9]) + + gauges.append([560,-71.96,41.04,0,1e9]) + gauges.append([9050,-74.06,40.61,0,1e9]) + gauges.append([5186,-73.26,40.62,0,1e9]) + + return rundata + # end of function setrun + # ---------------------- + + +#------------------- +def setgeo(rundata): +#------------------- + """ + Set GeoClaw specific runtime parameters. + For documentation see .... + """ + + try: + geo_data = rundata.geo_data + except: + print("*** Error, this rundata has no geo_data attribute") + raise AttributeError("Missing geo_data attribute") + + # == Physics == + geo_data.gravity = 9.81 + geo_data.coordinate_system = 2 + geo_data.earth_radius = 6367.5e3 + + # == Forcing Options + geo_data.coriolis_forcing = False + + # == Algorithm and Initial Conditions == + geo_data.sea_level = 0. + geo_data.dry_tolerance = 1.e-3 + geo_data.friction_forcing = True + geo_data.manning_coefficient =.025 + geo_data.friction_depth = 1e6 + + # Refinement settings + refinement_data = rundata.refinement_data + refinement_data.variable_dt_refinement_ratios = True + refinement_data.wave_tolerance = 0.05 + refinement_data.speed_tolerance = 6*[0.05] + + # == settopo.data values == + topo_data = rundata.topo_data + # for topography, append lines of the form + # [topotype, fname] + topofiles = topo_data.topofiles + + topodir = '../topography/' + topofiles.append([3, topodir + 'etopo1_-74.5_-71_40_41.4_1min.asc']) + + + + # == setdtopo.data values == + dtopo_data = rundata.dtopo_data + # for moving topography, append lines of the form : (<= 1 allowed for now!) + # [topotype, minlevel,maxlevel,fname] + dtopo_data.dtopofiles = [] + dtopo_data.dt_max_dtopo = 1.0 + + + # == setqinit.data values == + + + rundata.qinit_data.qinit_type = 0 + rundata.qinit_data.qinitfiles = [] + # for qinit perturbations, append lines of the form: (<= 1 allowed for now!) + # [minlev, maxlev, fname] + + # variable_eta_init newly added to QinitData: + rundata.qinit_data.variable_eta_init = False + + + # == fgmax.data values == + # MODIFIED FORMAT: + # Now specify a list of objects of class fgmax_tools.FGmaxGrid + # specifying any fgmax grids. + + rundata.fgmax_data.fgmax_grids = [] + rundata.fgmax_data.num_fgmax_val = 5 # Save depth, speed, hs, hss, hmin + + return rundata + # end of function setgeo + # ----------------------: + + + +if __name__ == '__main__': + # Set up run-time parameters and write all data files. + import sys + from clawpack.geoclaw import kmltools + rundata = setrun(*sys.argv[1:]) + rundata.write() + + kmltools.make_input_data_kmls(rundata) diff --git a/NewYorkBC/topography/fetch_etopo.py b/NewYorkBC/topography/fetch_etopo.py new file mode 100644 index 0000000..3913212 --- /dev/null +++ b/NewYorkBC/topography/fetch_etopo.py @@ -0,0 +1,51 @@ +""" +Simple script to download etopo1 topography/bathymetry data from + http://www.ngdc.noaa.gov/mgg/global/global.html + +The etopo1 data has 1-arcminute resolution, but you can request coarsening. +E.g. set resolution = 4./60. for 4-arcminute resolution. + +""" + +import os +from clawpack.geoclaw import topotools + +plot_topo = True + +# Set the limits of the domain and resolution: + +resolution = 1./60. # degrees + +# If user environment variable ETOPO is set to a valid path to a directory, +# downloaded data will be stored there. Allows sharing same data among +# various projects. + +try: + etopo_dir = os.environ['ETOPO'] + os.chdir(etopo_dir) # make sure it's a valid directory +except: + print("ETOPO not set or invalid directory, setting etopo_dir='.'") + etopo_dir = '.' + +extent = [-74.5,-71,40,41.4] +coarsen = 1 +topo = topotools.read_netcdf('etopo1', extent=extent, + coarsen=coarsen, verbose=True) + +name = 'etopo1_%s_%s_%s_%s_%smin' % tuple(extent + [coarsen]) +print('name = ',name) +fname = os.path.join(etopo_dir, name + '.asc') +topo.write(fname, topo_type=3, header_style='asc', + grid_registration='llcorner', Z_format='%.0f') + + + +if plot_topo: + # plot the topo and save as a png file... + import matplotlib.pyplot as plt + topo.plot() + plt.title('Topo file %s' % name) + fname = name + '.png' + fname = os.path.join(etopo_dir, fname) + plt.savefig(fname) + print('Created %s' % fname)