{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nTwo dimension dam break study case.\n===================================\n\nThe dam break is known to be difficult to solve numericaly due to its initial\ndiscontinuity.\n\nIt simulate the behavior of a column of fluid confined by wall. At the initial\ntime, the wall is removed, leading to a wave propagation.\n\nThe shallow water equation will modelise the fluid dynamic and is described\nby the equation below :\n\n\\begin{align}\\begin{cases}\n        \\frac{\\partial h}{\\partial t} + \\frac{\\partial h}{\\partial x} + \\frac{\\partial h}{\\partial y} &= 0 \\\\\n        \\frac{\\partial u}{\\partial t} + u\\,\\frac{\\partial u}{\\partial x} + v\\,\\frac{\\partial u}{\\partial y} &= -g\\,\\frac{\\partial h}{\\partial x}\\\\\n        \\frac{\\partial v}{\\partial t} + u\\,\\frac{\\partial v}{\\partial x} + v\\,\\frac{\\partial v}{\\partial y} &= -g\\,\\frac{\\partial h}{\\partial y}\n    \\end{cases}\\end{align}\n\nThis example will focus on 2D dam break, with wall at the boundaries (simulated\nas no-flux boundary) and a curved dam.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import skfdiff as sf\nimport numpy as np\nfrom scipy.ndimage import gaussian_filter\nimport pylab as plt"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Model definition\n----------------\n\nThe PDE has to be written. It is just the RHS of the model as defined earlier\nwritten as $\\frac{\\partial U}{\\partial t} = F(U)$.\nMoreover, in that specific case, we will also specify that the advective\nterms has to be discretized with an upwind scheme to mitigate the numerical\nerrors that will occur at the discontinuity. This is especially needed as\nwe do not have diffusion at all in the model that could smooth the numerical\nerrors.\n\nBoundary conditions\n+++++++++++++++++++\n\nThe boundary condition will be empty, and will fallback to scikit-fdiff\ndefault (no-flux condition at every boundaries).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "shallow_2D = sf.Model(\n    [\n        \"-dx(h * u) - dy(h * v)\",\n        \"-upwind(u, u, x, 2) - upwind(v, u, y, 2) - g * dxh\",\n        \"-upwind(u, v, x, 2) - upwind(v, v, y, 2) - g * dyh\",\n    ],\n    [\"h(x, y)\", \"u(x, y)\", \"v(x, y)\"],\n    parameters=[\"g\"],\n)"
      ]
    },
    {
      "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\", \"y\"), gaussian_filter(simul.fields[\"h\"], 0.46)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initial condition\n-----------------\n\nA quarter of a circle has been chosen as initial condition, with the fluid\ninside the bottom left of the domain being higher that the rest of the\ndomain.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = y = np.linspace(0, 10, 128)\nxx, yy = np.meshgrid(x, y, indexing=\"ij\")\nu = v = np.zeros_like(xx)\nh = np.where(xx ** 2 + yy ** 2 < 5 ** 2, 2, 1)\nfields = shallow_2D.Fields(x=x, y=y, u=u, v=v, h=h, g=9.81)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Simulation configuration\n------------------------\n\nThe only subtility is that the time-stepping will be desactivated. The\nproblem  does not need to have a variable time-step, and a properly chosen\none will garanty us a good resolution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "simul = shallow_2D.init_simulation(fields, 0.01, tmax=1, time_stepping=False)\nsimul.add_post_process(\"filter\", filter_instabilities)\ncontainer_shallow = simul.attach_container()\nt_max, fields_final = simul.run()\n\ncontainer_shallow.data.h[::20].plot(col=\"t\", col_wrap=3)\nplt.savefig(\"../docs/_static/2D_dam_break.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "![](../../_static/2D_dam_break.png)\n\n\n\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"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}