{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import S1DBedloadSolver as S1Dbed\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 second\n",
      "720000 second\n",
      "1440000 second\n",
      "2160000 second\n",
      "2880000 second\n",
      "3600000 second\n",
      "4320000 second\n",
      "5040000 second\n",
      "5760000 second\n",
      "6480000 second\n",
      "7200000 second\n",
      "Wall time: 2min 28s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "length = 10000.0\n",
    "dx = 100.0\n",
    "\n",
    "imax = int(length/dx) + 1\n",
    "dt = 10.0\n",
    "totalTime = 2000.1*3600.0\n",
    "outTimeStep = 200.0*3600.0\n",
    "RunUpTime = 3.0 * 3600.0\n",
    "hini = 1.0\n",
    "manning  = 0.03\n",
    "ib = 1.0/200.0\n",
    "ib2 = 1.0/200.0\n",
    "outputfilename = '1D_case3.json'\n",
    "\n",
    "# grain diameter classification\n",
    "screenclass = np.array( [ 2.0**1, 2.0**2, 2.0**3, 2.0**4, 2.0**5, 2.0**6, 2.0**7], dtype=float )/1000\n",
    "dsize = np.array( [ 10**(0.5*(np.log10(screenclass[i]) + np.log10(screenclass[i+1]))) for i in range(len(screenclass)-1) ], dtype=float )\n",
    "\n",
    "# percentage of grain size under exchange layer\n",
    "dmax = screenclass[-1]\n",
    "P = (screenclass/dmax)**1\n",
    "dratioStandard1 = P[1:] - P[:-1]\n",
    "dratioStandard1 = np.full_like(dsize, 1/len(dsize), dtype=float)\n",
    "dratioStandard = np.full( (imax, len(dsize) ), dratioStandard1, dtype=float)\n",
    "dratioStandard[0,:] = 0.0\n",
    "dratioStandard[0,-1] = 1.0\n",
    "\n",
    "# initial percentage of grain size in exchange layer\n",
    "dratio = np.copy(dratioStandard)\n",
    "\n",
    "# thickness of exchange layer \n",
    "hExlayer = dsize[-1]\n",
    "\n",
    "# Initial & Boundary condition\n",
    "B = np.full(imax, 1.0, dtype=float)\n",
    "A = hini*B\n",
    "Q = ib**0.5*(hini)**(5.0/3.0)/manning*B #normal flow\n",
    "zb = np.zeros(imax)\n",
    "for i in range(1,imax):\n",
    "    zb[i] = zb[i-1] + ib2*dx if i < 50 else zb[i-1] + ib*dx\n",
    "    \n",
    "zb = zb[::-1]\n",
    "# WLdown = 2.0\n",
    "# for i, (Ap, zbp) in enumerate(zip(A, zb)):\n",
    "#     if( Ap/B[i] + zbp < WLdown) : A[i] = (WLdown - zbp)*B[i]\n",
    "\n",
    "dAb = np.zeros(imax)\n",
    "Qup = Q[0]\n",
    "def Adown(time, Q, dzb, ib):\n",
    "#     return (WLdown - (dzb + zb[-1]))*B[-1]\n",
    "    return ( manning**2*Q**2/ib/B[-1]**2 )**0.3 * B[-1]\n",
    "\n",
    "def Qup(time):\n",
    "    return Q[0]\n",
    " \n",
    "S1Dbed.bedvariation(\n",
    "dx,dt,manning,totalTime,outTimeStep,RunUpTime\n",
    ",dsize ,dratioStandard ,dratio ,hExlayer ,A ,Q ,B ,zb ,dAb ,Qup ,Adown\n",
    ",outputfilename, screenclass\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.6.4"
  },
  "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": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "165px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}