{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 7: DG schemes via `DGMulti` solver"
   ],
   "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": [
    "`DGMulti` is a DG solver that allows meshes with simplex elements. The basic idea and\n",
    "implementation of this solver is explained in section \"Meshes\".\n",
    "Here, we want to give some examples and a quick overview about the options with `DGMulti`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We start with a simple example we already used in the tutorial about flux differencing.\n",
    "There, we implemented a simulation with `initial_condition_weak_blast_wave` for the\n",
    "2D compressible Euler equations `CompressibleEulerEquations2D` and used the DG formulation\n",
    "with flux differencing using volume flux `flux_ranocha` and surface flux `flux_lax_friedrichs`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we want to implement the equivalent example, only now using the `DGMulti` solver\n",
    "instead of `DGSEM`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "\n",
    "initial_condition = initial_condition_weak_blast_wave"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To use the Gauss-Lobatto nodes again, we choose `approximation_type=SBP()` in the solver.\n",
    "Since we want to start with a Cartesian domain discretized with squares, we use the element\n",
    "type `Quad()`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dg = DGMulti(polydeg = 3,\n",
    "             element_type = Quad(),\n",
    "             approximation_type = SBP(),\n",
    "             surface_flux = flux_lax_friedrichs,\n",
    "             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))\n",
    "\n",
    "cells_per_dimension = (32, 32)\n",
    "mesh = DGMultiMesh(dg,\n",
    "                   cells_per_dimension, # initial_refinement_level = 5\n",
    "                   coordinates_min=(-2.0, -2.0),\n",
    "                   coordinates_max=( 2.0,  2.0),\n",
    "                   periodicity=true)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,\n",
    "                                    boundary_conditions=boundary_condition_periodic)\n",
    "tspan = (0.0, 0.4)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "alive_callback = AliveCallback(alive_interval=10)\n",
    "analysis_callback = AnalysisCallback(semi, interval=100, uEltype=real(dg))\n",
    "callbacks = CallbackSet(analysis_callback, alive_callback);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Run the simulation with the same time integration algorithm as before."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = solve(ode, RDPK3SpFSAL49(), abstol=1.0e-6, reltol=1.0e-6,\n",
    "            callback=callbacks, save_everystep=false);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "pd = PlotData2D(sol)\n",
    "plot(pd)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plot(pd[\"rho\"])\n",
    "plot!(getmesh(pd))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This simulation is not as fast as the equivalent with `TreeMesh` since no special optimizations for\n",
    "quads (for instance tensor product structure) have been implemented. Figure 4 in [\"Efficient implementation of modern entropy stable\n",
    "and kinetic energy preserving discontinuous Galerkin methods for conservation laws\"](https://arxiv.org/abs/2112.10517)\n",
    "(2021) provides a nice runtime comparison between the different mesh types. On the other hand,\n",
    "the functions are more general and thus we have more option we can choose from."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Simulation with Gauss nodes\n",
    "For instance, we can change the approximation type of our simulation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "initial_condition = initial_condition_weak_blast_wave"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now use Gauss nodes instead of Gauss-Lobatto nodes which can be done for the element types\n",
    "`Quad()` and `Hex()`. Therefore, we set `approximation_type=GaussSBP()`. Alternatively, we\n",
    "can use a modal approach using the approximation type `Polynomial()`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dg = DGMulti(polydeg = 3,\n",
    "             element_type = Quad(),\n",
    "             approximation_type = GaussSBP(),\n",
    "             surface_flux = flux_lax_friedrichs,\n",
    "             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))\n",
    "\n",
    "cells_per_dimension = (32, 32)\n",
    "mesh = DGMultiMesh(dg,\n",
    "             cells_per_dimension, # initial_refinement_level = 5\n",
    "             coordinates_min=(-2.0, -2.0),\n",
    "             coordinates_max=( 2.0,  2.0),\n",
    "             periodicity=true)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,\n",
    "                              boundary_conditions=boundary_condition_periodic)\n",
    "tspan = (0.0, 0.4)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "alive_callback = AliveCallback(alive_interval=10)\n",
    "analysis_callback = AnalysisCallback(semi, interval=100, uEltype=real(dg))\n",
    "callbacks = CallbackSet(analysis_callback, alive_callback);\n",
    "\n",
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,\n",
    "            ode_default_options()..., callback=callbacks);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "pd = PlotData2D(sol)\n",
    "plot(pd)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Simulation with triangular elements\n",
    "Also, we can set another element type. We want to use triangles now."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "initial_condition = initial_condition_weak_blast_wave"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Since there is no direct equivalent to Gauss-Lobatto nodes on triangles, the approximation type\n",
    "`SBP()` now uses Gauss-Lobatto nodes on faces and special nodes in the interior of the triangular.\n",
    "More details can be found in the documentation of\n",
    "[StartUpDG.jl](https://jlchan.github.io/StartUpDG.jl/dev/RefElemData/#RefElemData-based-on-SBP-finite-differences)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dg = DGMulti(polydeg = 3,\n",
    "             element_type = Tri(),\n",
    "             approximation_type = SBP(),\n",
    "             surface_flux = flux_lax_friedrichs,\n",
    "             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))\n",
    "\n",
    "cells_per_dimension = (32, 32)\n",
    "mesh = DGMultiMesh(dg,\n",
    "                   cells_per_dimension, # initial_refinement_level = 5\n",
    "                   coordinates_min=(-2.0, -2.0),\n",
    "                   coordinates_max=( 2.0,  2.0),\n",
    "                   periodicity=true)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,\n",
    "                                    boundary_conditions=boundary_condition_periodic)\n",
    "tspan = (0.0, 0.4)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "alive_callback = AliveCallback(alive_interval=10)\n",
    "analysis_callback = AnalysisCallback(semi, interval=100, uEltype=real(dg))\n",
    "callbacks = CallbackSet(analysis_callback, alive_callback);\n",
    "\n",
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,\n",
    "            ode_default_options()..., callback=callbacks);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "pd = PlotData2D(sol)\n",
    "plot(pd)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plot(pd[\"rho\"])\n",
    "plot!(getmesh(pd))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Triangular meshes on non-Cartesian domains\n",
    "To use triangular meshes on a non-Cartesian domain, Trixi.jl uses the package [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl).\n",
    "The following example is based on [`elixir_euler_triangulate_pkg_mesh.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/dgmulti_2d/elixir_euler_triangulate_pkg_mesh.jl)\n",
    "and uses a pre-defined mesh from StartUpDG.jl."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We want to simulate the smooth initial condition `initial_condition_convergence_test`\n",
    "with source terms `source_terms_convergence_test` for the 2D compressible Euler equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "initial_condition = initial_condition_convergence_test\n",
    "source_terms = source_terms_convergence_test"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We create the solver `DGMulti` with triangular elements (`Tri()`) as before."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dg = DGMulti(polydeg = 3, element_type = Tri(),\n",
    "             approximation_type=Polynomial(),\n",
    "             surface_flux = flux_lax_friedrichs,\n",
    "             volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "StartUpDG.jl provides for instance a pre-defined Triangulate geometry for a rectangular domain with\n",
    "hole `RectangularDomainWithHole`. Other pre-defined Triangulate geometries are e.g., `SquareDomain`,\n",
    "`Scramjet`, and `CircularDomain`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "meshIO = StartUpDG.triangulate_domain(StartUpDG.RectangularDomainWithHole());"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The pre-defined Triangulate geometry in StartUpDG has integer boundary tags. With `DGMultiMesh`\n",
    "we assign boundary faces based on these integer boundary tags and create a mesh compatible with Trixi.jl."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "mesh = DGMultiMesh(dg, meshIO, Dict(:outer_boundary=>1, :inner_boundary=>2))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "boundary_condition_convergence_test = BoundaryConditionDirichlet(initial_condition)\n",
    "boundary_conditions = (; :outer_boundary => boundary_condition_convergence_test,\n",
    "                         :inner_boundary => boundary_condition_convergence_test)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg,\n",
    "                                    source_terms = source_terms,\n",
    "                                    boundary_conditions = boundary_conditions)\n",
    "\n",
    "tspan = (0.0, 0.2)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "alive_callback = AliveCallback(alive_interval=20)\n",
    "analysis_callback = AnalysisCallback(semi, interval=200, uEltype=real(dg))\n",
    "callbacks = CallbackSet(alive_callback, analysis_callback);\n",
    "\n",
    "sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),\n",
    "            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "pd = PlotData2D(sol)\n",
    "plot(pd[\"rho\"])\n",
    "plot!(getmesh(pd))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "For more information, please have a look in the [StartUpDG.jl documentation](https://jlchan.github.io/StartUpDG.jl/stable/)."
   ],
   "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\", \"StartUpDG\", \"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
}