simulation_UnSteadyCollocated_dx.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dad54267",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import numba"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b62ed7a",
   "metadata": {},
   "source": [
    "# Unsteady flow : Collocated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d171d6b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "@numba.jit(nopython=True, parallel=False)\n",
    "def UnSteadyflowCollocated(A, Q, Adown, Qup, zb, B, dx, dt, g, manning):\n",
    "    imax = len(A)\n",
    "    Anew, Qnew = np.zeros(imax), np.zeros(imax)\n",
    "    \n",
    "# continuous equation\n",
    "    for i in numba.prange(1, imax-1) : Anew[i] = A[i] - dt * ( Q[i] - Q[i-1] ) / dx\n",
    "        \n",
    "    Anew[0], Anew[-1] = Anew[1], Adown\n",
    "    \n",
    "# moumentum equation\n",
    "    for i in numba.prange(1,imax-1): \n",
    "        ip, ic, im = (i+1, i, i-1) \n",
    "        Cr1 = 0.5*( Q[ic]/A[ic] + Q[ip]/A[ip] )*dt/dx\n",
    "        Cr2 = 0.5*( Q[ic]/A[ic] + Q[im]/A[im] )*dt/dx\n",
    "        dHdx1 = ( Anew[ip]/B[ip] + zb[ip] - Anew[ic]/B[ic] - zb[ic] ) / dx\n",
    "        dHdx2 = ( Anew[ic]/B[ic] + zb[ic] - Anew[im]/B[im] - zb[im] ) / dx\n",
    "        dHdx = (1.0 - Cr1) * dHdx1 + Cr2 * dHdx2\n",
    "        \n",
    "        Qnew[ic] = Q[ic] - dt * ( Q[ic]**2/A[ic] - Q[im]**2/A[im] ) / dx  \\\n",
    "                   - dt * g * Anew[ic] * dHdx \\\n",
    "                   - dt * g * A[ic] * manning**2 * Q[ic]**2 / B[ic]**2 / ( A[ic]/B[ic] )**(10.0/3.0)\n",
    "                \n",
    "#     Qnew[0], Qnew[-1] = Qup, Qnew[-2]\n",
    "    Qnew[0] = Qup\n",
    "    \n",
    "# check downstream boundary Q \n",
    "    i = -1\n",
    "    ic, im = (i, i-1) \n",
    "    dHdx = ( Anew[ic]/B[ic] + zb[ic] - Anew[im]/B[im] - zb[im] ) / dx\n",
    "    \n",
    "    Qnew[ic] = Q[ic] - dt * ( Q[ic]**2/A[ic] - Q[im]**2/A[im] ) / dx  \\\n",
    "               - dt * g * Anew[ic] * dHdx \\\n",
    "               - dt * g * A[ic] * manning**2 * Q[ic]**2 / B[ic]**2 / ( A[ic]/B[ic] )**(10.0/3.0)\n",
    "            \n",
    "    return Anew, Qnew, np.abs(Anew - A).max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b089c07f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# df = pd.read_csv('zb2.csv',index_col=0)\n",
    "df = pd.read_csv('https://raw.githubusercontent.com/computational-sediment-hyd/1DSWEwithSmallerDx/main/zb2.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c6175238",
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 10.0\n",
    "n = 0.03\n",
    "ib = 1/400\n",
    "g = 9.8\n",
    "\n",
    "h0 = (q**2*n**2/ib)**0.3 #等流水深\n",
    "hc = (q**2/g)**(1/3) # 限界水深"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f20db158",
   "metadata": {},
   "source": [
    "## dx=0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bf96d81",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "dx = 0.5\n",
    "Ls1 = np.arange(0,10000.1,dx)\n",
    "Zip1 = np.interp(Ls1, df.L, df.Z)\n",
    "\n",
    "imax = len(Ls1)\n",
    "dt = 0.025\n",
    "totalTime = 3.0*3600.0\n",
    "manning = n\n",
    "\n",
    "# Initial & Boundary condition\n",
    "zb = Zip1.copy()\n",
    "zb = zb[::-1]\n",
    "\n",
    "B = np.full(imax, 1.0, dtype=float)\n",
    "A = h0*B\n",
    "Q = q*B\n",
    "\n",
    "Qup = Q[0]\n",
    "Adown = A[-1] \n",
    "\n",
    "tout = 600\n",
    "for it in range(int(totalTime/dt)):\n",
    "    A, Q, err = UnSteadyflowCollocated(A, Q, Adown, Qup, zb, B, dx, dt, g, manning)\n",
    "    if (dt*it)%tout==0: print((dt*it), err)\n",
    "        \n",
    "hc1 = A[::-1]/B[::-1] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b41cc71",
   "metadata": {},
   "source": [
    "## dx=10.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b29c0e41",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "dx = 10.0\n",
    "Ls2 = np.arange(0,10000.1,dx)\n",
    "Zip2 = np.interp(Ls2, df.L, df.Z)\n",
    "\n",
    "imax = len(Ls2)\n",
    "dt = 0.025\n",
    "totalTime = 3.0*3600.0\n",
    "manning = n\n",
    "\n",
    "# Initial & Boundary condition\n",
    "zb = Zip2.copy()\n",
    "zb = zb[::-1]\n",
    "\n",
    "B = np.full(imax, 1.0, dtype=float)\n",
    "A = h0*B\n",
    "Q = q*B\n",
    "\n",
    "Qup = Q[0]\n",
    "Adown = A[-1] \n",
    "\n",
    "tout = 600\n",
    "for it in range(int(totalTime/dt)):\n",
    "    A, Q, err = UnSteadyflowCollocated(A, Q, Adown, Qup, zb, B, dx, dt, g, manning)\n",
    "    if (dt*it)%tout==0: print((dt*it), err)\n",
    "        \n",
    "hc2 = A[::-1]/B[::-1] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2dc92f2",
   "metadata": {},
   "source": [
    "## dx=100.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a2a646ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "dx = 100.0\n",
    "Ls3 = np.arange(0,10000.1,dx)\n",
    "Zip3 = np.interp(Ls3, df.L, df.Z)\n",
    "\n",
    "imax = len(Ls3)\n",
    "dt = 0.025\n",
    "totalTime = 3.0*3600.0\n",
    "manning = n\n",
    "\n",
    "# Initial & Boundary condition\n",
    "zb = Zip3.copy()\n",
    "zb = zb[::-1]\n",
    "\n",
    "B = np.full(imax, 1.0, dtype=float)\n",
    "A = h0*B\n",
    "Q = q*B\n",
    "\n",
    "Qup = Q[0]\n",
    "Adown = A[-1] \n",
    "\n",
    "tout = 600\n",
    "for it in range(int(totalTime/dt)):\n",
    "    A, Q, err = UnSteadyflowCollocated(A, Q, Adown, Qup, zb, B, dx, dt, g, manning)\n",
    "    if (dt*it)%tout==0: print((dt*it), err)\n",
    "        \n",
    "hc3 = A[::-1]/B[::-1] "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5692dc0",
   "metadata": {},
   "source": [
    "# fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ca04555",
   "metadata": {},
   "outputs": [],
   "source": [
    "import holoviews as hv\n",
    "hv.extension('bokeh')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "86dfee7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = hv.Curve((df.L,df.Z), label='Zb dx=0.5m').options(show_grid=True, ylabel='elevation[m]', xlabel='distance[m]' ,color='k') \\\n",
    "* hv.Scatter((Ls2, Zip2), label='Zb dx=10m').options(color='r',size=6) \\\n",
    "* hv.Scatter((Ls3, Zip3), label='Zb dx=100m').options(color='g',size=6) \\\n",
    "* hv.Curve((Ls1, Zip1 + hc1), label='dx=0.5m' ).options(color='b') \\\n",
    "* hv.Curve((Ls2, Zip2 + hc2), label='dx=10m').options(color='r') \\\n",
    "* hv.Curve((Ls3, Zip3 + hc3), label='dx=100m').options(color='g')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2a9b939",
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.extension('bokeh')\n",
    "figo = fig.options(width=500, legend_position='bottom_right')\n",
    "# d = hv.save(figo, 'USCdx.html')\n",
    "figo "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50a6b2db",
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.extension('bokeh')\n",
    "figo2 = fig.options(xlim=(3050, 3350), ylim=(7, 12), width=700, height=600, legend_position='right')\n",
    "# d = hv.save(figo2, 'USCdx2.html')\n",
    "figo2"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}