{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 10: Adding a new scalar conservation law"
   ],
   "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": [
    "If you want to use Trixi.jl for your own research, you might be interested in\n",
    "a new physics model that's not already included in Trixi.jl. In this tutorial,\n",
    "we will implement the cubic conservation law\n",
    "$$\n",
    "\\partial_t u(t,x) + \\partial_x u(t,x)^3 = 0\n",
    "$$\n",
    "in a periodic domain in one space dimension. In Trixi.jl, such a mathematical model\n",
    "is encoded as a subtype of `Trixi.AbstractEquations`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Basic setup"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi\n",
    "\n",
    "struct CubicEquation <: Trixi.AbstractEquations{1 #= number of spatial dimensions =#,\n",
    "                                                1 #= number of primary variables, i.e. scalar =#};\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We create `CubicEquation` as an empty `struct` since we do not use any parameters\n",
    "for this equation. Other models could bundle arbitrary parameters, e.g., the\n",
    "ideal gas constant for the compressible Euler equations."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we define the physical flux `f(u) = u^3` using the calling structure\n",
    "used in Trixi.jl."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "Trixi.flux(u, orientation, equation::CubicEquation) = u.^3\n",
    "Trixi.varnames(_, ::CubicEquation) = (\"scalar\",)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In Trixi.jl, the conserved variables `u` are usually passed as `SVector`s of variables\n",
    "at a single physical location. Hence, we must use `u.^3` instead of the scalar\n",
    "operation `u^3`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "That's already enough to run a simple simulation with a standard DGSEM discretization\n",
    "using the non-dissipative central flux at interfaces."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "\n",
    "# Create a simulation setup\n",
    "equation = CubicEquation()\n",
    "\n",
    "initial_condition_sine(x, t, equation::CubicEquation) = SVector(sinpi(x[1]))\n",
    "\n",
    "mesh = TreeMesh(-1.0, 1.0, # min/max coordinates\n",
    "                initial_refinement_level=4,\n",
    "                n_cells_max=10^4)\n",
    "\n",
    "solver = DGSEM(3 #= polynomial degree =#, flux_central)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equation, initial_condition_sine, solver)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We wrap the return value of the `initial_condition_sine` inside an `SVector` since that's the approach\n",
    "used in Trixi.jl also for systems of equations. We need to index the spatial coordinate `x[1]`,\n",
    "since it is an `SVector` with one component. In multiple space dimensions, all spatial coordinates\n",
    "are passed together."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Next, we create an `ODEProblem` from the SciML/DifferentialEquations ecosystem.\n",
    "We can solve this ODE numerically using any time integration method,\n",
    "e.g., `SSPRK43` from [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl).\n",
    "Before, we set up a callback to summarize the simulation setup."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "# Create ODE problem with given time span\n",
    "tspan = (0.0, 0.09)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "summary_callback = SummaryCallback()\n",
    "callbacks = CallbackSet(summary_callback)\n",
    "\n",
    "# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks\n",
    "sol = solve(ode, SSPRK43();\n",
    "            ode_default_options()..., callback=callbacks);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "That's it, you ran your first simulation using your new equation with Trixi.jl! Now, we can plot\n",
    "the solution at the final time using Plots.jl."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "You can already see that discontinuities will develop and oscillations start to\n",
    "occur around steep parts of the wave. That's expected from our central discretization.\n",
    "To avoid these issues, we need to use dissipative numerical fluxes (approximate\n",
    "Riemann solvers) at interfaces."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Advanced setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Thus, we add a Godunov's flux for our cubic equation. That is easy for this equation\n",
    "since the wave speed `f'(u) = 3u^2` is always non-negative."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@inline Trixi.flux_godunov(u_ll, u_rr, orientation, equation::CubicEquation) = flux(u_ll, orientation, equation)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Let's run the example again but with a dissipative numerical flux at interfaces.\n",
    "`remake` will recreate the semidiscretization we used before and only change\n",
    "selected parameters, in this case the `solver`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "# A new setup with dissipation\n",
    "semi = remake(semi, solver=DGSEM(3, flux_godunov))\n",
    "ode = semidiscretize(semi, tspan)\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...)\n",
    "plot!(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "You can see that there are fewer oscillations, in particular around steep edges.\n",
    "Now let's increase the final time (and also the spatial resolution)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "# A larger final time: Nonclassical shocks develop (you can even increase the refinement to 12)\n",
    "semi = remake(semi, mesh=TreeMesh(-1.0, 1.0, initial_refinement_level=8, n_cells_max=10^5))\n",
    "ode = semidiscretize(semi, (0.0, 0.5) #= tspan =#)\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...)\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "You can observe that nonclassical shocks develop and are stable under grid refinement,\n",
    "e.g. for `initial_refinement_level=12`. In this case, these nonclassical shocks\n",
    "can be avoided by using an entropy-dissipative semidiscretization. Thus, we need\n",
    "to define an entropy-conservative numerical flux"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "@inline function Trixi.flux_ec(u_ll, u_rr, orientation, equation::CubicEquation)\n",
    "  return SVector(0.25 * (u_ll[1]^3 + u_ll[1]^2 * u_rr[1] + u_ll[1] * u_rr[1]^2 + u_rr[1]^3))\n",
    "end"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "and use a `VolumeIntegralFluxDifferencing` instead of the standard\n",
    "`VolumeIntegralWeakForm` in the DGSEM."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "# Let's use a provably entropy-dissipative semidiscretization\n",
    "semi = remake(semi, solver=DGSEM(3, flux_godunov, VolumeIntegralFluxDifferencing(flux_ec)))\n",
    "ode = semidiscretize(semi, (0.0, 0.5))\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...);\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Possible next steps could be\n",
    "- to define `Trixi.max_abs_speeds(u, equations::CubicEquation) = 3 * u[1]^2`\n",
    "  to use CFL-based time step control via a `StepsizeCallback`\n",
    "- to define quantities of interest like `Trixi.entropy(u, equations::CubicEquation) = u[1]^2`\n",
    "  and integrate them in a simulation using the `AnalysisCallback`\n",
    "- to experiment with shock-capturing volume integrals `VolumeIntegralShockCapturingHG`\n",
    "  and adaptive mesh refinement `AMRCallback`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For further reading, Trixi.jl provides another example on adding a scalar equation. In the\n",
    "[elixir about the KPP problem](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_2d_dgsem/elixir_kpp.jl),\n",
    "the 2D scalar \"KPP equation\" from [Kurganov, Petrova, Popov (2007)](https://doi.org/10.1137/040614189) is\n",
    "implemented."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Summary of the code"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "To sum up, here is the complete code that we used (without the callbacks since these create a\n",
    "lot of unnecessary output in the doctests of this tutorial).\n",
    "In addition, we create the `struct` inside the new module `CubicConservationLaw`. That\n",
    "ensures that we can re-create `struct`s defined therein without having to restart Julia."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "# Define new physics\n",
    "module CubicConservationLaw\n",
    "\n",
    "using Trixi\n",
    "\n",
    "struct CubicEquation <: Trixi.AbstractEquations{1 #= number of spatial dimensions =#,\n",
    "                                                1 #= number of primary variables, i.e. scalar =#}\n",
    "end\n",
    "\n",
    "@inline Trixi.flux(u, orientation, equation::CubicEquation) = u.^3\n",
    "Trixi.varnames(_, ::CubicEquation) = (\"scalar\",)\n",
    "\n",
    "@inline Trixi.flux_godunov(u_ll, u_rr, orientation, equation::CubicEquation) = flux(u_ll, orientation, equation)\n",
    "@inline function Trixi.flux_ec(u_ll, u_rr, orientation, equation::CubicEquation)\n",
    "  return SVector(0.25 * (u_ll[1]^3 + u_ll[1]^2 * u_rr[1] + u_ll[1] * u_rr[1]^2 + u_rr[1]^3))\n",
    "end\n",
    "\n",
    "end # module\n",
    "\n",
    "\n",
    "# Create a simulation setup\n",
    "import .CubicConservationLaw\n",
    "using Trixi\n",
    "using OrdinaryDiffEq\n",
    "using Plots\n",
    "\n",
    "equation = CubicConservationLaw.CubicEquation()\n",
    "\n",
    "initial_condition_sine(x, t, equation::CubicConservationLaw.CubicEquation) = SVector(sinpi(x[1]))\n",
    "\n",
    "mesh = TreeMesh(-1.0, 1.0, # min/max coordinates\n",
    "                initial_refinement_level=4,\n",
    "                n_cells_max=10^4)\n",
    "\n",
    "solver = DGSEM(3 #= polynomial degree =#, flux_central)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equation, initial_condition_sine, solver)\n",
    "\n",
    "# Create ODE problem with given time span\n",
    "tspan = (0.0, 0.1)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...)\n",
    "plot(sol)\n",
    "\n",
    "\n",
    "# A new setup with dissipation\n",
    "semi = remake(semi, solver=DGSEM(3, flux_godunov))\n",
    "ode = semidiscretize(semi, tspan)\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...)\n",
    "plot!(sol)\n",
    "\n",
    "\n",
    "# A larger final time: Nonclassical shocks develop (you can even increase the refinement to 12)\n",
    "semi = remake(semi, mesh=TreeMesh(-1.0, 1.0, initial_refinement_level=8, n_cells_max=10^5))\n",
    "ode = semidiscretize(semi, (0.0, 0.5))\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...)\n",
    "plot(sol)\n",
    "\n",
    "\n",
    "# Let's use a provably entropy-dissipative semidiscretization\n",
    "semi = remake(semi, solver=DGSEM(3, flux_godunov, VolumeIntegralFluxDifferencing(flux_ec)))\n",
    "ode = semidiscretize(semi, (0.0, 0.5))\n",
    "sol = solve(ode, SSPRK43(); ode_default_options()...)\n",
    "plot(sol)"
   ],
   "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.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}