{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nOne dimension shallow water, steady lake validation case.\n=========================================================\n\nA steady lake is modelised by shallow water equation.\n\nThe model reads as\n\n\\begin{align}\\begin{cases}\n        \\frac{\\partial h}{\\partial t} + \\frac{\\partial h}{\\partial x} &= 0 \\\\\n        \\frac{\\partial u}{\\partial t} + u\\,\\frac{\\partial u}{\\partial x} &= -g\\,\\frac{\\partial h}{\\partial x}\n    \\end{cases}\\end{align}\n\nThe initial condition is a uniform fluid altitude on a curvy bottom. We expect the fluid altitude to remain\nconstant as well as the total volume of fluid.\n\nThe finite difference being unable ensure quantity conservation, this very simple case\ncan fail as the numerical error stacks, leading to wave that appears and propagate.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import pylab as pl\nimport numpy as np\nfrom skfdiff import Model\nfrom scipy.integrate import trapz\n\nshallow = Model(\n    [\"-dxq\", \"-(g*h*dxh + 2*q*dxq/h - q**2*dxh/h**2 + g * h * dxZ - k * dxxh) * eps\"],\n    [\"h(x)\", \"q(x)\"],\n    [\"g\", \"k\", \"eps(x)\", \"Z(x)\"],\n    boundary_conditions=\"periodic\",\n)\n\nx, dx = np.linspace(-10, 8, 500, retstep=True)\nZ = x ** 2 * np.sin(x) + 3 * x + 80\n\nh = 80 - Z\n\nini_fields = shallow.Fields(x=x, h=h, q=x * 0, Z=Z, g=9.81, k=0, eps=np.ones(x.size))"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This hook ensure proper computation in dry places.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def dry_hook(t, fields):\n    fields[\"eps\"] = \"x\", np.where(fields.h <= 1, 0, 1)\n    fields[\"q\"] = \"x\", np.where(fields.h <= 0, 0, fields.q)\n    return fields\n\n\nsimul = shallow.init_simulation(\n    fields=ini_fields,\n    hook=dry_hook,\n    t=0,\n    dt=60,\n    tmax=3600 * 3,\n    id=\"steady_lake\",\n    tol=1e-1,\n)\n\ncontainer = simul.attach_container()\nsimul.run()\n\ndata = container.data.copy()\ndata[\"h\"] = data[\"h\"].where(data.h > 0)\ndata[\"q\"] = data[\"q\"].where(data.h > 0)\ndata[\"eta\"] = data.h + data.Z\ndata[\"V\"] = \"t\", trapz(data.h.where(data.h >= 0, other=0), x=data.x)\n\nfig = pl.figure(figsize=(5, 3))\n\nax = pl.subplot2grid((2, 2), (0, 0))\npl.sca(ax)\ndata[\"eta\"].isel(t=0).plot(color=\"navy\")\npl.fill_between(\n    data.x, data[\"eta\"].isel(t=0), data[\"Z\"].isel(t=0), color=\"navy\", alpha=0.3\n)\ndata[\"Z\"].isel(t=0).plot(color=\"black\")\npl.xlim(data[\"x\"].min(), data[\"x\"].max())\npl.title(\"$t=0$ (hours)\")\n\nax = pl.subplot2grid((2, 2), (0, 1))\npl.sca(ax)\ndata[\"eta\"].isel(t=-1).plot(color=\"navy\")\ndata[\"Z\"].isel(t=-1).plot(color=\"black\")\npl.fill_between(\n    data.x, data[\"eta\"].isel(t=-1), data[\"Z\"].isel(t=-1), color=\"navy\", alpha=0.3\n)\npl.xlim(data[\"x\"].min(), data[\"x\"].max())\npl.title(\"$t=3$ (hours)\")\n\nax = pl.subplot2grid((2, 2), (1, 0), colspan=2)\npl.sca(ax)\npl.plot(data[\"t\"] / 3600, data[\"V\"], color=\"black\")\npl.ylim(242.7621835042269, 242.7621835042272)\npl.xlim((data[\"t\"] / 3600).min(), (data[\"t\"] / 3600).max())\npl.title(r\"Total volume ($V \\approx 242.76$, $\\Delta V \\approx 3.1e-13$)\")\npl.xlabel(\"t (hour)\")\npl.ylabel(\"V\")\n\npl.tight_layout()\npl.show()"
      ]
    }
  ],
  "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.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}