CDMATH/tests/doc/2DWaveSystemVF_stationary/Convergence_WaveSystem_Upwind_SQUARE.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Upwind scheme for the Wave System"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Wave System on the square"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider the following Wave system with periodic boundary conditions\n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{array}{l}\n",
    "\\partial_t p + c^2\\nabla\\cdot\\vec q = 0\\\\\n",
    "\\partial_t \\vec q + \\quad\\vec\\nabla p = 0\n",
    "\\end{array}\\right..\n",
    "$$\n",
    "\n",
    "The wave system can be written in matrix form \n",
    "$$\n",
    "\\partial_t\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    " p \\\\\n",
    " \\vec q\n",
    "\\end{array}\\right)\n",
    "+\n",
    "\\left(\\begin{array}{cc}\n",
    " 0   & c^2 \\nabla \\cdot \\\\ \n",
    " \\vec\\nabla & 0\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    " p \\\\ \n",
    " \\vec q\n",
    "\\end{array}\\right)=0\n",
    "$$\n",
    "\n",
    "In $d$ space dimensions the wave system is an hyperbolic system of $d+1$ equations\n",
    "$$\n",
    "\\partial_t U +\\sum_{i=1}^d A_i\\partial_{x_i} U=0,\\quad U={}^t(p,\\vec q)\n",
    "$$\n",
    "where the jacobian matrix is\n",
    "$$\n",
    "A(\\vec n)=\\sum_{i=1}^d n_i A_i =\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    " 0              & c^2 {}^t\\vec n \\\\ \n",
    " \\vec n & 0\n",
    "\\end{array}\\right),\\quad \\vec n\\in\\mathbb{R}^d.\n",
    "$$\n",
    "has $d+1$ eigenvalues $-c,0,\\dots,0,c$.\n",
    "\n",
    "The wave system also takes the conservative form\n",
    "$$\n",
    "\\partial_t U + \\nabla\\cdot F(U)=0,\n",
    "$$\n",
    "where the flux matrix $F$ is defined by\n",
    "$$\n",
    "F(U)\\vec n=A(\\vec n)U,\\quad \\vec n\\in\\mathbb{R}.\n",
    "$$\n",
    "\n",
    "On the square domain $\\Omega= [0,1]\\times [0,1]$ we consider the initial data $(p_0,\\vec q_0)$ \n",
    "\n",
    "$$\n",
    "\\left\\{\\begin{array}{l}\n",
    "p_0(x,y)=constant\\\\\n",
    "q_{x0}(x,y)= \\sin(\\pi x) \\cos(\\pi y)\\\\\n",
    "q_{y0}(x,y)=-\\sin(\\pi y) \\cos(\\pi x)\n",
    "\\end{array}\\right..\n",
    "$$  \n",
    "\n",
    "The initial data $(p_0,q_x,q_y)$ is a stationary solution of the wave system.\n",
    "\n",
    "$$\n",
    "||\\vec q_0||^2=\\sin(\\pi x)^2 \\cos(\\pi y)^2+\\sin(\\pi y)^2 \\cos(\\pi x)^2=\\frac{1}{2}\\sin(\\pi (x + y))^2+\\frac{1}{2}\\sin(\\pi (x-y))^2\n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Upwind scheme for the Wave System"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The domain $\\Omega$ is decomposed into cells $C_i$.\n",
    "\n",
    "- $|C_i|$ is the measure of the cell $C_i$.\n",
    "\n",
    "- $f_{ij}$ is the interface between two cells $C_i$ and $C_j$. \n",
    "\n",
    "- $s_{ij}$ is the measure of the interface $f_{ij}$.\n",
    "\n",
    "- $d_{ij}$ is the distance between the centers of mass of the two cells $C_i$ and $C_j$.\n",
    "\n",
    "The semi-discrete colocated finite volume equation is\n",
    "$$\n",
    "\\partial_t U + \\frac{1}{|C_i|} \\sum s_{ij}F_{ij}=0,\n",
    "$$\n",
    "where\n",
    "$U_i$ is the approximation of $U$ in the cell $C_i$.\n",
    "\n",
    "$F_{ij}$ is a numerical approximation of the outward normal interfacial flux from cell $i$ to cell $j$ usually in the upwind form\n",
    "$$\n",
    "F_{ij}=\\frac{F(U_i)+F(U_j)}{2}\\vec n-D(\\vec n)\\frac{U_i-U_j}{2}.\n",
    "$$\n",
    "In the case of the upwind scheme the upwind matrix is\n",
    "$$\n",
    "D_{upwind}(\\vec n)=|A(\\vec n)| = \n",
    "c\n",
    "\\left(\n",
    "\\begin{array}{cc}\n",
    "1 & 0 \\\\\n",
    "0 & \\vec n\\otimes\\vec n\n",
    "\\end{array}\\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some bibliographical remarks\n",
    "\n",
    "Most of the literature on finite volume methods for hyperbolic equations focus on scalar equation whilst we are interested in systems of equation. However the following references do apply to the wave system and the upwind scheme.\n",
    "\n",
    "- The upwind scheme converges on general meshes with order 1/2  \n",
    "  *J.-P. Vila, P. Villedieu, Convergence de la méthode des volumes finis pour les systèmes de Friedrichs, C. R. Acad. Sci. Paris, t. 325, Série I, 1997.*\n",
    "\n",
    "Stationary regimes are of great importance in several application. In particular in the wave system the stationary momentum should satisfy $\\partial_x q_x+\\partial_y q_y=0$. However at the discrete level, the stationary regime can be very different.\n",
    "\n",
    "- Study of the kernel of the discrete divergence, which corresponds to the discrete stationary regimes of the upwind scheme  \n",
    "  *S. Dellacherie, P. Omnes, F. Rieper, The influence of cell geometry on the Godunov scheme applied to the linear wave equation, J. Comp. Phys., 229(14), pp. 5315-5338, 2010*\n",
    "    - On rectangular meshes the discrete kernel is $\\partial_x q_x=\\partial_y q_y=0$ which is a much smaller space than $\\partial_x q_x+\\partial_y q_y=0$.\n",
    "    - On triangular meshes the discrete kernel is much larger. In practical simulation we observe convergence towards the correct stationary regime.\n",
    "\n",
    "- Relationships between stationarity preservation, low Mach precision and vorticity preservation  \n",
    "  *Wasilij Barsukow: Stationarity preserving schemes for multi-dimensional linear systems, 2017, accepted in Math. Comp*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The script"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "#Condition initiale\n",
    "pressure_field, velocity_field = initial_conditions_wave_system(my_mesh)\n",
    "\n",
    "#Pas de temps\n",
    "dt = cfl * dx_min / c0\n",
    "\n",
    "#Matrice des systèmes linéaires\n",
    "divMat=computeDivergenceMatrix(my_mesh,nbVoisinsMax,dt,test_bc)\n",
    "\n",
    "# Construction du vecteur inconnu\n",
    "Un=cdmath.Vector(nbCells*(dim+1))\n",
    "for k in range(nbCells):\n",
    "    Un[k*(dim+1)+0] =     pressure_field[k]\n",
    "    Un[k*(dim+1)+1] =rho0*velocity_field[k,0]\n",
    "    Un[k*(dim+1)+2] =rho0*velocity_field[k,1]\n",
    "\n",
    "\n",
    "# Création du solveur linéaire\n",
    "LS=cdmath.LinearSolver(divMat,Un,iterGMRESMax, precision, \"GMRES\",\"ILU\")\n",
    "\n",
    "# Time loop\n",
    "while (it<ntmax and time <= tmax and not isStationary):\n",
    "    LS.setSndMember(Un)\n",
    "    Un=LS.solve();\n",
    "    Un.writeVTK\n",
    "    \n",
    "# Automatic postprocessing :  save 2D picture and plot diagonal data\n",
    "#===========================\n",
    "diag_data=VTK_routines.Extract_field_data_over_line_to_numpyArray(my_ResultField,[0,1,0],[1,0,0], resolution)\n",
    "plt.legend()\n",
    "plt.xlabel('Position on diagonal line')\n",
    "plt.ylabel('Value on diagonal line')\n",
    "if len(sys.argv) >1 :\n",
    "    plt.title('Plot over diagonal line for finite volumes \\n for Wave system on a 2D square '+my_mesh.getName())\n",
    "    plt.plot(curv_abs, diag_data, label= str(nbCells)+ ' cells mesh')\n",
    "    plt.savefig(\"FiniteVolumes2D_square_ResultField_\"+str(nbCells)+ '_cells'+\"_PlotOverDiagonalLine.png\")\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical results for upwind scheme on cartesian meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cartesian meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mesh 1 | mesh 2 | mesh 3\n",
    "     - | -    - | -\n",
    "![](2DWaveSystemUpwind/squareWithSquares_2.png) | ![](2DWaveSystemUpwind/squareWithSquares_3.png)  | ![](2DWaveSystemUpwind/squareWithSquares_4.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Velocity initial data (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwind15x15_velocity_initial.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwind31x31_velocity_initial.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwind51x51_velocity_initial.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary velocity (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwind15x15_velocity_Stat.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwind31x31_velocity_Stat.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwind51x51_velocity_Stat.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary velocity (x component)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwind15x15_velocityX_Stat.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwind31x31_velocityX_Stat.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwind51x51_velocityX_Stat.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary velocity (y component)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwind15x15_velocityY_Stat.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwind31x31_velocityY_Stat.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwind51x51_velocityY_Stat.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence on velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](2DWaveSystemUpwind/meshSquareWithSquares_Velocity_2DWaveSystemUpwind_ConvergenceCurve.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical results for upwind scheme on triangular meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Triangular meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mesh 1 | mesh 2 | mesh 3\n",
    "     - | -    - | -\n",
    "![](2DWaveSystemUpwind/squareWithTriangles_2.png) | ![](2DWaveSystemUpwind/squareWithTriangles_3.png)  | ![](2DWaveSystemUpwind/squareWithTriangles_4.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Velocity initial data (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithTriangles224_velocity_initial.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithTriangles934_velocity_initial.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithTriangles6422_velocity_initial.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary velocity  (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithTriangles224_velocity_Stat.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithTriangles934_velocity_Stat.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithTriangles6422_velocity_Stat.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence on stationary velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](2DWaveSystemUpwind/squareWithTriangles_Velocity_2DWaveSystemUpwind_Triangles_ConvergenceCurve.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical results for upwind scheme on checkerboard meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Checkerboard meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mesh 1 | mesh 2 | mesh 3\n",
    "     - | -    - | -\n",
    "![](2DWaveSystemUpwind/checkerboard_17x17.png) | ![](2DWaveSystemUpwind/checkerboard_33x33.png)  | ![](2DWaveSystemUpwind/checkerboard_65x65.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Velocity initial data (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithCheckerboard204_velocity_initial.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithCheckerboard724_velocity_initial.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithCheckerboard2724_velocity_initial.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary velocity (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithCheckerboard204_velocity_Stat.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithCheckerboard724_velocity_Stat.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithCheckerboard2724_velocity_Stat.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence on stationary velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](2DWaveSystemUpwind/squareWithCheckerboard_Velocity_2DWaveSystemUpwind_Checkerboard_ConvergenceCurve.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numerical results for upwind scheme on flat triangle meshes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flat cross triangle meshes (from a $(n,n^2)$ rectangular grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "mesh 1 | mesh 2 | mesh 3\n",
    "     - | -    - | -\n",
    "![](2DWaveSystemUpwind/squareWithCrossTriangles_00.png) | ![](2DWaveSystemUpwind/squareWithFlatTriangles_0.png)  | ![](2DWaveSystemUpwind/squareWithFlatTriangles_1.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Velocity initial data (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithFlatCrossTriangles500_velocity_initial.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithFlatCrossTriangles5324_velocity_initial.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithFlatCrossTriangles37044_velocity_initial.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stationary velocity (magnitude)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "result 1 | result 2 | result 3\n",
    "       - | -      - | -\n",
    "![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithFlatCrossTriangles500_velocity_Stat.png) | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithFlatCrossTriangles5324_velocity_Stat.png)  | ![](2DWaveSystemUpwind/WaveSystem2DUpwindsquareWithFlatCrossTriangles37044_velocity_Stat.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convergence on stationary velocity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![](2DWaveSystemUpwind/squareWithFlatCrossTriangles_Velocity_2DWaveSystemUpwind_ConvergenceCurve.png)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}