1D Unsteady Flow Noncontiuous ver. using jupyter .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous Equation\n",
    "\n",
    "$$\n",
    "    \\begin{align}\n",
    "        \\frac{\\partial A}{\\partial t}+\\frac{\\partial Q}{\\partial x} = 0\n",
    "    \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def conEq(A, Q, dt,dx, hmin):\n",
    "    nmax=len(A)\n",
    "    An = np.array([0.0 for i in range(nmax)])\n",
    "\n",
    "    for i in range(nmax):\n",
    "    \n",
    "        if i == nmax-1: # wall boundary\n",
    "            Qfp = 0.0\n",
    "        else:\n",
    "            Qc, Qip = Q[i], Q[i+1]\n",
    "            Qfp = Qc  if Qc > 0.0 and Qip > 0.0 else \\\n",
    "                 (Qip if Qc < 0.0 and Qip < 0.0 else 0.5*Qc+0.5*Qip )\n",
    "\n",
    "        if i == 0: # wall boundary\n",
    "            Qfm = 0.0\n",
    "        else:\n",
    "            Qc, Qim = Q[i], Q[i-1]\n",
    "            Qfm = Qim if Qc > 0.0 and Qim > 0.0 else \\\n",
    "                 (Qc  if Qc < 0.0 and Qim < 0.0 else 0.5*Qc+0.5*Qim )\n",
    "\n",
    "        An[i] = A[i] - dt * (Qfp - Qfm) / dx\n",
    "\n",
    "        if An[i] < hmin : An[i] = hmin\n",
    " \n",
    "    return An"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Momentum Equation\n",
    "\n",
    "$$\n",
    "    \\begin{align}\n",
    "        \\frac{\\partial Q}{\\partial t}+\\frac{\\partial }{\\partial x}\\left(\\frac{Q^2}{A}\\right)+gA\\frac{\\partial H}{\\partial x}+gAi_e = 0\n",
    "    \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "\n",
    "def momentEq(A, Q, An, zb, dt,dx, hmin):\n",
    "    gravity = 9.8\n",
    "    cManing = 0.02\n",
    "    nmax=len(Q)\n",
    "\n",
    "    Qn = np.array([0.0 for i in range(nmax)])\n",
    "\n",
    "    for i in range(1, nmax-1):\n",
    "        if An[i] <= hmin :\n",
    "            Qn[i] = 0.0\n",
    "        else:\n",
    "            Qc, Qip, Qim ,Qihp ,Qihm = Q[i], Q[i+1], Q[i-1], 0.5*(Q[i]+Q[i+1]), 0.5*(Q[i]+Q[i-1])\n",
    "            Ac, Aip, Aim ,Aihp ,Aihm = A[i], A[i+1], A[i-1], 0.5*(A[i]+A[i+1]), 0.5*(A[i]+A[i-1])\n",
    "            Hnc, Hnip, Hnim ,Hnihp ,Hnihm = An[i]+zb[i], An[i+1]+zb[i+1], An[i-1]+zb[i-1], \\\n",
    "            0.5*( An[i]+zb[i]+An[i+1]+zb[i+1] ), 0.5*( An[i]+zb[i]+An[i-1]+zb[i-1] )         \n",
    "\n",
    "            if (zb[i] > Hnip) or (zb[i+1] > Hnc):  #add\n",
    "                Q2byAfp = 0.0\n",
    "            else:\n",
    "                Q2byAfp = 0.0  if  Aip < hmin else \\\n",
    "                          Qc**2/Ac     if Qc > 0.0 and Qip > 0.0 else \\\n",
    "                         (Qip**2/Aip if Qc < 0.0 and Qip < 0.0 else \\\n",
    "                          Qihp**2/Aihp)\n",
    "\n",
    "            if (zb[i] > Hnim) or (zb[i-1] > Hnc): #add\n",
    "                Q2byAfm = 0.0\n",
    "            else:\n",
    "                Q2byAfm =  0.0  if  Aim < hmin else \\\n",
    "                          Qim**2/Aim   if Qc > 0.0 and Qim > 0.0 else \\\n",
    "                         (Qc**2/Ac     if Qc < 0.0 and Qim < 0.0 else \\\n",
    "                          Qihm**2/Aihm)           \n",
    "\n",
    "\n",
    "            if (zb[i] > Hnip) or (zb[i+1] > Hnc):  #add\n",
    "                dHdx =  ( Hnc-Hnim )/dx\n",
    "            elif (zb[i] > Hnim) or (zb[i-1] > Hnc): #add\n",
    "                dHdx =  ( Hnip-Hnc )/dx\n",
    "            else:\n",
    "\n",
    "                Vc = Qc/Ac\n",
    "                Vip = Qip/Aip if Aip > hmin else 0.0\n",
    "                Vim = Qim/Aim if Aim > hmin else 0.0\n",
    "\n",
    "                if   Qc > 0.0 and Qip > 0.0 and Qim > 0.0: \n",
    "                    Cr1, Cr2 = 0.5*( abs(Vc) + abs(Vip) )*dt/dx, 0.5*( abs(Vc) + abs(Vim) )*dt/dx\n",
    "                    dHdx1, dHdx2 = ( Hnip-Hnc )/dx, ( Hnc-Hnim )/dx\n",
    "\n",
    "                elif Qc < 0.0 and Qip < 0.0 and Qim < 0.0:\n",
    "                    Cr1, Cr2 = 0.5*( abs(Vc) + abs(Vim) )*dt/dx, 0.5*( abs(Vc) + abs(Vip) )*dt/dx\n",
    "                    dHdx1, dHdx2 = ( Hnc-Hnim )/dx, ( Hnip-Hnc )/dx\n",
    "\n",
    "                else:\n",
    "                    Cr1 = Cr2 = 0.5*( abs( 0.5*(Vc+Vip) ) + abs( 0.5*(Vc+Vim) ) )*dt/dx\n",
    "                    dHdx1 = dHdx2 = ( Hnihp - Hnihm )/dx\n",
    "\n",
    "                w1, w2 = 1.0 - Cr1, Cr2    \n",
    "                dHdx = w1 * dHdx1 + w2 * dHdx2\n",
    "            \n",
    "\n",
    "            Qn[i] = Q[i] - dt * (Q2byAfp - Q2byAfm) / dx \\\n",
    "                         - dt * gravity * An[i] * dHdx \\\n",
    "                         - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
    "\n",
    "#refrection Boundary\n",
    "    i = 0\n",
    "    Qc, Qip ,Qihp = Q[i], Q[i+1], 0.5*(Q[i]+Q[i+1])\n",
    "    Ac, Aip ,Aihp = A[i], A[i+1], 0.5*(A[i]+A[i+1])\n",
    "    Hnc, Hnip ,Hnihp = An[i]+zb[i], An[i+1]+zb[i+1], 0.5*(An[i]+zb[i]+An[i+1]+zb[i+1])         \n",
    "\n",
    "    Q2byAfp = Qc**2/Ac     if Qc > 0.0 and Qip > 0.0 else \\\n",
    "             (Qip**2/Aip   if Qc < 0.0 and Qip < 0.0 else \\\n",
    "              Qihp**2/Aihp)\n",
    "\n",
    "    Qn[i] = Q[i] - dt * (2.0*Q2byAfp) / dx \\\n",
    "                 - dt * gravity * An[i] * ((Hnip - Hnc)/dx) \\\n",
    "                 - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
    "        \n",
    "    i = nmax-1\n",
    "    Qc, Qim, Qihm = Q[i], Q[i-1], 0.5*(Q[i]+Q[i-1])\n",
    "    Ac, Aim, Aihm = A[i], A[i-1], 0.5*(A[i]+A[i-1])\n",
    "    Hnc, Hnim, Hnihm = An[i]+zb[i], An[i-1]+zb[i-1], 0.5*(An[i]+zb[i]+An[i-1]+zb[i-1])         \n",
    "\n",
    "    Q2byAfm = Qim**2/Aim   if Qc > 0.0 and Qim > 0.0 else \\\n",
    "             (Qc**2/Ac     if Qc < 0.0 and Qim < 0.0 else \\\n",
    "              Qihm**2/Aihm)           \n",
    "    \n",
    "    Qn[i] = Q[i] - dt * (-2.0*Q2byAfm) / dx \\\n",
    "                 - dt * gravity * An[i] * ((Hnc - Hnim)/dx) \\\n",
    "                 - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
    "\n",
    "    return Qn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulation(A, Q, zb, dt,dx,  hmin):\n",
    "    An = conEq(A, Q, dt,dx, hmin)\n",
    "    Qn = momentEq(A, Q, An, zb, dt,dx, hmin)\n",
    "    return An, Qn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#%%timeit -r1\n",
    "#%matplotlib nbagg\n",
    "%matplotlib inline\n",
    "\n",
    "hmin = 10.0**(-4)\n",
    "hu = 0.4\n",
    "hd = hmin\n",
    "nmax = 150\n",
    "ng = 20\n",
    "dx = 0.1\n",
    "dt = 0.0005\n",
    "dtout= 0.2\n",
    "tmax = 40.0\n",
    "\n",
    "Q =  np.array([0.0 for i in range(nmax)])\n",
    "# zb = np.array([0.0 for i in range(nmax)])\n",
    "# zb = np.array([0.0 if i < 40 else (i-40)*dx*0.1 for i in range(nmax)])\n",
    "zb = np.array([0.0 if i < 50 else 0.5 for i in range(nmax)])\n",
    "Lx = np.array([i*dx for i in range(nmax)])\n",
    "\n",
    "# Initial Condition\n",
    "A =  np.array([hd for i in range(nmax)])\n",
    "A[120:] = hu\n",
    "\n",
    "#gragh \n",
    "fig, ax = plt.figure(figsize=(8, 6), dpi=80), plt.subplot()\n",
    "\n",
    "def timeGenerater():\n",
    "    it, itout = 0, 1\n",
    "    isTiter, isOut = True, False\n",
    "    yield it*dt, isTiter, isOut\n",
    "\n",
    "    while True:\n",
    "        it += 1\n",
    "        if it*dt >= itout*dtout:\n",
    "            itout += 1\n",
    "            isOut = True\n",
    "        else:\n",
    "            isOut = False\n",
    "\n",
    "        if it*dt >= tmax: isTiter = False\n",
    "\n",
    "        yield it*dt, isTiter, isOut\n",
    "\n",
    "        \n",
    "tgen = timeGenerater() \n",
    "t, isTiter, isOut = next(tgen)\n",
    "\n",
    "ims=[]\n",
    "\n",
    "fr = ax.plot(Lx, A+zb, 'r')\n",
    "zbf = ax.plot(Lx, zb, 'k')\n",
    "\n",
    "title = ax.set_title( 't = ' + str(round(t,3)) + ', Volume = '  + str(sum(A*dx)) )\n",
    "ax.set_ylim(0, 1.0, 0.1)\n",
    "ax.set_xlim(0, dx*nmax, 1)\n",
    "plt.tight_layout\n",
    "plt.savefig( str( round(t*1000) ).zfill(6) + \".png\" , bbox_inches='tight')\n",
    "ax.cla()\n",
    "\n",
    "while isTiter:\n",
    "    A, Q = simulation(A, Q, zb, dt,dx, hmin)\n",
    "    t, isTiter, isOut = next(tgen)\n",
    "      \n",
    "    if isOut:\n",
    "        title = ax.set_title('t = ' + str(round(t,3))+', Volume = '+ str(sum(A*dx)) )\n",
    "        fr = ax.plot(Lx, A+zb,'r' )\n",
    "        zbf = ax.plot(Lx, zb, 'k' )\n",
    "        ax.set_ylim(0, 1.0, 0.1)\n",
    "        ax.set_xlim(0, dx*nmax, 1)\n",
    "        plt.tight_layout\n",
    "        plt.savefig( str( round(t*1000) ).zfill(6) + \".png\" , bbox_inches='tight')\n",
    "        ax.cla()\n",
    "\n",
    "        \n",
    "#ims.append(fr) #imagemagick\n",
    "#ani = animation.ArtistAnimation(fig, ims)\n",
    "#ani.save('out.gif', writer='imagemagick', fps=5) "
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.0"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}