2D Unsteady Flow using jupyter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.animation as animation\n",
    "from mpl_toolkits.mplot3d import Axes3D"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Continuous Equation\n",
    "\n",
    "$$\n",
    "    \\begin{align}\n",
    "        \\frac{\\partial h}{\\partial t}+\\frac{\\partial q_x}{\\partial x} +\\frac{\\partial q_y}{\\partial y} = 0\n",
    "    \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def conEq(dep, qx, qy, dt, dx, dy, hmin):\n",
    "    nxmax, nymax = len(dep), len(dep[0])\n",
    "    depn = np.array([[ 0.0 for j in range(nymax)] for i in range(nxmax)])\n",
    "\n",
    "    fluxp = lambda qc, qp : qc  if qc > 0.0 and qp > 0.0 else \\\n",
    "                           (qp  if qc < 0.0 and qp < 0.0 else 0.5*qc+0.5*qp )\n",
    "    fluxm = lambda qc, qm : qm  if qc > 0.0 and qm > 0.0 else \\\n",
    "                           (qc  if qc < 0.0 and qm < 0.0 else 0.5*qc+0.5*qm )\n",
    "\n",
    "    for i in range(nxmax):\n",
    "        for j in range(nymax):\n",
    "            qxfp = 0.0 if i == nxmax-1 else fluxp(qx[i][j], qx[i+1][j  ])\n",
    "            qxfm = 0.0 if i == 0       else fluxm(qx[i][j], qx[i-1][j  ])\n",
    "            qyfp = 0.0 if j == nymax-1 else fluxp(qy[i][j], qy[i  ][j+1])\n",
    "            qyfm = 0.0 if j == 0       else fluxm(qy[i][j], qy[i  ][j-1])\n",
    "\n",
    "            depn[i][j] = dep[i][j] - dt*(qxfp - qxfm)/dx - dt*(qyfp - qyfm)/dy\n",
    "            if depn[i][j] < hmin : depn[i][j] = hmin\n",
    "\n",
    "    return depn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Momentum Equation\n",
    "\n",
    "$$\n",
    "    \\begin{align}\n",
    "        \\frac{\\partial q_x}{\\partial t}+\\frac{\\partial u q_x}{\\partial x}+\\frac{\\partial v q_x}{\\partial y}+gh\\frac{\\partial H}{\\partial x}+\\frac{\\tau_{0x}}{\\rho} \n",
    "        - \\nu_t h \\left(\\frac{\\partial^2 u}{\\partial x^2}+\\frac{\\partial^2 u}{\\partial y^2} \\right)= 0 \\\\\n",
    "        \\frac{\\partial q_y}{\\partial t}+\\frac{\\partial u q_y}{\\partial x}+\\frac{\\partial v q_y}{\\partial y}+gh\\frac{\\partial H}{\\partial y}+\\frac{\\tau_{0y}}{\\rho}- \\nu_t h \\left(\\frac{\\partial^2 v}{\\partial x^2}+\\frac{\\partial^2 v}{\\partial y^2} \\right)\n",
    "        = 0\n",
    "    \\end{align}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def momentEq(dep, qx, qy, depn, zb, dt, dx, dy, hmin, direction):\n",
    "    #direction = 1:x, 2:y\n",
    "    gravity = 9.8\n",
    "    cManing = 0.012\n",
    "\n",
    "    #q = qx if direction == 1 else qy\n",
    "    if direction == 1 :\n",
    "        q = qx\n",
    "    elif direction == 2 :\n",
    "        q = qy\n",
    "    else:\n",
    "        print('error')\n",
    "        \n",
    "    nxmax, nymax = len(q), len(q[0])\n",
    "    qn = np.array([[ 0.0 for j in range(nymax)] for i in range(nxmax)])\n",
    "\n",
    "    fluxp = lambda qc,qp,Vc,Vp : qc*Vc if Vc > 0.0 and Vp > 0.0 else \\\n",
    "                                 (qp*Vp if Vc < 0.0 and Vp < 0.0 else 0.5*(qc+qp)*0.5*(Vc+Vp))\n",
    "    fluxm = lambda qc,qm,Vc,Vm : qm*Vm if Vc > 0.0 and Vm > 0.0 else \\\n",
    "                                 (qc*Vc if Vc < 0.0 and Vm < 0.0 else 0.5*(qc+qm)*0.5*(Vc+Vm))\n",
    "\n",
    "\n",
    "\n",
    "    for i in range(nxmax):\n",
    "        for j in range(nymax):    \n",
    "            if depn[i][j] <= hmin*2 :\n",
    "                qn[i][j] = 0.0\n",
    "            else:\n",
    "\n",
    "            # advection term\n",
    "                if i == nxmax-1 : #refrection\n",
    "                    qc, qm  =   q[i][j], q[i-1][j  ]\n",
    "                    uc, um  =  qx[i][j]/dep[i][j], qx[i-1][j  ]/dep[i-1][j  ]                                \n",
    "                    uqfm = fluxm(qc,qm,uc,um)\n",
    "                    uqfp = -uqfm\n",
    "                elif i == 0 :    #refrection\n",
    "                    qc,qp = q[i][j], q[i+1][j]\n",
    "                    uc,up = qx[i][j]/dep[i][j],qx[i+1][j]/dep[i+1][j]\n",
    "                    uqfp = fluxp(qc,qp,uc,up)\n",
    "                    uqfm = -uqfp\n",
    "                else :\n",
    "                    qc,qp,qm = q[i][j], q[i+1][j], q[i-1][j]\n",
    "                    uc,up,um = qx[i][j]/dep[i][j],qx[i+1][j]/dep[i+1][j],qx[i-1][j]/dep[i-1][j]\n",
    "                    uqfp = fluxp(qc,qp,uc,up)\n",
    "                    uqfm = fluxm(qc,qm,uc,um)\n",
    "                if j == nymax-1 : #refrection\n",
    "                    qc,qm =  q[i][j], q[i][j-1]\n",
    "                    vc,vm = qy[i][j]/dep[i][j], qy[i  ][j-1]/dep[i  ][j-1]\n",
    "                    vqfm = fluxm(qc,qm,vc,vm)\n",
    "                    vqfp = -vqfm\n",
    "                elif j == 0 :    #refrection\n",
    "                    qc,qp =  q[i][j], q[i][j+1]\n",
    "                    vc,vp = qy[i][j]/dep[i][j], qy[i  ][j+1]/dep[i  ][j+1]\n",
    "                    vqfp = fluxp(qc,qp,vc,vp)\n",
    "                    vqfm = -vqfp\n",
    "                else :\n",
    "                    qc,qp,qm =  q[i][j], q[i][j+1], q[i][j-1]\n",
    "                    vc,vp,vm = qy[i][j]/dep[i][j], qy[i  ][j+1]/dep[i  ][j+1], qy[i  ][j-1]/dep[i  ][j-1]\n",
    "                    vqfp = fluxp(qc,qp,vc,vp)\n",
    "                    vqfm = fluxm(qc,qm,vc,vm)\n",
    "\n",
    "    # pressure & gravity term\n",
    "                if direction ==1 :\n",
    "                    ib, pb, mb, delta,idr,jdr = i, nxmax-1, 0, dx, 1, 0 \n",
    "                else : \n",
    "                    ib, pb, mb, delta,idr,jdr = j, nymax-1, 0, dy, 0, 1\n",
    "                    \n",
    "                if ib == pb :\n",
    "                    Hc, Hm = depn[i][j]+zb[i][j], depn[i-idr][j-jdr]+zb[i-idr][j-jdr]\n",
    "                    dHdx = (Hc-Hm)/delta  \n",
    "                elif ib == mb :\n",
    "                    Hc, Hp = depn[i][j]+zb[i][j], depn[i+idr][j+jdr]+zb[i+idr][j+jdr]\n",
    "                    dHdx = (Hp-Hc)/delta  \n",
    "                else :\n",
    "                    Vc, Vp, Vm = q[i][j]/dep[i][j], q[i+idr][j+jdr]/dep[i+idr][j+jdr], q[i-idr][j-jdr]/dep[i-idr][j-jdr]\n",
    "                    Hc, Hp, Hm = depn[i][j]+zb[i][j], depn[i+idr][j+jdr]+zb[i+idr][j+jdr], depn[i-idr][j-jdr]+zb[i-idr][j-jdr]\n",
    "\n",
    "                    if Vc > 0.0 and Vp > 0.0 and Vm > 0.0: \n",
    "                        Cr1, Cr2 = 0.5*(abs(Vc)+abs(Vp))*dt/delta, 0.5*(abs(Vc)+abs(Vm))*dt/delta\n",
    "                        dHdx1, dHdx2 = (Hp-Hc)/delta, (Hc-Hm)/delta\n",
    "                    elif Vc < 0.0 and Vp < 0.0 and Vm < 0.0:\n",
    "                        Cr1, Cr2 = 0.5*(abs(Vc)+abs(Vm))*dt/delta, 0.5*(abs(Vc)+abs(Vp))*dt/delta\n",
    "                        dHdx1, dHdx2 = (Hc-Hm)/delta, (Hp-Hc)/delta           \n",
    "                    else:\n",
    "                        Cr1 = Cr2 = 0.5*(abs(0.5*(Vc+Vp))+abs(0.5*(Vc+Vm)))*dt/delta\n",
    "                        dHdx1 = dHdx2 = (0.5*(Hc+Hp) - 0.5*(Hc+Hm)) / delta\n",
    "\n",
    "                    w1, w2 = 1-Cr1**0.5, Cr2**0.5\n",
    "                    dHdx = w1 * dHdx1 + w2 * dHdx2   \n",
    "\n",
    "#viscous sublayer\n",
    "                Cf = gravity*cManing**2.0/dep[i][j]**(1.0/3.0) \n",
    "                u, v = qx[i][j]/dep[i][j],qy[i][j]/dep[i][j] \n",
    "                Vnorm = np.sqrt(u**2.0+v**2.0) \n",
    "                if direction ==1 :\n",
    "                    Vis = Cf * Vnorm * u   #Vnorm**2*u/Vnorm\n",
    "                else:\n",
    "                    Vis = Cf * Vnorm * v\n",
    "#turbulence\n",
    "#                kenergy = 2.07*Cf*Vnorm**2\n",
    "                nut = 0.4/6.0*dep[i][j]*np.sqrt(Cf)*np.abs(Vnorm)\n",
    "    \n",
    "                if i == nxmax-1 :\n",
    "                    turbX = nut * (-(q[i][j]/dep[i][j]-q[i-1][j]/dep[i-1][j]))/dx**2.0\n",
    "                elif i == 0 :\n",
    "                    turbX = nut * ((q[i+1][j]/dep[i+1][j]-q[i][j]/dep[i][j]))/dx**2.0\n",
    "                else :\n",
    "                    turbX = nut * ((q[i+1][j]/dep[i+1][j]-q[i][j]/dep[i][j]) \\\n",
    "                                 -(q[i][j]/dep[i][j]-q[i-1][j]/dep[i-1][j]) \\\n",
    "                                 )/dx**2.0\n",
    "\n",
    "                if j == nymax-1 :\n",
    "                    turbY = nut * ( 0.0 - (q[i][j]/dep[i][j]-q[i][j-1]/dep[i][j-1]))/dy**2.0\n",
    "                elif j==0 :\n",
    "                    turbY = nut * ((q[i][j+1]/dep[i][j+1]-q[i][j]/dep[i][j]))/dy**2.0\n",
    "                else :\n",
    "                    turbY = nut * ((q[i][j+1]/dep[i][j+1]-q[i][j]/dep[i][j]) \\\n",
    "                                 -(q[i][j]/dep[i][j]-q[i][j-1]/dep[i][j-1]) \\\n",
    "                                 )/dy**2.0\n",
    "                    \n",
    "                qn[i][j] = q[i][j] - dt * (uqfp - uqfm) / dx \\\n",
    "                                   - dt * (vqfp - vqfm) / dy \\\n",
    "                                   - dt * gravity * depn[i][j] * dHdx \\\n",
    "                                   - dt * Vis \\\n",
    "                                   + dt * dep[i][j] * (turbX+turbY)\n",
    "\n",
    "    return qn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def simulation(dep, qx, qy, zb, dt, dx, dy, hmin):\n",
    "    depn = conEq(dep, qx, qy, dt, dx, dy, hmin)\n",
    "    qxn = momentEq(dep, qx, qy, depn, zb, dt, dx, dy, hmin, 1)\n",
    "    qyn = momentEq(dep, qx, qy, depn, zb, dt, dx, dy, hmin, 2)\n",
    "\n",
    "    nxmax, nymax = len(depn), len(depn[0])\n",
    "    \n",
    "    CFLmin = 0.0\n",
    "    for i in range(nxmax):\n",
    "        for j in range(nymax):    \n",
    "            if depn[i][j] > hmin :\n",
    "                CFLmin = np.max( \\\n",
    "                                [((np.abs(qxn[i][j]/depn[i][j])+np.sqrt(9.8*depn[i][j]))/dx \\\n",
    "                                + (np.abs(qyn[i][j]/depn[i][j])+np.sqrt(9.8*depn[i][j]))/dy)*dt \\\n",
    "                                ,CFLmin])\n",
    "    return depn, qxn, qyn, CFLmin     "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Main routine"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%timeit -r1\n",
    "\n",
    "#%matplotlib nbagg\n",
    "#%matplotlib inline\n",
    "\n",
    "hmin = 10.0**(-4)\n",
    "nmax = 101\n",
    "ng = 50\n",
    "hu = 1.0\n",
    "hd = hmin\n",
    "dx = 0.1\n",
    "dy = 0.1\n",
    "dt = 0.005\n",
    "dtout= 0.1\n",
    "tmax = 20.0\n",
    "\n",
    "qx = np.array([[0.0 for j in range(nmax)] for i in range(nmax)])\n",
    "qy = np.array([[0.0 for j in range(nmax)] for i in range(nmax)])\n",
    "zb = np.array([[0.0 for j in range(nmax)] for i in range(nmax)])\n",
    "\n",
    "Lx = np.array([i*dx for i in range(nmax)])\n",
    "Ly = np.array([i*dy for i in range(nmax)])\n",
    "\n",
    "# Initial Condition\n",
    "dep = np.array([[hd for j in range(nmax)] for i in range(nmax)])\n",
    "for i in range(40,61):\n",
    "    for j in range(40,61):\n",
    "        dep[i][j] = hu\n",
    "\n",
    "#gragh \n",
    "fig, ax = plt.figure(figsize=(8, 6), dpi=150), plt.subplot(projection='3d')\n",
    "#fig, ax = plt.figure(figsize=(4, 3), dpi=72), plt.subplot(projection='3d')\n",
    "\n",
    "ax.set_zlim3d(0, 1.0)\n",
    "\n",
    "X, Y = np.meshgrid(Lx, Ly)\n",
    "deps = list(map(list,zip(*dep)))\n",
    "surf = ax.plot_surface(X, Y, deps, cstride = 1, rstride = 1,linewidth=0.1)\n",
    "fig.suptitle('t = 0' )\n",
    "plt.savefig('H000.png')\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",
    "nout= 1\n",
    "while isTiter:\n",
    "    dep, qx, qy, CFLmin = simulation(dep, qx, qy, zb, dt, dx, dy, hmin)\n",
    "    t, isTiter, isOut = next(tgen)\n",
    "    if isOut:\n",
    "        surf.remove()\n",
    "        deps = list(map(list,zip(*dep)))\n",
    "        surf = ax.plot_surface(X, Y, deps, cstride = 1, rstride = 1,linewidth=0.1)\n",
    "        fig.suptitle('t = ' + str(round(t,2)) +' CFLmax =' + str(CFLmin))\n",
    "        plt.savefig('H'+str(nout).zfill(3)+'.png')\n",
    "        nout += 1\n",
    "\n",
    "#plt.show()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}