{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 19: Differentiable programming"
   ],
   "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": [
    "[Julia and its ecosystem provide some tools for differentiable programming](https://sinews.siam.org/Details-Page/scientific-machine-learning-how-julia-employs-differentiable-programming-to-do-it-best).\n",
    "Trixi.jl is designed to be flexible, extendable, and composable with Julia's growing ecosystem for\n",
    "scientific computing and machine learning. Thus, the ultimate goal is to have fast implementations\n",
    "that allow automatic differentiation (AD) without too much hassle for users. If some parts do not\n",
    "meet these requirements, please feel free to open an issue or propose a fix in a PR."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In the following, we will walk through some examples demonstrating how to differentiate through\n",
    "Trixi.jl."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Forward mode automatic differentiation"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Trixi.jl integrates well with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl)\n",
    "for forward mode AD."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Computing the Jacobian\n",
    "The high-level interface to compute the Jacobian this way is `jacobian_ad_forward`.\n",
    "First, we load the required packages and compute the Jacobian of a semidiscretization\n",
    "of the compressible Euler equations, a system of nonlinear conservation laws."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, LinearAlgebra, Plots\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "\n",
    "solver = DGSEM(3, flux_central)\n",
    "mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)\n",
    "\n",
    "J = jacobian_ad_forward(semi);\n",
    "size(J)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we compute the eigenvalues of the Jacobian."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "λ = eigvals(J)\n",
    "scatter(real.(λ), imag.(λ), label=\"central flux\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As you can see here, the maximal real part is close to zero."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "relative_maximum = maximum(real, λ) / maximum(abs, λ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Interestingly, if we add dissipation by switching to the `flux_lax_friedrichs`\n",
    "at the interfaces, the maximal real part of the eigenvalues increases."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(3, flux_lax_friedrichs)\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)\n",
    "\n",
    "J = jacobian_ad_forward(semi)\n",
    "λ = eigvals(J)\n",
    "\n",
    "scatter!(real.(λ), imag.(λ), label=\"Lax-Friedrichs flux\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Although the maximal real part is still somewhat small, it's larger than for\n",
    "the purely central discretization."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "relative_maximum = maximum(real, λ) / maximum(abs, λ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "However, we should be careful when using this analysis, since the eigenvectors\n",
    "are not necessarily well-conditioned."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "λ, V = eigen(J)\n",
    "condition_number = cond(V)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In one space dimension, the situation is a bit different."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "equations = CompressibleEulerEquations1D(1.4)\n",
    "\n",
    "solver = DGSEM(3, flux_central)\n",
    "mesh = TreeMesh((-1.0,), (1.0,), initial_refinement_level=6, n_cells_max=10^5)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)\n",
    "\n",
    "J = jacobian_ad_forward(semi)\n",
    "\n",
    "λ = eigvals(J)\n",
    "\n",
    "scatter(real.(λ), imag.(λ), label=\"central flux\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, the maximal real part is basically zero to machine accuracy."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "relative_maximum = maximum(real, λ) / maximum(abs, λ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Moreover, the eigenvectors are not as ill-conditioned as in 2D."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "λ, V = eigen(J)\n",
    "condition_number = cond(V)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "If we add dissipation, the maximal real part is still approximately zero."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(3, flux_lax_friedrichs)\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)\n",
    "\n",
    "J = jacobian_ad_forward(semi)\n",
    "λ = eigvals(J)\n",
    "\n",
    "scatter!(real.(λ), imag.(λ), label=\"Lax-Friedrichs flux\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As you can see from the plot generated above, the maximal real part is still\n",
    "basically zero to machine precision."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "relative_maximum = maximum(real, λ) / maximum(abs, λ)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's check the condition number of the eigenvectors."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "λ, V = eigen(J)\n",
    "\n",
    "condition_number = cond(V)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the condition number of the eigenvector matrix increases but is\n",
    "still smaller than for the example in 2D."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Computing other derivatives"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "It is also possible to compute derivatives of other dependencies using AD in Trixi.jl. For example,\n",
    "you can compute the gradient of an entropy-dissipative semidiscretization with respect to the\n",
    "ideal gas constant of the compressible Euler equations as described in the following. This example\n",
    "is also available as the elixir\n",
    "[`examples/special_elixirs/elixir_euler_ad.jl`](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/special_elixirs/elixir_euler_ad.jl)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "First, we create a semidiscretization of the compressible Euler equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, LinearAlgebra, ForwardDiff\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "\n",
    "\"\"\"\n",
    "    initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)\n",
    "\n",
    "The classical isentropic vortex test case of\n",
    "- Chi-Wang Shu (1997)\n",
    "  Essentially Non-Oscillatory and Weighted Essentially Non-Oscillatory\n",
    "  Schemes for Hyperbolic Conservation Laws\n",
    "  [NASA/CR-97-206253](https://ntrs.nasa.gov/citations/19980007543)\n",
    "\"\"\"\n",
    "function initial_condition_isentropic_vortex(x, t, equations::CompressibleEulerEquations2D)\n",
    "  inicenter = SVector(0.0, 0.0) # initial center of the vortex\n",
    "  iniamplitude = 5.0            # size and strength of the vortex\n",
    "\n",
    "  rho = 1.0  # base flow\n",
    "  v1 = 1.0\n",
    "  v2 = 1.0\n",
    "  vel = SVector(v1, v2)\n",
    "  p = 25.0\n",
    "\n",
    "  rt = p / rho                      # ideal gas equation\n",
    "  t_loc = 0.0\n",
    "\n",
    "  cent = inicenter + vel*t_loc      # shift advection of center to handle periodic BC, but only for v1 = v2 = 1.0\n",
    "  cent = x - cent                   # distance to center point\n",
    "  cent = SVector(-cent[2], cent[1])\n",
    "\n",
    "  r2 = cent[1]^2 + cent[2]^2\n",
    "  du = iniamplitude / (2*π) * exp(0.5 * (1 - r2)) # vel. perturbation\n",
    "  dtemp = -(equations.gamma - 1) / (2 * equations.gamma * rt) * du^2 # isentropic\n",
    "\n",
    "  rho = rho * (1 + dtemp)^(1 / (equations.gamma - 1))\n",
    "  vel = vel + du * cent\n",
    "  v1, v2 = vel\n",
    "  p = p * (1 + dtemp)^(equations.gamma / (equations.gamma - 1))\n",
    "\n",
    "  prim = SVector(rho, v1, v2, p)\n",
    "  return prim2cons(prim, equations)\n",
    "end\n",
    "\n",
    "mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)\n",
    "\n",
    "solver = DGSEM(3, flux_lax_friedrichs, VolumeIntegralFluxDifferencing(flux_ranocha))\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_isentropic_vortex, solver)\n",
    "\n",
    "u0_ode = Trixi.compute_coefficients(0.0, semi)\n",
    "size(u0_ode)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we compute the Jacobian using `ForwardDiff.jacobian`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "J = ForwardDiff.jacobian((du_ode, γ) -> begin\n",
    "    equations_inner = CompressibleEulerEquations2D(first(γ))\n",
    "    semi_inner = Trixi.remake(semi, equations=equations_inner, uEltype=eltype(γ))\n",
    "    Trixi.rhs!(du_ode, u0_ode, semi_inner, 0.0)\n",
    "end, similar(u0_ode), [1.4]); # γ needs to be an `AbstractArray`\n",
    "\n",
    "round.(extrema(J), sigdigits=2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that we create a semidiscretization `semi` at first to determine the state `u0_ode` around\n",
    "which we want to perform the linearization. Next, we wrap the RHS evaluation inside a closure\n",
    "and pass that to `ForwardDiff.jacobian`. There, we need to make sure that the internal caches\n",
    "are able to store dual numbers from ForwardDiff.jl by setting `uEltype` appropriately. A similar\n",
    "approach is used by `jacobian_ad_forward`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note that the ideal gas constant does not influence the semidiscrete rate of change of the\n",
    "density, as demonstrated by"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "norm(J[1:4:end])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Here, we used some knowledge about the internal memory layout of Trixi.jl, an array of structs\n",
    "with the conserved variables as fastest-varying index in memory."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Differentiating through a complete simulation"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "It is also possible to differentiate through a complete simulation. As an example, let's differentiate\n",
    "the total energy of a simulation using the linear scalar advection equation with respect to the\n",
    "wave number (frequency) of the initial data."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq, ForwardDiff, Plots\n",
    "\n",
    "function energy_at_final_time(k) # k is the wave number of the initial condition\n",
    "    equations = LinearScalarAdvectionEquation2D(1.0, -0.3)\n",
    "    mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)\n",
    "    solver = DGSEM(3, flux_lax_friedrichs)\n",
    "    initial_condition = (x, t, equation) -> begin\n",
    "            x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)\n",
    "            return SVector(sinpi(k * sum(x_trans)))\n",
    "    end\n",
    "    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,\n",
    "                                               uEltype=typeof(k))\n",
    "    ode = semidiscretize(semi, (0.0, 1.0))\n",
    "    sol = solve(ode, BS3(), save_everystep=false)\n",
    "    Trixi.integrate(energy_total, sol.u[end], semi)\n",
    "end\n",
    "\n",
    "k_values = range(0.9, 1.1, length=101)\n",
    "\n",
    "plot(k_values, energy_at_final_time.(k_values), label=\"Energy\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "You see a plot of a curve that resembles a parabola with local maximum around `k = 1.0`.\n",
    "Why's that? Well, the domain is fixed but the wave number changes. Thus, if the wave number is\n",
    "not chosen as an integer, the initial condition will not be a smooth periodic function in the\n",
    "given domain. Hence, the dissipative surface flux (`flux_lax_friedrichs` in this example)\n",
    "will introduce more dissipation. In particular, it will introduce more dissipation for \"less smooth\"\n",
    "initial data, corresponding to wave numbers `k` further away from integers."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can compute the discrete derivative of the energy at the final time with respect to the wave\n",
    "number `k` as follows."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "round(ForwardDiff.derivative(energy_at_final_time, 1.0), sigdigits=2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This is rather small and we can treat it as zero in comparison to the value of this derivative at\n",
    "other wave numbers `k`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "dk_values = ForwardDiff.derivative.((energy_at_final_time,), k_values);\n",
    "\n",
    "plot(k_values, dk_values, label=\"Derivative\")"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "If you remember basic calculus, a sufficient condition for a local maximum is that the first derivative\n",
    "vanishes and the second derivative is negative. We can also check this discretely."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "second_derivative = round(ForwardDiff.derivative(\n",
    "        k -> Trixi.ForwardDiff.derivative(energy_at_final_time, k), 1.0),\n",
    "      sigdigits=2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Having seen this application, let's break down what happens step by step."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function energy_at_final_time(k) # k is the wave number of the initial condition\n",
    "    equations = LinearScalarAdvectionEquation2D(1.0, -0.3)\n",
    "    mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)\n",
    "    solver = DGSEM(3, flux_lax_friedrichs)\n",
    "    initial_condition = (x, t, equation) -> begin\n",
    "        x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)\n",
    "        return SVector(sinpi(k * sum(x_trans)))\n",
    "    end\n",
    "    semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,\n",
    "                                               uEltype=typeof(k))\n",
    "    ode = semidiscretize(semi, (0.0, 1.0))\n",
    "    sol = solve(ode, BS3(), save_everystep=false)\n",
    "    Trixi.integrate(energy_total, sol.u[end], semi)\n",
    "end\n",
    "\n",
    "k = 1.0\n",
    "round(ForwardDiff.derivative(energy_at_final_time, k), sigdigits=2)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "When calling `ForwardDiff.derivative(energy_at_final_time, k)` with `k=1.0`, ForwardDiff.jl\n",
    "will basically use the chain rule and known derivatives of existing basic functions\n",
    "to calculate the derivative of the energy at the final time with respect to the\n",
    "wave number `k` at `k0 = 1.0`. To do this, ForwardDiff.jl uses dual numbers, which\n",
    "basically store the result and its derivative w.r.t. a specified parameter at the\n",
    "same time. Thus, we need to make sure that we can treat these `ForwardDiff.Dual`\n",
    "numbers everywhere during the computation. Fortunately, generic Julia code usually\n",
    "supports these operations. The most basic problem for a developer is to ensure\n",
    "that all types are generic enough, in particular the ones of internal caches."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The first step in this example creates some basic ingredients of our simulation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "equations = LinearScalarAdvectionEquation2D(1.0, -0.3)\n",
    "mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=3, n_cells_max=10^4)\n",
    "solver = DGSEM(3, flux_lax_friedrichs);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "These do not have internal caches storing intermediate values of the numerical\n",
    "solution, so we do not need to adapt them. In fact, we could also define them\n",
    "outside of `energy_at_final_time` (but would need to take care of globals or\n",
    "wrap everything in another function)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we define the initial condition"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "initial_condition = (x, t, equation) -> begin\n",
    "    x_trans = Trixi.x_trans_periodic_2d(x - equation.advection_velocity * t)\n",
    "    return SVector(sinpi(k * sum(x_trans)))\n",
    "end;"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "as a closure capturing the wave number `k` passed to `energy_at_final_time`.\n",
    "If you call `energy_at_final_time(1.0)`, `k` will be a `Float64`. Thus, the\n",
    "return values of `initial_condition` will be `SVector`s of `Float64`s. When\n",
    "calculating the `ForwardDiff.derivative`, `k` will be a `ForwardDiff.Dual` number.\n",
    "Hence, the `initial_condition` will return `SVector`s of `ForwardDiff.Dual`\n",
    "numbers."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The semidiscretization `semi` uses some internal caches to avoid repeated allocations\n",
    "and speed up the computations, e.g. for numerical fluxes at interfaces. Thus, we\n",
    "need to tell Trixi.jl to allow `ForwardDiff.Dual` numbers in these caches. That's what\n",
    "the keyword argument `uEltype=typeof(k)` in"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,\n",
    "                                    uEltype=typeof(k));"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "does. This is basically the only part where you need to modify your standard Trixi.jl\n",
    "code to enable automatic differentiation. From there on, the remaining steps"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ode = semidiscretize(semi, (0.0, 1.0))\n",
    "sol = solve(ode, BS3(), save_everystep=false)\n",
    "round(Trixi.integrate(energy_total, sol.u[end], semi), sigdigits=5)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "do not need any modifications since they are sufficiently generic (and enough effort\n",
    "has been spend to allow general types inside these calls)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Propagating errors using Measurements.jl"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "[![Error bars by Randall Munroe](https://imgs.xkcd.com/comics/error_bars.png)](https://xkcd.com/2110/)\n",
    "\"Error bars\" by Randall Munroe, linked from https://xkcd.com/2110"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Similar to AD, Trixi.jl also allows propagating uncertainties using linear error propagation\n",
    "theory via [Measurements.jl](https://github.com/JuliaPhysics/Measurements.jl).\n",
    "As an example, let's create a system representing the linear advection equation\n",
    "in 1D with an uncertain velocity. Then, we create a semidiscretization using a\n",
    "sine wave as initial condition, solve the ODE, and plot the resulting uncertainties\n",
    "in the primitive variables."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq, Measurements, Plots, LaTeXStrings\n",
    "\n",
    "equations = LinearScalarAdvectionEquation1D(1.0 ± 0.1)\n",
    "\n",
    "mesh = TreeMesh((-1.0,), (1.0,), n_cells_max=10^5, initial_refinement_level=5)\n",
    "\n",
    "solver = DGSEM(3)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,\n",
    "                                    solver, uEltype=Measurement{Float64})\n",
    "\n",
    "ode = semidiscretize(semi, (0.0, 1.5))\n",
    "\n",
    "sol = solve(ode, BS3(), save_everystep=false);\n",
    "\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "You should see a plot where small error bars are shown around the extrema\n",
    "and larger error bars are shown in the remaining parts.\n",
    "This result is in accordance with expectations. Indeed, the uncertain propagation\n",
    "speed will affect the extrema less since the local variation of the solution is\n",
    "relatively small there. In contrast, the local variation of the solution is large\n",
    "around the turning points of the sine wave, so the uncertainties will be relatively\n",
    "large there."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "All this is possible due to allowing generic types and having good abstractions\n",
    "in Julia that allow packages to work together seamlessly."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Finite difference approximations"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Trixi.jl provides the convenience function `jacobian_fd` to approximate the Jacobian\n",
    "via central finite differences."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, LinearAlgebra\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "\n",
    "solver = DGSEM(3, flux_central)\n",
    "\n",
    "mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_density_wave, solver)\n",
    "\n",
    "J_fd = jacobian_fd(semi)\n",
    "\n",
    "J_ad = jacobian_ad_forward(semi)\n",
    "\n",
    "relative_difference = norm(J_fd - J_ad) / size(J_fd, 1)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "This discrepancy is of the expected order of magnitude for central finite difference approximations."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Linear systems"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "When a linear PDE is discretized using a linear scheme such as a standard DG method,\n",
    "the resulting semidiscretization yields an affine ODE of the form"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "$$\n",
    "\\partial_t u(t) = A u(t) + b,\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "where `A` is a linear operator (\"matrix\") and `b` is a vector. Trixi.jl allows you\n",
    "to obtain this linear structure in a matrix-free way by using `linear_structure`.\n",
    "The resulting operator `A` can be used in multiplication, e.g. `mul!` from\n",
    "LinearAlgebra, converted to a sparse matrix using `sparse` from SparseArrays,\n",
    "or converted to a dense matrix using `Matrix` for detailed eigenvalue analyses.\n",
    "For example,"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, LinearAlgebra, Plots\n",
    "\n",
    "equations = LinearScalarAdvectionEquation2D(1.0, -0.3)\n",
    "\n",
    "solver = DGSEM(3, flux_lax_friedrichs)\n",
    "\n",
    "mesh = TreeMesh((-1.0, -1.0), (1.0, 1.0), initial_refinement_level=2, n_cells_max=10^5)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test, solver)\n",
    "\n",
    "A, b = linear_structure(semi)\n",
    "\n",
    "size(A), size(b)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we compute the eigenvalues of the linear operator."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "λ = eigvals(Matrix(A))\n",
    "\n",
    "scatter(real.(λ), imag.(λ))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As you can see here, the maximal real part is close to machine precision."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "λ = eigvals(Matrix(A))\n",
    "relative_maximum = maximum(real, λ) / maximum(abs, λ)"
   ],
   "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\", \"ForwardDiff\"],\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
}