{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nOne dimension shallow water, dam break case.\n============================================\n\nThis validation example use the shallow water equation to\nmodel a dam sudden break. The two part of the domain have different\nfluid depth and are separated by a well. A a certain instant, the wall\ndisappear, leading to a discontinuity wave in direction of the lower depth,\nand a rarefaction wave in direction of the higher depth.\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\nand the results are validated on the Randall J. LeVeque book (LeVeque, R. (2002).\nFinite Volume Methods for Hyperbolic Problems (Cambridge Texts in Applied Mathematics).\nCambridge: Cambridge University Press. doi:10.1017/CBO9780511791253).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import pylab as pl\nimport pandas as pd\nfrom skfdiff import Model, Simulation\nimport numpy as np\nfrom scipy.signal import savgol_filter\n\nshallow_water = Model([\"-dx(h * u)\", \"-upwind(u, u, x, 1) - dxh\"], [\"h(x)\", \"u(x)\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Filter post-process\n-------------------\n\nAs the discontinuity will be still harsh to handle by a \"simple\" finite\ndifference solver (a finite volume solver is usually the tool you will need),\nwe will use a gaussian filter that will be applied to the fluid height. This\nwill smooth the oscillation (generated by numerical errors). This can be seen\nas a way to introduce numerical diffusion to the system, often done by adding\na diffusion term in the model. The filter has to be carefully tuned (the same\nway an artificial diffusion has its coefficient diffusion carefully chosen)\nto smooth the numerical oscillation without affecting the physical behavior\nof the simulation.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def filter_instabilities(simul):\n    simul.fields[\"h\"] = (\"x\",), savgol_filter(simul.fields[\"h\"], 21, 4)\n\n\nx, dx = np.linspace(-5, 5, 1000, retstep=True)\nh = np.where(x < 0, 3, 1)\nu = x * 0\n\ninit_fields = shallow_water.Fields(x=x, h=h, u=u)\n\nsimul = Simulation(\n    shallow_water,\n    t=0,\n    dt=0.02,\n    tmax=2,\n    fields=init_fields,\n    time_stepping=False,\n    id=\"dambreak\",\n)\n\nsimul.add_post_process(\"filter\", filter_instabilities)\n\ncontainer = simul.attach_container()\n\nsimul.run()\n\ndata = container.data.sel(t=[0, 0.5, 2], method=\"nearest\")\nfig, axs = pl.subplots(2, 2, sharex=\"all\", figsize=(5.5, 2.5))\nfor i, t in enumerate(data.t):\n    if i == 1:\n        continue\n    if i == 2:\n        pl.sca(axs[i - 1, 0])\n    else:\n        pl.sca(axs[i, 0])\n    data.sel(t=t).h.plot(color=\"black\", label=\"Sol.\")\n    ref_data = pd.read_csv(\"valid_randall/dam_h%i.csv\" % i)\n    pl.scatter(ref_data.x, ref_data.h, color=\"red\", marker=\".\", label=\"Ref.\")\n    pl.title(\"\")\n    pl.ylabel(r\"$h$\")\n    pl.xlim(-5, 5)\n    pl.ylim(0.5, 4)\n    if i == 0:\n        pl.legend()\n    if i == 2:\n        pl.sca(axs[i - 1, 1])\n    else:\n        pl.sca(axs[i, 1])\n    (data.sel(t=t).h * data.sel(t=t).u).plot(color=\"black\", label=\"Sol.\")\n    ref_data = pd.read_csv(\"valid_randall/dam_hu%i.csv\" % i)\n    pl.scatter(ref_data.x, ref_data.h, color=\"red\", marker=\".\", label=\"Ref.\")\n    pl.title(\"\")\n    pl.ylabel(r\"$u\\,h$\")\n    pl.xlim(-5, 5)\n    pl.ylim(-0.5, 1.5)\n    pl.tight_layout()\n\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
}