{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a65b7df9",
   "metadata": {},
   "source": [
    "# TP transport\n",
    "Dans ce TP nous nous intéressons à l'équation de transport\n",
    "$$\n",
    "\\partial_t q = -\\mathsf{div}(\\mathbf{U}\\,q)\n",
    "$$\n",
    "écrite sous forme flux, où $q$ est la concentration de tracer et $\\mathbf{U}$ est la vitesse transportante. La méthode numérique utilise un schéma en temps et une discrétisation du flux $\\mathbf{U}q$.\n",
    "\n",
    "Le but de ce TP est de mieux comprendre\n",
    " - comment les choix numériques influencent une solution\n",
    " - ce qu’est un schéma dispersif vs un schéma dissipatif\n",
    " - le rôle de la déformation (strain)\n",
    " - comment se manifeste une instabilité numérique\n",
    " - ce qu’est une formulation en flux\n",
    " - les différentes façons de calculer un flux\n",
    " - comment on garantit la conservation d’un traceur\n",
    " - la différence entre les termes dissipation et mélange"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3dbae08f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import fluids2d as f2d\n",
    "from tracer_advection import *\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bbc64ff",
   "metadata": {},
   "source": [
    "Commençons par un cas simple: la rotation solide dans un domaine carré"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e18fc39b",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "\n",
    "param.integrator = \"rk3\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "q0 = model.state.q.copy()\n",
    "xy = model.mesh.xy()\n",
    "\n",
    "plt.figure()\n",
    "plt.pcolor(*xy, q0)\n",
    "plt.colorbar()\n",
    "plt.title(\"Distribution initiale de traceur\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df37ca48",
   "metadata": {},
   "source": [
    "On peut à présenter lancer l'intégration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad972087",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ffa7598",
   "metadata": {},
   "source": [
    "L'évolution du traceur au cours du temps doit vous sembler raisonable. Raisonable mais pas tout à fait conforme à une vraie rotation solide. C'est à cause de la forme du domaine. Le traceur devrait tourner de manière uniforme, sans se déformer. Ici il se déforme dans les coins. Pour s'en sortir on va refaire l'expérience dans un domaine circulaire afin qu'il n'y ait **aucune déformation** dans le champ de vitesse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23a6f17d",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "\n",
    "param.integrator = \"rk3\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae3f05ac",
   "metadata": {},
   "source": [
    "Cette fois nous avons bien une rotation solide. La solution semble propre, sauf près des bords où le traceur semble \"baver\". La solution est propre car par défaut le code utilise le meilleur combo numérique: une discrétisation WENO5z pour les flux et un schéma RK3 en temps. Près des bords, comme il y a moins de points disponibles pour faire l'interpolation, les flux sont calculés avec une discrétisation WENO3z et upwind 1er ordre. C'est ce qui est responsable du moins bon comportement.\n",
    "\n",
    "Refaisons l'expérience avec cette fois\n",
    "  - un schéma Leap-Frog en temps\n",
    "  - des flux centrés d'ordre 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7f95057c",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 2\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"centered\"\n",
    "param.integrator = \"LFRA\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94f8788e",
   "metadata": {},
   "source": [
    "La solution est maintenant **bruitée**. Ce bruit est dû au caractère **dispersif** de la discrétisation centrée. Ce bruit est un artefact numérique. Il n'a pas lieu d'être. Non seulement la solution est bruitée mais les valeurs du traceurs qui devraient rester comprises dans l'intervale $[0,1]$ sont maintenant en dehors de cet intervale"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "809ee798",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(model.state.q.ravel(),np.linspace(-1,2,61),density=True)\n",
    "plt.xlabel(\"q\")\n",
    "plt.ylabel(\"P(q)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b912f65e",
   "metadata": {},
   "source": [
    "Pour éliminer ce bruit on peut\n",
    "   1) ajouter un terme de diffusion explicite\n",
    "   2) utiliser une discrétisation amont (upwind), donc décentrée\n",
    "\n",
    "C'est cette 2e approche que nous tester avec le combo: schéma Euler avant + upwind 1er ordre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12c4034b",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 2\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"upwind\"\n",
    "param.integrator = \"ef\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d108821",
   "metadata": {},
   "source": [
    "Il n'y a plus de bruit ... mais il y a maintenant trop de diffusion. C'est le problème avec ce combo. L'ordre est trop bas. Passons à l'ordre 3 en espace."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "16c1d703",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 3\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"upwind\"\n",
    "param.integrator = \"ef\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d7c73610",
   "metadata": {},
   "source": [
    "Cette fois ce combo est **instable**, quelque soit le nombre de courant (`param.cfl`). Il nous faut revenir au schéma RK3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a5e9ee87",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 3\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"upwind\"\n",
    "param.integrator = \"rk3\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60e39503",
   "metadata": {},
   "source": [
    "Cette fois c'est bien stable, le champ semble lisse et sans bruit. En y regardant de plus près, il y a quand même des traces de dispersion. Les valeurs du traceur dépassent à nouveau les bornes. C'est moins fort que pour un schéma centré mais c'est quand même présent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5f77409",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.hist(model.state.q.ravel(),np.linspace(-1,2,61),density=True)\n",
    "plt.xlabel(\"q\")\n",
    "plt.ylabel(\"P(q)\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bfefc36f",
   "metadata": {},
   "source": [
    "Pour éliminer ce comportement il nous faut prendre un **schéma nonlinéaire**. Avec `Fluids2d` on l'obtient en utilisant des schémas WENO pour les flux (combo par défaut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6dccd152",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 6\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"weno\"\n",
    "param.integrator = \"rk3\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"bodyrotation\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c8621fc",
   "metadata": {},
   "source": [
    "Attention, dans les vrais écoulements géophysiques, il y a toujours de la déformation. La déformation est responsable de l'étirement du traceur, ce qui crée de la filamentation. Ces filaments ne peuvent pas être de plus petite échelle que la maille. L'écoulement produit donc naturellement de la petite échelle ! C'est le phénomène de cascade directe pour la variance de traceur.\n",
    "\n",
    "Regardons ce phénomène en prenant l'écoulement dû à un vortex: la vorticité est localisée au centre, nulle ailleurs (au lieu d'être uniforme dans le cas de la rotation solide)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09db8edc",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 6\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"weno\"\n",
    "param.integrator = \"rk3\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"vortex\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59725486",
   "metadata": {},
   "source": [
    "La diffusion induite par le numérique correspond donc à un phénomène physique souhaitable.\n",
    "\n",
    "Pour finir, regardons ce qu'un schéma centré se comporte dans ce cas."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db9f27af",
   "metadata": {},
   "outputs": [],
   "source": [
    "param = f2d.Param()\n",
    "\n",
    "param.model = \"advection\"\n",
    "param.nx = 100\n",
    "param.ny = 100\n",
    "param.tend = 100\n",
    "param.maxite = 5_000\n",
    "param.cfl = 0.9\n",
    "param.nplot = 5\n",
    "param.animation = True\n",
    "param.plotvar = \"q\"\n",
    "param.maxorder = 2\n",
    "param.RAgamma = 0\n",
    "param.compflux = \"centered\"\n",
    "param.integrator = \"LFRA\"\n",
    "\n",
    "model = f2d.Model(param)\n",
    "x, y = model.mesh.xy()\n",
    "r2 = (x-0.5)**2+(y-0.5)**2\n",
    "model.mesh.msk[r2 > 0.5**2] = 0\n",
    "\n",
    "model.mesh.finalize()\n",
    "\n",
    "set_initial_velocity(model, flow=\"vortex\")\n",
    "\n",
    "set_initial_tracer(model)\n",
    "\n",
    "model.run()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b617bac",
   "metadata": {},
   "source": [
    "Les filaments ne sont pas dissipés. A la place, le schéma fabrique du bruit. **Propriété extraordinaire** ce bruit peut être complètement défait! Renversons le champ de vitesse et intégrons à partir de cette situation bruitée."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f734d99d",
   "metadata": {},
   "outputs": [],
   "source": [
    "fliptime(model)\n",
    "model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b4ddde4",
   "metadata": {},
   "source": [
    "Le bruit a disparu! Ce combo est réversible. La variance de traceur est conservée. Sur le papier ça peut sembler une propriété désirable. En pratique, on voit qu'**il faut du mélange**. On peut avoir ce mélange en ajoutant un terme de diffusion ou ... en utilisant des flux upwind (amont) qui embarquent avec eux la diffusion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55b01c56",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fluids2d",
   "language": "python",
   "name": "fluids2d"
  },
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}