{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 15: Structured mesh with curvilinear mapping"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note:** To improve responsiveness via caching, the notebooks are updated only once a week. They are only\n",
    "available for the latest stable release of Trixi.jl at the time of caching."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we want to introduce another mesh type of [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).\n",
    "More precisely, this tutorial is about the curved mesh type `StructuredMesh` supporting\n",
    "curved meshes."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Creating a curved mesh\n",
    "There are two basic options to define a curved `StructuredMesh` in Trixi.jl. You can\n",
    "implement curves for the domain boundaries, or alternatively, set up directly the complete\n",
    "transformation mapping. We now present one short example each."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Mesh defined by domain boundary curves\n",
    "Both examples are based on a semdiscretization of the 2D compressible Euler equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "using Trixi\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We start with a pressure perturbation at `(xs, 0.0)` as initial condition."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function initial_condition_pressure_perturbation(x, t, equations::CompressibleEulerEquations2D)\n",
    "  xs = 1.5 # location of the initial disturbance on the x axis\n",
    "  w = 1/8 # half width\n",
    "  p = exp(-log(2) * ((x[1]-xs)^2 + x[2]^2)/w^2) + 1.0\n",
    "  v1 = 0.0\n",
    "  v2 = 0.0\n",
    "  rho = 1.0\n",
    "\n",
    "  return prim2cons(SVector(rho, v1, v2, p), equations)\n",
    "end\n",
    "initial_condition = initial_condition_pressure_perturbation"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Initialize every boundary as a `boundary_condition_slip_wall`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "boundary_conditions = boundary_condition_slip_wall"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The approximation setup is an entropy-stable split-form DG method with `polydeg=4`. We are using\n",
    "the two fluxes `flux_ranocha` and `flux_lax_friedrichs`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(polydeg=4, surface_flux=flux_lax_friedrichs,\n",
    "               volume_integral=VolumeIntegralFluxDifferencing(flux_ranocha))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We want to define a circular cylinder as physical domain. It contains an inner semicircle with\n",
    "radius `r0` and an outer semicircle of radius `r1`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![](https://user-images.githubusercontent.com/74359358/159492083-1709510f-8ba4-4416-9fb1-e2ed2a11c62c.png)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The domain boundary curves with curve parameter in $[-1,1]$ are sorted as shown in the sketch.\n",
    "They always are orientated from negative to positive coordinate, such that the corners have to\n",
    "fit like this $f_1(+1) = f_4(-1)$, $f_3(+1) = f_2(-1)$, etc."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In our case we can define the domain boundary curves as follows:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "r0 = 0.5 # inner radius\n",
    "r1 = 5.0 # outer radius\n",
    "f1(xi)  = SVector( r0 + 0.5 * (r1 - r0) * (xi + 1), 0.0) # right line\n",
    "f2(xi)  = SVector(-r0 - 0.5 * (r1 - r0) * (xi + 1), 0.0) # left line\n",
    "f3(eta) = SVector(r0 * cos(0.5 * pi * (eta + 1)), r0 * sin(0.5 * pi * (eta + 1))) # inner circle\n",
    "f4(eta) = SVector(r1 * cos(0.5 * pi * (eta + 1)), r1 * sin(0.5 * pi * (eta + 1))) # outer circle"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We create a curved mesh with 16 x 16 elements. The defined domain boundary curves are passed as a tuple."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "cells_per_dimension = (16, 16)\n",
    "mesh = StructuredMesh(cells_per_dimension, (f1, f2, f3, f4), periodicity=false)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, we define the simulation with endtime `T=3` with `semi`, `ode` and `callbacks`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,\n",
    "                                    boundary_conditions=boundary_conditions)\n",
    "\n",
    "tspan = (0.0, 3.0)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "analysis_interval = 100\n",
    "analysis_callback = AnalysisCallback(semi, interval=analysis_interval)\n",
    "\n",
    "alive_callback = AliveCallback(analysis_interval=analysis_interval)\n",
    "\n",
    "stepsize_callback = StepsizeCallback(cfl=0.9)\n",
    "\n",
    "callbacks = CallbackSet(analysis_callback,\n",
    "                        alive_callback,\n",
    "                        stepsize_callback);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Running the simulation"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),\n",
    "            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback\n",
    "            save_everystep=false, callback=callbacks);\n",
    "\n",
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "pd = PlotData2D(sol)\n",
    "plot(pd[\"p\"])\n",
    "plot!(getmesh(pd))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Mesh directly defined by the transformation mapping\n",
    "As mentioned before, you can also define the domain for a `StructuredMesh` by directly setting up\n",
    "a transformation mapping. Here, we want to present a nice mapping, which is often used to test\n",
    "free-stream preservation. Exact free-stream preservation is a crucial property of any numerical\n",
    "method on curvilinear grids. The mapping is a reduced 2D version of the mapping described in\n",
    "[Rueda-Ramírez et al. (2021), p.18](https://arxiv.org/abs/2012.12040)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "using Trixi\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As mentioned, this mapping is used for testing free-stream preservation. So, we use a constant\n",
    "initial condition."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "initial_condition = initial_condition_constant\n",
    "\n",
    "solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We define the transformation mapping with variables in $[-1, 1]$ as described in\n",
    "Rueda-Ramírez et al. (2021), p.18 (reduced to 2D):"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function mapping(xi_, eta_)\n",
    "  # Transform input variables between -1 and 1 onto [0,3]\n",
    "  xi = 1.5 * xi_ + 1.5\n",
    "  eta = 1.5 * eta_ + 1.5\n",
    "\n",
    "  y = eta + 3/8 * (cos(1.5 * pi * (2 * xi - 3)/3) *\n",
    "                   cos(0.5 * pi * (2 * eta - 3)/3))\n",
    "\n",
    "  x = xi + 3/8 * (cos(0.5 * pi * (2 * xi - 3)/3) *\n",
    "                  cos(2 * pi * (2 * y - 3)/3))\n",
    "\n",
    "  return SVector(x, y)\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Instead of a tuple of boundary functions, the `mesh` now has the mapping as its parameter."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "cells_per_dimension = (16, 16)\n",
    "mesh = StructuredMesh(cells_per_dimension, mapping)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)\n",
    "\n",
    "tspan = (0.0, 1.0)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "analysis_callback = AnalysisCallback(semi, interval=250)\n",
    "\n",
    "stepsize_callback = StepsizeCallback(cfl=0.8)\n",
    "\n",
    "callbacks = CallbackSet(analysis_callback,\n",
    "                        stepsize_callback)\n",
    "\n",
    "sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),\n",
    "            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback\n",
    "            save_everystep=false, callback=callbacks);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we want to verify the free-stream preservation property and plot the mesh. For the verification,\n",
    "we calculate the absolute difference of the first conservation variable density `u[1]` and `1.0`.\n",
    "To plot this error and the mesh, we are using the visualization feature `ScalarPlotData2D`,\n",
    "explained in visualization."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "error_density = let u = Trixi.wrap_array(sol.u[end], semi)\n",
    "  abs.(u[1, :, :, :] .- 1.0) # density, x, y, elements\n",
    "end\n",
    "pd = ScalarPlotData2D(error_density, semi)\n",
    "\n",
    "using Plots\n",
    "plot(pd, title=\"Error in density\")\n",
    "plot!(getmesh(pd))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We observe that the errors in the variable `density` are at the level of machine accuracy.\n",
    "Moreover, the plot shows the mesh structure resulting from our transformation mapping."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Of course, you can also use other mappings as for instance shifts by $(x, y)$\n",
    "```julia\n",
    "mapping(xi, eta) = SVector(xi + x, eta + y)\n",
    "```\n",
    "or rotations with a rotation matrix $T$\n",
    "```julia\n",
    "mapping(xi, eta) = T * SVector(xi, eta).\n",
    "```"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For more curved mesh mappings, please have a look at some\n",
    "[elixirs for `StructuredMesh`](https://github.com/trixi-framework/Trixi.jl/tree/main/examples).\n",
    "For another curved mesh type, there is a tutorial about Trixi.jl's\n",
    "unstructured mesh type [`UnstructuredMesh2D`] and its use of the\n",
    "[High-Order Hex-Quad Mesh (HOHQMesh) generator](https://github.com/trixi-framework/HOHQMesh),\n",
    "created and developed by David Kopriva."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Package versions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "These results were obtained using the following versions."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using InteractiveUtils\n",
    "versioninfo()\n",
    "\n",
    "using Pkg\n",
    "Pkg.status([\"Trixi\", \"OrdinaryDiffEq\", \"Plots\"],\n",
    "           mode=PKGMODE_MANIFEST)"
   ],
   "metadata": {},
   "execution_count": null
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}