{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 3: Introduction to DG methods"
   ],
   "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": [
    "This tutorial is about how to set up a simple way to approximate the solution of a hyperbolic partial\n",
    "differential equation. First, we will implement a basic and naive algorithm. Then, we will use predefined\n",
    "features from [Trixi.jl](https://github.com/trixi-framework/Trixi.jl) to show how you can use Trixi.jl on your own."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will implement the scalar linear advection equation in 1D with the advection velocity $1$.\n",
    "$$\n",
    "u_t + u_x = 0,\\; \\text{for} \\;t\\in \\mathbb{R}^+, x\\in\\Omega=[-1,1]\n",
    "$$\n",
    "We define the domain $\\Omega$ by setting the boundaries."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "coordinates_min = -1.0 # minimum coordinate\n",
    "coordinates_max = 1.0  # maximum coordinate"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We assume periodic boundaries and the following initial condition."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "initial_condition_sine_wave(x) = 1.0 + 0.5 * sin(pi * x)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## The discontinuous Galerkin collocation spectral element method (DGSEM)\n",
    "### i. Discretization of the physical domain\n",
    "To improve precision we want to approximate the solution on small parts of the physical domain.\n",
    "So, we split the domain $\\Omega=[-1, 1]$ into elements $Q_l$ of length $dx$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "n_elements = 16 # number of elements\n",
    "\n",
    "dx = (coordinates_max - coordinates_min) / n_elements # length of one element"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To make the calculation more efficient and storing less information, we transform each element\n",
    "$Q_l$ with center point $x_l$ to a reference element $E=[-1, 1]$\n",
    "$$\n",
    "Q_l=\\Big[x_l-\\frac{dx}{2}, x_l+\\frac{dx}{2}\\Big] \\underset{x(\\xi)}{\\overset{\\xi(x)}{\\rightleftarrows}} [-1, 1].\n",
    "$$\n",
    "So, for every element the transformation from the reference domain to the physical domain is defined by\n",
    "$$\n",
    "x(\\xi) = x_l + \\frac{dx}{2} \\xi,\\; \\xi\\in[-1, 1]\n",
    "$$\n",
    "Therefore,\n",
    "$$\n",
    "\\begin{align*}\n",
    "u &= u(x(\\xi), t) \\\\\n",
    "u_x &= u_\\xi \\frac{d\\xi}{dx} \\\\[3pt]\n",
    "\\frac{d\\xi}{dx} &= (x_\\xi)^{-1} = \\frac{2}{dx} =: J^{-1}. \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "Here, $J$ is the Jacobian determinant of the transformation."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Using this transformation, we can transform our equation for each element $Q_l$.\n",
    "$$\n",
    "\\frac{dx}{2} u_t^{Q_l} + u_\\xi^{Q_l} = 0 \\text{, for }t\\in\\mathbb{R}^+,\\; \\xi\\in[-1, 1]\n",
    "$$\n",
    "Here, $u_t^{Q_l}$ and $u_\\xi^{Q_l}$ denote the time and spatial derivatives of the solution on the element $Q_l$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### ii. Polynomial approach\n",
    "Now, we want to approximate the solution in each element $Q_l$ by a polynomial of degree $N$. Since we transformed\n",
    "the equation, we can use the same polynomial approach for the reference coordinate $\\xi\\in[-1, 1]$ in every\n",
    "physical element $Q_l$. This saves a lot of resources by reducing the amount of calculations needed\n",
    "and storing less information."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For DGSEM we choose [Lagrange basis functions](https://en.wikipedia.org/wiki/Lagrange_polynomial)\n",
    "$\\{l_j\\}_{j=0}^N$ as our polynomial basis of degree $N$ in $[-1, 1]$.\n",
    "The solution in element $Q_l$ can be approximated by\n",
    "$$\n",
    "u(x(\\xi), t)\\big|_{Q_l} \\approx u^{Q_l}(\\xi, t) = \\sum_{j=0}^N u_j^{Q_l}(t) l_j(\\xi)\n",
    "$$\n",
    "with $N+1$ coefficients $\\{u_j^{Q_l}\\}_{j=0}^N$.\n",
    "By construction the Lagrange basis has some useful advantages. This basis is defined by $N+1$ nodes, which\n",
    "fulfill a Kronecker property at the exact same nodes. Let $\\{\\xi_i\\}_{i=0}^N$ be these nodes.\n",
    "$$\n",
    "l_j(\\xi_i) = \\delta_{i,j} =\n",
    "\\begin{cases}\n",
    "1, & \\text{if } i=j \\\\\n",
    "0, & \\text{else.}\n",
    "\\end{cases}\n",
    "$$\n",
    "Because of this property, the polynomial coefficients are exact the values of $u^{Q_l}$ at the nodes\n",
    "$$\n",
    "u^{Q_l}(\\xi_i, t) = \\sum_{j=0}^N u_j^{Q_l}(t) \\underbrace{l_j(\\xi_i)}_{=\\delta_{ij}} = u_i^{Q_l}(t).\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we want to select the nodes $\\{\\xi_i\\}_{i=0}^N$, which we use for the construction of the Lagrange\n",
    "polynomials. We choose the $N+1$ Gauss-Lobatto nodes, which are used for the\n",
    "[Gaussian-Lobatto quadrature](https://mathworld.wolfram.com/LobattoQuadrature.html).\n",
    "These always contain the boundary points at $-1$ and $+1$ and are well suited as interpolation nodes.\n",
    "The corresponding weights will be referred to as $\\{w_j\\}_{j=0}^N$.\n",
    "In Trixi.jl the basis with Lagrange polynomials on Gauss-Lobatto nodes is already defined."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi\n",
    "polydeg = 3 #= polynomial degree = N =#\n",
    "basis = LobattoLegendreBasis(polydeg)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Gauss-Lobatto nodes are"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "nodes = basis.nodes"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "with the corresponding weights"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "weights = basis.weights"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To illustrate how you can integrate using numerical quadrature with this Legendre-Gauss-Lobatto nodes,\n",
    "we give an example for $f(x)=x^3$. Since $f$ is of degree $3$, a polynomial interpolation with $N=3$ is exact.\n",
    "Therefore, the integral on $[-1, 1]$ can be calculated by\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\int_{-1}^1 f(x) dx &= \\int_{-1}^1 \\Big( \\sum_{j=0}^3 f(\\xi_j)l_j(x) \\Big) dx\n",
    "= \\sum_{j=0}^3 f(\\xi_j) \\int_{-1}^1 l_j(x)dx \\\\\n",
    "&=: \\sum_{j=0}^3 f(\\xi_j) w_j\n",
    "= \\sum_{j=0}^3 \\xi_j^3 w_j\n",
    "\\end{align*}\n",
    "$$\n",
    "Let's use our nodes and weights for $N=3$ and plug in"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "integral = sum(nodes.^3 .* weights)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Using this polynomial approach leads to the equation\n",
    "$$\n",
    "\\frac{dx}{2} \\dot{u}^{Q_l}(\\xi, t) + u^{Q_l}(\\xi, t)' = 0\n",
    "$$\n",
    "with $\\dot{u}=\\frac{\\partial}{\\partial t}u$ and $u'=\\frac{\\partial}{\\partial x}u$.\n",
    "To approximate the solution, we need to get the polynomial coefficients $\\{u_j^{Q_l}\\}_{j=0}^N$\n",
    "for every element $Q_l$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "After defining all nodes, we can implement the spatial coordinate $x$ and its initial value $u0 = u(t_0)$\n",
    "for every node."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "x = Matrix{Float64}(undef, length(nodes), n_elements)\n",
    "for element in 1:n_elements\n",
    "    x_l = coordinates_min + (element - 1) * dx + dx/2\n",
    "    for i in eachindex(nodes)\n",
    "        ξ = nodes[i] # nodes in [-1, 1]\n",
    "        x[i, element] = x_l + dx/2 * ξ\n",
    "    end\n",
    "end\n",
    "\n",
    "u0 = initial_condition_sine_wave.(x)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To have a look at the initial sinus curve, we plot it."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(vec(x), vec(u0), label=\"initial condition\", legend=:topleft)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "### iii. Variational formulation\n",
    "After defining the equation and initial condition, we want to implement an algorithm to\n",
    "approximate the solution."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "From now on, we only write $u$ instead of $u^{Q_l}$ for simplicity, but consider that all the following\n",
    "calculation only concern one element.\n",
    "Multiplying the new equation with the smooth Lagrange polynomials $\\{l_i\\}_{i=0}^N$ (test functions)\n",
    "and integrating over the reference element $E=[-1,1]$, we get the variational formulation of our\n",
    "transformed partial differential equation for $i=0,...,N$:\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\int_{-1}^1 \\Big( \\frac{dx}{2} \\dot{u}(\\xi, t) + u'(\\xi, t) \\Big) l_i(\\xi)d\\xi\n",
    "  &= \\underbrace{\\frac{dx}{2} \\int_{-1}^1 \\dot{u}(\\xi, t) l_i(\\xi)d\\xi}_{\\text{Term I}} + \\underbrace{\\int_{-1}^1 u'(\\xi, t) l_i(\\xi)d\\xi}_{\\text{Term II}} = 0\n",
    "\\end{align*}\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We deal with the two terms separately. We write $\\int_{-1, N}^1 \\;\\cdot\\; d\\xi$ for the approximation\n",
    "of the integral using numerical quadrature with $N+1$ basis points. We use the Gauss-Lobatto nodes\n",
    "again. The numerical scalar product $\\langle\\cdot, \\cdot\\rangle_N$ is defined by\n",
    "$\\langle f, g\\rangle_N := \\int_{-1, N}^1 f(\\xi) g(\\xi) d\\xi$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Term I:\n",
    "In the following calculation we approximate the integral numerically with quadrature on the Gauss-Lobatto\n",
    "nodes $\\{\\xi_i\\}_{i=0}^N$ and then use the Kronecker property of the Lagrange polynomials. This approach\n",
    "of using the same nodes for the interpolation and quadrature is called collocation.\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\frac{dx}{2} \\int_{-1}^1 \\dot{u}(\\xi, t) l_i(\\xi)d\\xi\n",
    "&\\approx \\frac{dx}{2} \\int_{-1, N}^1 \\dot{u}(\\xi, t) l_i(\\xi)d\\xi \\\\\n",
    "&= \\frac{dx}{2} \\sum_{k=0}^N \\underbrace{\\dot{u}(\\xi_k, t)}_{=\\dot{u}_k(t)} \\underbrace{l_i(\\xi_k)}_{=\\delta_{k,i}}w_k \\\\\n",
    "&= \\frac{dx}{2} \\dot{u}_i(t) w_i\n",
    "\\end{align*}\n",
    "$$\n",
    "We define the Legendre-Gauss-Lobatto (LGL) mass matrix $M$ and by the Kronecker property follows:\n",
    "$$\n",
    "M_{ij} = \\langle l_j, l_i\\rangle_N = \\delta_{ij} w_j,\\; i,j=0,...,N.\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using LinearAlgebra\n",
    "M = diagm(weights)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we can write the integral with this new matrix.\n",
    "$$\n",
    "\\frac{dx}{2} \\int_{-1, N}^1 \\dot{u}(\\xi, t) \\underline{l}(\\xi)d\\xi = \\frac{dx}{2} M \\underline{\\dot{u}}(t),\n",
    "$$\n",
    "where $\\underline{\\dot{u}} = (\\dot{u}_0, ..., \\dot{u}_N)^T$ and $\\underline{l}$ respectively."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note:** Since the LGL quadrature with $N+1$ nodes is exact up to functions of degree $2N-1$ and\n",
    "$\\dot{u}(\\xi, t) l_i(\\xi)$ is of degree $2N$, in general the following holds\n",
    "$$\n",
    "\\int_{-1}^1 \\dot{u}(\\xi, t) l_i(\\xi) d\\xi \\neq \\int_{-1, N}^1 \\dot{u}(\\xi, t) l_i(\\xi) d\\xi.\n",
    "$$\n",
    "With an exact integration the mass matrix would be dense. Choosing numerical integrating and quadrature\n",
    "with the exact same nodes (collocation) leads to the sparse and diagonal mass matrix $M$. This\n",
    "is called mass lumping and has the big advantage of an easy invertation of the matrix."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Term II:\n",
    "We use spatial partial integration for the second term:\n",
    "$$\n",
    "\\int_{-1}^1 u'(\\xi, t) l_i(\\xi) d\\xi = [u l_i]_{-1}^1 - \\int_{-1}^1 u l_i'd\\xi\n",
    "$$\n",
    "The resulting integral can be solved exactly with LGL quadrature since the polynomial is now\n",
    "of degree $2N-1$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Again, we split the calculation in two steps."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Surface term\n",
    "As mentioned before, we approximate the solution with a polynomial in every element. Therefore, in\n",
    "general the value of this approximation at the interfaces between two elements is not unique. To solve\n",
    "this problem we introduce the idea of the numerical flux $u^*$, which will give an exact value at\n",
    "the interfaces. One of many different approaches and definitions for the calculation of the\n",
    "numerical flux we will deal with in 4. Numerical flux.\n",
    "$$\n",
    "[u l_i]_{-1}^1 = u^*\\big|^1 l_i(+1) - u^*\\big|_{-1} l_i(-1)\n",
    "$$\n",
    "Since the Gauss-Lobatto nodes contain the element boundaries $-1$ and $+1$, we can use the\n",
    "Kronecker property of $l_i$ for the calculation of $l_i(-1)$ and $l_i(+1)$.\n",
    "$$\n",
    "[u \\underline{l}]_{-1}^1 = u^*\\big|^1 \\left(\\begin{array}{c} 0 \\\\ \\vdots \\\\ 0 \\\\ 1 \\end{array}\\right)\n",
    "- u^*\\big|_{-1} \\left(\\begin{array}{c} 1 \\\\ 0 \\\\ \\vdots \\\\ 0\\end{array}\\right)\n",
    "= B \\underline{u}^*(t)\n",
    "$$\n",
    "with the boundary matrix\n",
    "$$\n",
    "B = \\begin{pmatrix}\n",
    "-1 & 0 & \\cdots & 0\\\\\n",
    "0 & 0 & \\cdots & 0\\\\\n",
    "\\vdots & \\vdots & 0 & 0\\\\\n",
    "0 & \\cdots & 0 & 1\n",
    "\\end{pmatrix}\n",
    "\\qquad\\text{and}\\qquad\n",
    "\\underline{u}^*(t) = \\left(\\begin{array}{c} u^*\\big|_{-1} \\\\ 0 \\\\ \\vdots \\\\ 0 \\\\ u^*\\big|^1\\end{array}\\right).\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "B = diagm([-1; zeros(polydeg - 1); 1])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Volume term\n",
    "As mentioned before, the new integral can be solved exact since the function inside is of degree $2N-1$.\n",
    "$$\n",
    "- \\int_{-1}^1 u l_i'd\\xi = - \\int_{-1, N}^1 u l_i' d\\xi\n",
    "= - \\sum_{k=0}^N u(\\xi_k, t) l_i'(\\xi_k) w_k\n",
    "= - \\sum_{k=0}^N u_k(t) D_{ki} w_k\n",
    "$$\n",
    "where $D$ is the derivative matrix defined by $D_{ki} = l_i'(\\xi_k)$ for $i,k=0,...,N$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "D = basis.derivative_matrix"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To show why this matrix is called the derivative matrix, we go back to our example $f(x)=x^3$.\n",
    "We calculate the derivation of $f$ at the Gauss-Lobatto nodes $\\{\\xi_k\\}_{k=0}^N$ with $N=8$.\n",
    "$$\n",
    "f'|_{x=\\xi_k} = \\Big( \\sum_{j=0}^8 f(\\xi_j) l_j(x) \\Big)'|_{x=\\xi_k} = \\sum_{j=0}^8 f(\\xi_j) l_j'(\\xi_k)\n",
    "= \\sum_{j=0}^8 f(\\xi_j) D_{kj}\n",
    "$$\n",
    "for $k=0,...,N$ and therefore, $\\underline{f}' = D \\underline{f}$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "basis_N8 = LobattoLegendreBasis(8)\n",
    "plot(vec(x), x -> 3 * x^2, label=\"f'\", lw=2)\n",
    "scatter!(basis_N8.nodes, basis_N8.derivative_matrix * basis_N8.nodes.^3, label=\"Df\", lw=3)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Combining the volume term for every $i=0,...,N$ results in\n",
    "$$\n",
    "\\int_{-1}^1 u \\underline{l'} d\\xi = - D^T M \\underline{u}(t)\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Putting all parts together we get the following equation for the element $Q_l$\n",
    "$$\n",
    "\\frac{dx}{2} M \\underline{\\dot{u}}(t) = - B \\underline{u}^*(t) + D^T M \\underline{u}(t)\n",
    "$$\n",
    "or equivalent\n",
    "$$\n",
    "\\underline{\\dot{u}}^{Q_l}(t) = \\frac{2}{dx} \\Big[ - M^{-1} B \\underline{u}^{{Q_l}^*}(t) + M^{-1} D^T M \\underline{u}^{Q_l}(t)\\Big].\n",
    "$$\n",
    "This is called the weak form of the DGSEM."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note:** For every element $Q_l$ we get a system of $N+1$ ordinary differential equations to\n",
    "calculate $N+1$ coefficients. Since the numerical flux $u^*$ is depending on extern values at\n",
    "the interfaces, the equation systems of adjacent elements are weakly linked."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### iv. Numerical flux\n",
    "As mentioned above, we still have to handle the problem of different values at the same point at\n",
    "the interfaces. This happens with the ideas of the numerical flux $f^*(u)=u^*$. The role of $f^*$\n",
    "might seem minor in this simple example, but is important for more complicated problems.\n",
    "There are two values at the same spatial coordinate. Let's say we are looking at the interface between\n",
    "the elements $Q_l$ and $Q_{l+1}$, while both elements got $N+1$ nodes as defined before. We call\n",
    "the first value of the right element $u_R=u_0^{Q_{l+1}}$ and the last one of the left element\n",
    "$u_L=u_N^{Q_l}$. So, for the value of the numerical flux on that interface the following holds\n",
    "$$\n",
    "u^* = u^*(u_L, u_R).\n",
    "$$\n",
    "These values are interpreted as start values of a so-called Riemann problem. There are many\n",
    "different (approximate) Riemann solvers available and useful for different problems. We will\n",
    "use the local Lax-Friedrichs flux."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "surface_flux = flux_lax_friedrichs"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The only missing ingredient is the flux calculation at the boundaries $-1$ and $+1$.\n",
    "$$\n",
    "u^{{Q_{first}}^*}\\big|_{-1} = u^{{Q_{first}}^*}\\big|_{-1}(u^{bound}(-1), u_R)\n",
    "\\quad\\text{and}\\quad\n",
    "u^{{Q_{last}}^*}\\big|^1 = u^{{Q_{last}}^*}\\big|^1(u_L, u^{bound}(1))\n",
    "$$\n",
    "The boundaries are periodic, which means that the last value of the last element $u^{Q_{last}}_N$\n",
    "is used as $u_L$ at the first interface and accordingly for the other boundary."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we implement a function, that calculates $\\underline{\\dot{u}}^{Q_l}$ for the given matrices,\n",
    "$\\underline{u}$ and $\\underline{u}^*$."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function rhs!(du, u, x, t)\n",
    "    # Reset du and flux matrix\n",
    "    du .= zero(eltype(du))\n",
    "    flux_numerical = copy(du)\n",
    "\n",
    "    # Calculate interface and boundary fluxes, $u^* = (u^*|_{-1}, 0, ..., 0, u^*|^1)^T$\n",
    "    # Since we use the flux Lax-Friedrichs from Trixi.jl, we have to pass some extra arguments.\n",
    "    # Trixi.jl needs the equation we are dealing with and an additional `1`, that indicates the\n",
    "    # first coordinate direction.\n",
    "    equations = LinearScalarAdvectionEquation1D(1.0)\n",
    "    for element in 2:n_elements-1\n",
    "        # left interface\n",
    "        flux_numerical[1, element] = surface_flux(u[end, element-1], u[1, element], 1, equations)\n",
    "        flux_numerical[end, element-1] = flux_numerical[1, element]\n",
    "        # right interface\n",
    "        flux_numerical[end, element] = surface_flux(u[end, element], u[1, element+1], 1, equations)\n",
    "        flux_numerical[1, element+1] = flux_numerical[end, element]\n",
    "    end\n",
    "    # boundary flux\n",
    "    flux_numerical[1, 1] = surface_flux(u[end, end], u[1, 1], 1, equations)\n",
    "    flux_numerical[end, end] = flux_numerical[1, 1]\n",
    "\n",
    "    # Calculate surface integrals, $- M^{-1} * B * u^*$\n",
    "    for element in 1:n_elements\n",
    "        du[:, element] -= (M \\ B) * flux_numerical[:, element]\n",
    "    end\n",
    "\n",
    "    # Calculate volume integral, $+ M^{-1} * D^T * M * u$\n",
    "    for element in 1:n_elements\n",
    "        flux = u[:, element]\n",
    "        du[:, element] += (M \\ transpose(D)) * M * flux\n",
    "    end\n",
    "\n",
    "    # Apply Jacobian from mapping to reference element\n",
    "    for element in 1:n_elements\n",
    "        du[:, element] *= 2 / dx\n",
    "    end\n",
    "\n",
    "    return nothing\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Combining all definitions and the function that calculates the right-hand side, we define the ODE and\n",
    "solve it until `t=2` with OrdinaryDiffEq's `solve` function and the Runge-Kutta method `RDPK3SpFSAL49()`,\n",
    "which is optimized for discontinuous Galerkin methods and hyperbolic PDEs. We set some common\n",
    "error tolerances `abstol=1.0e-6, reltol=1.0e-6` and pass `save_everystep=false` to avoid saving intermediate\n",
    "solution vectors in memory."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "tspan = (0.0, 2.0)\n",
    "ode = ODEProblem(rhs!, u0, tspan, x)\n",
    "\n",
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6, ode_default_options()...)\n",
    "\n",
    "plot(vec(x), vec(sol.u[end]), label=\"solution at t=$(tspan[2])\", legend=:topleft, lw=3)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Alternative Implementation based on Trixi.jl\n",
    "Now, we implement the same example. But this time, we directly use the functionality that Trixi.jl\n",
    "provides."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq, Plots"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "First, define the equation with a advection_velocity of `1`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "advection_velocity = 1.0\n",
    "equations = LinearScalarAdvectionEquation1D(advection_velocity)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, create a DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux.\n",
    "The implementation of the basis and the numerical flux is now already done."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will now create a mesh with 16 elements for the physical domain `[-1, 1]` with periodic boundaries.\n",
    "We use Trixi.jl's standard mesh `TreeMesh`. Since it's limited to hypercube domains, we\n",
    "choose `2^4=16` elements. The mesh type supports AMR, that' why `n_cells_max` has to be set, even\n",
    "if we don't need AMR here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "coordinates_min = -1.0 # minimum coordinate\n",
    "coordinates_max = 1.0  # maximum coordinate\n",
    "mesh = TreeMesh(coordinates_min, coordinates_max,\n",
    "                initial_refinement_level=4, # number of elements = 2^4\n",
    "                n_cells_max=30_000) # set maximum capacity of tree data structure (only needed for AMR)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "A semidiscretization collects data structures and functions for the spatial discretization.\n",
    "In Trixi.jl, an initial condition has the following parameter structure and is of the type `SVector`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "initial_condition_sine_wave(x, t, equations) = SVector(1.0 + 0.5 * sin(pi * sum(x - equations.advection_velocity * t)))\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sine_wave, solver)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Again, combining all definitions and the function that calculates the right-hand side, we define the ODE and\n",
    "solve it until `t=2` with OrdinaryDiffEq's `solve` function and the Runge-Kutta method `RDPK3SpFSAL49()`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tspan = (0.0, 2.0)\n",
    "ode_trixi  = semidiscretize(semi, tspan)\n",
    "\n",
    "sol_trixi = solve(ode_trixi, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6, ode_default_options()...);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We add a plot of the new approximated solution to the one calculated before."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "plot!(sol_trixi, label=\"solution at t=$(tspan[2]) with Trixi.jl\", legend=:topleft, linestyle=:dash, lw=2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Summary of the code\n",
    "To sum up, here is the complete code that we used."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Raw implementation"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "# basis: Legendre-Gauss-Lobatto\n",
    "using Trixi, LinearAlgebra, OrdinaryDiffEq, Plots\n",
    "polydeg = 3 #= polynomial degree =#\n",
    "basis = LobattoLegendreBasis(polydeg)\n",
    "nodes = basis.nodes # Gauss-Lobatto nodes in [-1, 1]\n",
    "D = basis.derivative_matrix\n",
    "M = diagm(basis.weights) # mass matrix\n",
    "B = diagm([-1; zeros(polydeg - 1); 1])\n",
    "\n",
    "# mesh\n",
    "coordinates_min = -1.0 # minimum coordinate\n",
    "coordinates_max = 1.0  # maximum coordinate\n",
    "n_elements      = 16   # number of elements\n",
    "\n",
    "dx = (coordinates_max - coordinates_min) / n_elements # length of one element\n",
    "\n",
    "x = Matrix{Float64}(undef, length(nodes), n_elements)\n",
    "for element in 1:n_elements\n",
    "    x_l = -1 + (element - 1) * dx + dx/2\n",
    "    for i in eachindex(nodes) # basis points in [-1, 1]\n",
    "        ξ = nodes[i]\n",
    "        x[i, element] = x_l + dx/2 * ξ\n",
    "    end\n",
    "end\n",
    "\n",
    "# initial condition\n",
    "initial_condition_sine_wave(x) = 1.0 + 0.5 * sin(pi * x)\n",
    "u0 = initial_condition_sine_wave.(x)\n",
    "\n",
    "plot(vec(x), vec(u0), label=\"initial condition\", legend=:topleft)\n",
    "\n",
    "# flux Lax-Friedrichs\n",
    "surface_flux = flux_lax_friedrichs\n",
    "\n",
    "# rhs! method\n",
    "function rhs!(du, u, x, t)\n",
    "    # reset du\n",
    "    du .= zero(eltype(du))\n",
    "    flux_numerical = copy(du)\n",
    "\n",
    "    # calculate interface and boundary fluxes\n",
    "    equations = LinearScalarAdvectionEquation1D(1.0)\n",
    "    for element in 2:n_elements-1\n",
    "        # left interface\n",
    "        flux_numerical[1, element] = surface_flux(u[end, element-1], u[1, element], 1, equations)\n",
    "        flux_numerical[end, element-1] = flux_numerical[1, element]\n",
    "        # right interface\n",
    "        flux_numerical[end, element] = surface_flux(u[end, element], u[1, element+1], 1, equations)\n",
    "        flux_numerical[1, element+1] = flux_numerical[end, element]\n",
    "    end\n",
    "    # boundary flux\n",
    "    flux_numerical[1, 1] = surface_flux(u[end, end], u[1, 1], 1, equations)\n",
    "    flux_numerical[end, end] = flux_numerical[1, 1]\n",
    "\n",
    "    # calculate surface integrals\n",
    "    for element in 1:n_elements\n",
    "        du[:, element] -= (M \\ B) * flux_numerical[:, element]\n",
    "    end\n",
    "\n",
    "    # calculate volume integral\n",
    "    for element in 1:n_elements\n",
    "        flux = u[:, element]\n",
    "        du[:, element] += (M \\ transpose(D)) * M * flux\n",
    "    end\n",
    "\n",
    "    # apply Jacobian from mapping to reference element\n",
    "    for element in 1:n_elements\n",
    "        du[:, element] *= 2 / dx\n",
    "    end\n",
    "\n",
    "    return nothing\n",
    "end\n",
    "\n",
    "# create ODE problem\n",
    "tspan = (0.0, 2.0)\n",
    "ode = ODEProblem(rhs!, u0, tspan, x)\n",
    "\n",
    "# solve\n",
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6, ode_default_options()...)\n",
    "\n",
    "plot(vec(x), vec(sol.u[end]), label=\"solution at t=$(tspan[2])\", legend=:topleft, lw=3)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Alternative Implementation based on Trixi.jl"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq, Plots\n",
    "\n",
    "# equation with a advection_velocity of `1`.\n",
    "advection_velocity = 1.0\n",
    "equations = LinearScalarAdvectionEquation1D(advection_velocity)\n",
    "\n",
    "# create DG solver with flux lax friedrichs and LGL basis\n",
    "solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)\n",
    "\n",
    "# distretize domain with `TreeMesh`\n",
    "coordinates_min = -1.0 # minimum coordinate\n",
    "coordinates_max = 1.0 # maximum coordinate\n",
    "mesh = TreeMesh(coordinates_min, coordinates_max,\n",
    "                initial_refinement_level=4, # number of elements = 2^4\n",
    "                n_cells_max=30_000)\n",
    "\n",
    "# create initial condition and semidiscretization\n",
    "initial_condition_sine_wave(x, t, equations) = SVector(1.0 + 0.5 * sin(pi * sum(x - equations.advection_velocity * t)))\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_sine_wave, solver)\n",
    "\n",
    "# solve\n",
    "tspan = (0.0, 2.0)\n",
    "ode_trixi  = semidiscretize(semi, tspan)\n",
    "sol_trixi  = solve(ode_trixi, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6, ode_default_options()...);\n",
    "\n",
    "plot!(sol_trixi, label=\"solution at t=$(tspan[2]) with Trixi.jl\", legend=:topleft, linestyle=:dash, lw=2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "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.9.4"
  },
  "kernelspec": {
   "name": "julia-1.9",
   "display_name": "Julia 1.9.4",
   "language": "julia"
  }
 },
 "nbformat": 4
}