{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Rotating double gyre\n",
    "\n",
    "## Description\n",
    "\n",
    "The rotating double gyre model was introduced by\n",
    "[Mosovsky & Meiss](https://doi.org/10.1137/100794110). It can be derived from\n",
    "the stream function\n",
    "$$\n",
    "\\psi(x,y,t)=(1−s(t))\\psi_P +s(t)\\psi_F \\\\ \\psi_P (x, y) = \\sin(2\\pi x) \\sin(\\pi y) \\\\ \\psi_F (x, y) = \\sin(\\pi x) \\sin(2\\pi y)\n",
    "$$\n",
    "where $s$ is (usually taken to be) a cubic interpolating function satisfying\n",
    "$s(0) = 0$ and $s(1) = 1$. It therefore interpolates two double-gyre-type\n",
    "velocity fields, from horizontally to vertically arranged counter-rotating gyres.\n",
    "The corresponding velocity field can be conveniently defined using the\n",
    "`@velo_from_stream` macro from [`StreamMacros.jl`](https://github.com/CoherentStructures/StreamMacros.jl):"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Distributed\n",
    "nprocs() == 1 && addprocs()\n",
    "\n",
    "@everywhere using StreamMacros\n",
    "const rot_double_gyre = @velo_from_stream Ψ_rot_dgyre begin\n",
    "    st          = heaviside(t)*heaviside(1-t)*t^2*(3-2*t) + heaviside(t-1)\n",
    "    heaviside(x)= 0.5*(sign(x) + 1)\n",
    "    Ψ_P         = sin(2π*x)*sin(π*y)\n",
    "    Ψ_F         = sin(π*x)*sin(2π*y)\n",
    "    Ψ_rot_dgyre = (1-st) * Ψ_P + st * Ψ_F\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "![](https://raw.githubusercontent.com/natschil/misc/db22aeef/images/double_gyre.gif)\n",
    "\n",
    "## FEM-Based Methods\n",
    "\n",
    "The following code demonstrates how to use these methods."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using CoherentStructures\n",
    "LL, UR = (0.0, 0.0), (1.0, 1.0)\n",
    "ctx, _ = regularTriangularGrid((50, 50), LL, UR)\n",
    "\n",
    "A = x -> mean_diff_tensor(rot_double_gyre, x, [0.0, 1.0], 1.e-10, tolerance= 1.e-4)\n",
    "K = assembleStiffnessMatrix(ctx, A)\n",
    "M = assembleMassMatrix(ctx)\n",
    "λ, v = CoherentStructures.get_smallest_eigenpairs(-K, M, 6);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This velocity field is given by the `rot_double_gyre` function. The third\n",
    "argument to `mean_diff_tensor` is a vector of time instances at which we evaluate\n",
    "(and subsequently average) the pullback diffusion tensors. The fourth parameter\n",
    "is the step size δ used for the finite-difference scheme, `tolerance` is passed\n",
    "to the ODE solver from [DifferentialEquations.jl](http://juliadiffeq.org/). In\n",
    "the above, `A(x)` approximates the mean diffusion tensor given by\n",
    "\n",
    "$$\n",
    "A(x) = \\sum_{t \\in \\mathcal T}(D\\Phi^t(x))^{-1} (D\\Phi^t x)^{-T}\n",
    "$$\n",
    "\n",
    "The eigenfunctions saved in `v` approximate those of $\\Delta^{dyn}$"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "import Plots\n",
    "res = [plot_u(ctx, v[:,i], 100, 100, colorbar=:none, clim=(-3,3)) for i in 1:6];\n",
    "fig = Plots.plot(res..., margin=-10Plots.px)\n",
    "\n",
    "Plots.plot(fig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Looking at the spectrum, there appears a gap after the third eigenvalue."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "spectrum_fig = Plots.scatter(1:6, real.(λ))\n",
    "\n",
    "Plots.plot(spectrum_fig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can use the [Clustering.jl](https://github.com/JuliaStats/Clustering.jl) package\n",
    "to compute coherent structures from the first two nontrivial eigenfunctions:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Clustering\n",
    "\n",
    "ctx2, _ = regularTriangularGrid((200, 200))\n",
    "v_upsampled = sample_to(v, ctx, ctx2)\n",
    "\n",
    "numclusters=2\n",
    "res = kmeans(permutedims(v_upsampled[:,2:numclusters+1]), numclusters + 1)\n",
    "u = kmeansresult2LCS(res)\n",
    "res = Plots.plot([plot_u(ctx2, u[:,i], 200, 200, c=:viridis, cbar=:none) for i in 1:3]...)\n",
    "\n",
    "Plots.plot(res)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Geodesic vortices"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we demonstrate how to calculate black-hole vortices, see\n",
    "Geodesic elliptic material vortices for references and details."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@everywhere using CoherentStructures, OrdinaryDiffEq\n",
    "q = 21\n",
    "const tspan = range(0., stop=1., length=q)\n",
    "nx = ny = 101\n",
    "xmin, xmax, ymin, ymax = 0.0, 1.0, 0.0, 1.0\n",
    "xspan = range(xmin, stop=xmax, length=nx)\n",
    "yspan = range(ymin, stop=ymax, length=ny)\n",
    "P = tuple.(xspan, permutedims(yspan))\n",
    "const δ = 1.e-6\n",
    "mcg_tensor = u -> av_weighted_CG_tensor(rot_double_gyre, u, tspan, δ; tolerance=1e-6, solver=Tsit5())\n",
    "\n",
    "C̅ = pmap(mcg_tensor, P; batch_size=ceil(Int, length(P)/nprocs()^2))\n",
    "p = LCSParameters(0.5)\n",
    "vortices, singularities = ellipticLCS(C̅, xspan, yspan, p; outermost=true)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The results are then visualized as follows."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "trace = tensor_invariants(C̅)[5]\n",
    "fig = plot_vortices(vortices, singularities, (xmin, ymin), (xmax, ymax);\n",
    "    bg=trace, xspan=xspan, yspan=yspan, title=\"DBS field and transport barriers\", showlabel=true, clims=(0,5))\n",
    "Plots.plot(fig)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "---\n",
    "\n",
    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.5"
  },
  "kernelspec": {
   "name": "julia-1.8",
   "display_name": "Julia 1.8.5",
   "language": "julia"
  }
 },
 "nbformat": 4
}