1D Unsteady Flow Discontiuous ver. using jupyter .ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "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": [
    "# judge discontinous boundary\n",
    " - set wall boundary : search discontinous boundary\n",
    " - simply check whether the water surface is conected between cells　"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def isDiscontinous(zbp, zbm, Ap, Am):\n",
    "    if (( zbp < zbm ) and ((zbp + Ap) < zbm)): \n",
    "        isInterBound = True\n",
    "    elif (( zbp > zbm ) and ( zbp > (zbm + Am))): \n",
    "        isInterBound = True \n",
    "    else :\n",
    "        isInterBound = False\n",
    "\n",
    "    return isInterBound"
   ]
  },
  {
   "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": null,
   "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",
    "    judgeDiscontinous = True\n",
    "\n",
    "    flux = np.array([0.0 for i in range(nmax+1)])\n",
    "\n",
    "#set flux\n",
    "    for i in range(1, nmax):\n",
    "        Qp, Qm = Q[i], Q[i-1]\n",
    "        flux[i] = Qm  if Qp >= 0.0 and Qm >= 0.0 else (Qp if Qp <= 0.0 and Qm <= 0.0 else 0.5*Qp+0.5*Qm )\n",
    "\n",
    "#add calculation of discontinous boundary\n",
    "        if judgeDiscontinous:\n",
    "            if isDiscontinous( zb[i], zb[i-1], A[i], A[i-1] ) :\n",
    "                Qp, Qm = Q[i], Q[i-1]\n",
    "                if(zb[i] < zb[i-1]):\n",
    "                    flux[i] = Qm  if Qm >= 0.0  else  0.0\n",
    "                elif(zb[i] > zb[i-1]):\n",
    "                    flux[i] = Qp  if Qp <= 0.0  else  0.0\n",
    "\n",
    "# set boundary flux\n",
    "    flux[0], flux[nmax] = 0.0, 0.0\n",
    "    \n",
    "    for i in range(nmax):\n",
    "    \n",
    "        An[i] = A[i] - dt * (flux[i+1] - flux[i]) / dx\n",
    "        if An[i] <= hmin :\n",
    "            # if An[i] <= hmin*0.99 : print('out',i , A[i], A[i+1], An[i], Q[i], flux[i+1], flux[i])\n",
    "            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": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def momentEq(A, Q, An, zb, dt,dx, hmin):\n",
    "    gravity = 9.8\n",
    "    cManing = 0.02\n",
    "    nmax=len(Q)\n",
    "    judgeDiscontinous = True\n",
    "    \n",
    "    Qn = np.array([0.0 for i in range(nmax)])\n",
    "    flux = np.array([0.0 for i in range(nmax+1)])\n",
    "\n",
    "    def caldHdx(Qc,Qip,Qim,Ac,Aip,Aim,Hnc,Hnip,Hnim):\n",
    "\n",
    "        Qihp, Qihm = 0.5*(Qc+Qip), 0.5*(Qc+Qim)\n",
    "        Aihp, Aihm = 0.5*(Ac+Aip), 0.5*(Ac+Aim)\n",
    "        Hnihp, Hnihm = 0.5*(Hnc+Hnip), 0.5*(Hnc+Hnim)\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",
    "        return dHdx\n",
    "\n",
    "#set flux\n",
    "    for i in range(1, nmax):\n",
    "        Qp, Qm, Qh  = Q[i], Q[i-1], 0.5*(Q[i] + Q[i-1])\n",
    "        Ap, Am, Ah  = A[i], A[i-1], 0.5*(A[i] + A[i-1])\n",
    "\n",
    "        flux[i]= Qm**2/Am if Qp >= 0.0 and Qm >= 0.0 else (Qp**2/Ap  if Qp <= 0.0 and Qm <= 0.0 else Qh**2/Ah )\n",
    "\n",
    "#add calculation of discontinous boundary\n",
    "#In discontinous boundary like a weir, momentum flux sets 0 for free wall in upper and for refrection boundary in downer as a result\n",
    "        if judgeDiscontinous:\n",
    "            if isDiscontinous( zb[i], zb[i-1], An[i], An[i-1] ) : flux[i] = 0.0\n",
    "\n",
    "# set boundary flux : refrection\n",
    "            flux[0], flux[nmax] = 0.0, 0.0\n",
    "\n",
    "\n",
    "    for i in range(nmax):\n",
    "        if An[i] <= hmin :\n",
    "            Qn[i] = 0.0\n",
    "        else:\n",
    "\n",
    "            if i == 0: #refrection\n",
    "                dHdx = caldHdx( Q[i], Q[i+1], -Q[i] \\\n",
    "                            , A[i], A[i+1], A[i] \\\n",
    "                            , An[i]+zb[i], An[i+1]+zb[i+1], An[i]+zb[i] )\n",
    "            elif i == nmax-1:\n",
    "                dHdx = caldHdx( Q[i], -Q[i], Q[i-1] \\\n",
    "                            , A[i], A[i], A[i-1] \\\n",
    "                            , An[i]+zb[i], An[i]+zb[i], An[i-1]+zb[i-1] )\n",
    "            else:\n",
    "                dHdx = caldHdx( Q[i], Q[i+1], Q[i-1] \\\n",
    "                            , A[i], A[i+1], A[i-1] \\\n",
    "                            , An[i]+zb[i], An[i+1]+zb[i+1], An[i-1]+zb[i-1] )\n",
    "\n",
    "# set wall boundary : search discontinous boundary\n",
    "            if judgeDiscontinous :\n",
    "                #boundConp and boundConm is 0:normal, 1:wall, 2:free fall \n",
    "                if i == nmax-1 :\n",
    "                    boundConp = 1\n",
    "                elif isDiscontinous( zb[i+1], zb[i], An[i+1], An[i] ) :\n",
    "                    boundConp = 1 if( zb[i] < zb[i+1] ) else 2 \n",
    "                else:\n",
    "                    boundConp = 0\n",
    "\n",
    "                if i == 0 :\n",
    "                    boundConm = 1\n",
    "                elif isDiscontinous( zb[i], zb[i-1], An[i], An[i-1] ):\n",
    "                    boundConm = 1 if( zb[i] < zb[i-1] ) else 2 \n",
    "                else:\n",
    "                    boundConm = 0\n",
    "\n",
    "                if boundConp == 0 and boundConm == 0 :\n",
    "                    dHdx = dHdx\n",
    "                elif boundConp == 1 and boundConm == 1 :\n",
    "                    dHdx = 0.0\n",
    "                elif boundConp == 2 and boundConm == 2 :\n",
    "                    dHdx = 0.0\n",
    "                elif boundConp == 1 : #wall\n",
    "                    dHdx = caldHdx( Q[i], -Q[i], Q[i-1] \\\n",
    "                            , A[i], A[i], A[i-1] \\\n",
    "                            , An[i]+zb[i], An[i]+zb[i], An[i-1]+zb[i-1] )\n",
    "                elif boundConp == 2 : #free wall\n",
    "                    dHdx =  ( 0.0 - 0.5*(An[i]+An[i-1]) )/dx + ( zb[i]-zb[i-1] )/dx\n",
    "                elif boundConm == 1 : #wall\n",
    "                    dHdx = caldHdx( Q[i], Q[i+1], -Q[i] \\\n",
    "                            , A[i], A[i+1], A[i] \\\n",
    "                            , An[i]+zb[i], An[i+1]+zb[i+1], An[i]+zb[i] )\n",
    "                elif boundConm == 2 : #free wall\n",
    "                    dHdx =  ( 0.5*(An[i]+An[i+1]) - 0.0 )/dx + ( zb[i+1] - zb[i] )/dx\n",
    "                else :\n",
    "                    print('error')\n",
    "\n",
    "            advection = (flux[i+1]- flux[i]) / dx\n",
    "\n",
    "            Qn[i] = Q[i] - dt * advection - dt * gravity * An[i] * dHdx - dt * gravity * cManing**2 * Q[i] * abs(Q[i]) / A[i] / A[i]**(4 / 3)\n",
    "\n",
    "    return Qn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": [
    "%matplotlib inline\n",
    "\n",
    "hmin = 10.0**(-5)\n",
    "hu, hd = 0.6, hmin\n",
    "nmax, ng = 100, 20\n",
    "dx, dt = 0.1, 0.0005\n",
    "dtout= 0.4\n",
    "tmax = 60.0\n",
    "\n",
    "Q =  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.4 if i > 80 else 0.0 for i in range(nmax)])\n",
    "zb[:20] =0.15 \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[90:] = 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",
    "tgen = timeGenerater() \n",
    "t, isTiter, isOut = next(tgen)\n",
    "\n",
    "fr, zbf = ax.plot(Lx, A+zb, 'r'), ax.plot(Lx, zb, 'k')\n",
    "title = ax.set_title( 't = ' + str(round(t, 3)) + ', Volume = '  + str(sum(A*dx)) )\n",
    "yl, xl = ax.set_ylim(0, 1.0, 0.1), 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",
    "        print('t=',t)\n",
    "        title = ax.set_title('t = ' + str(round(t,3))+', Volume = '+ str(sum(A*dx)) )\n",
    "        eng = (Q/A)**2/2/9.8+A+zb\n",
    "        fr, zbf = ax.plot(Lx, A+zb, 'r'), ax.plot(Lx, zb, 'k')\n",
    "        yl, xl = ax.set_ylim(0, 1.0, 0.1), 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()"
   ]
  }
 ],
 "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
}