{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\nTwo dimension acoustic waves simulation.\n========================================\n\nThank to Matteo Ravasi (@mrava87) for this example.\n\nThis example show the propagation of a 2D acoustic wave across\ntwo different medium with different sound speed / density.\n\nThe wave start as a source as a signal in time at a given location\nin space. This source is inserted via a hook function.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\n\nfrom scipy.ndimage import gaussian_filter\nfrom skfdiff import Model, Simulation\n\nfrom collections import defaultdict\n\n\ndef hook(t, fields):\n    f0 = 40\n    t0 = 0.01\n    q = np.zeros_like(fields.P)\n    q[20, ny // 2] = (1 - 2 * (np.pi * f0 * (t - t0)) ** 2) * np.exp(\n        -(np.pi * f0 * (t - t0)) ** 2\n    )\n    fields[\"q\"] = (\"x\", \"y\"), gaussian_filter(q, sigma=2)\n    return fields\n\n\nbc = defaultdict(lambda: (\"dirichlet\", \"dirichlet\"))\n\nmodel = Model(\n    [\n        \"-upwind(1/rho, P, x, 1)\",\n        \"-upwind(1/rho, P, y, 1)\",\n        \"upwind(-1/k, Vx, x, 1) + upwind(-1/k, Vy, y, 1) + q\",\n    ],\n    [\"Vx(x, y)\", \"Vy(x, y)\", \"P(x, y)\"],\n    parameters=[\"rho(x, y)\", \"k(x, y)\", \"q(x, y)\"],\n    boundary_conditions=bc,\n)\n\n\nnx = ny = 256\nx = np.linspace(-200, 200, nx)\ny = np.linspace(-200, 200, ny)\nxx, yy = np.meshgrid(x, y, indexing=\"ij\")\n\ndt = 1e-3\ntmax = 0.15\nU = np.zeros((nx, ny))\nq = np.zeros((nx, ny))\nc = 1800 * np.ones((nx, ny))\nc[nx // 3 :] = 3000\nrho = 1000 * np.ones((nx, ny))\nrho[nx // 3 :] = 4000\nk = 1 / (c ** 2 * rho)\n\n# We fill the fields container\ninitial_fields = model.Fields(x=x, y=y, P=U, Vx=U, Vy=U, q=q, k=k, rho=rho)\n\n# We initialize the simulation\nsimulation = Simulation(\n    model,\n    initial_fields,\n    dt=dt,\n    tol=1e-1,\n    tmax=tmax,\n    hook=hook,\n    solver=\"iteratif\",\n    time_stepping=True,\n)\n\ncontainer_acoustic = simulation.attach_container()\nsimulation.run()\n\nfacet = container_acoustic.data.P[::30].T.plot(\n    col=\"t\", col_wrap=3, cmap=\"seismic\", vmin=-3e-5, vmax=3e-5\n)\nfor ax in facet.axes.flatten():\n    left_domain = plt.Polygon(\n        [\n            (x.min(), y.min()),\n            (x[nx // 3], y.min()),\n            (x[nx // 3], y.max()),\n            (x.min(), y.max()),\n        ],\n        color=\"blue\",\n        alpha=0.2,\n    )\n    right_domain = plt.Polygon(\n        [\n            (x[nx // 3], y.min()),\n            (x.max(), y.min()),\n            (x.max(), y.max()),\n            (x[nx // 3], y.max()),\n        ],\n        color=\"red\",\n        alpha=0.2,\n    )\n    ax.add_patch(left_domain)\n    ax.add_patch(right_domain)\nplt.savefig(\"../docs/_static/2D_acoustic.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "![](../../_static/2D_acoustic.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
}