{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 2: Behind the scenes of a simulation setup"
   ],
   "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 will guide you through a simple Trixi.jl setup (\"elixir\"), giving an overview of\n",
    "what happens in the background during the initialization of a simulation. While the setup\n",
    "described herein does not cover all details, it involves relatively stable parts of Trixi.jl that\n",
    "are unlikely to undergo significant changes in the near future. The goal is to clarify some of\n",
    "the more fundamental, *technical* concepts that are applicable to a variety of\n",
    "(also more complex) configurations."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Trixi.jl follows the [method of lines](http://www.scholarpedia.org/article/Method_of_lines) concept for solving partial differential equations (PDEs).\n",
    "Firstly, the PDEs are reduced to a (potentially huge) system of\n",
    "ordinary differential equations (ODEs) by discretizing the spatial derivatives. Subsequently,\n",
    "these generated ODEs may be solved with methods available in OrdinaryDiffEq.jl or those specifically\n",
    "implemented in Trixi.jl. The following steps elucidate the process of transitioning from PDEs to\n",
    "ODEs within the framework of Trixi.jl."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Basic setup"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Import essential libraries and specify an equation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi, OrdinaryDiffEq\n",
    "equations = LinearScalarAdvectionEquation2D((-0.2, 0.7))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Generate a spatial discretization using a `TreeMesh` with a pre-coarsened set of cells."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "coordinates_min = (-2.0, -2.0)\n",
    "coordinates_max = (2.0, 2.0)\n",
    "\n",
    "coarsening_patches = ((type = \"box\", coordinates_min = [0.0, -2.0],\n",
    "                       coordinates_max = [2.0, 0.0]),)\n",
    "\n",
    "mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 2,\n",
    "                n_cells_max = 30_000,\n",
    "                coarsening_patches = coarsening_patches)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The created `TreeMesh` looks like the following:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![TreeMesh_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/d5ef76ee-8246-4730-a692-b472c06063a3)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Instantiate a `DGSEM` solver with a user-specified polynomial degree. The solver\n",
    "will define `polydeg + 1` [Gauss-Lobatto nodes](https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss%E2%80%93Lobatto_rules) and their associated weights within\n",
    "the reference interval $[-1, 1]$ in each spatial direction. These nodes will be subsequently\n",
    "used to approximate solutions on each leaf cell of the `TreeMesh`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(polydeg = 3)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Gauss-Lobatto nodes with `polydeg = 3`:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![Gauss-Lobatto_nodes_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/1d894611-801e-4f75-bff0-d77ca1c672e5)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Overview of the `SemidiscretizationHyperbolic` type"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "At this stage, all necessary components for configuring the spatial discretization are in place.\n",
    "The remaining task is to combine these components into a single structure that will be used\n",
    "throughout the entire simulation process. This is where `SemidiscretizationHyperbolic`\n",
    "comes into play."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition_convergence_test,\n",
    "                                    solver)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The constructor for the `SemidiscretizationHyperbolic` object calls numerous sub-functions to\n",
    "perform the necessary initialization steps. A brief description of the key sub-functions is\n",
    "provided below."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `init_elements(leaf_cell_ids, mesh, equations, dg.basis, RealT, uEltype)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  The fundamental elements for approximating the solution are the leaf\n",
    "  cells. The solution is constructed as a polynomial of the degree specified in the `DGSEM`\n",
    "  solver in each spatial direction on each leaf cell. This polynomial approximation is evaluated\n",
    "  at the Gauss-Lobatto nodes mentioned earlier. The `init_elements` function extracts\n",
    "  these leaf cells from the `TreeMesh`, assigns them the label \"elements\", records their\n",
    "  coordinates, and maps the Gauss-Lobatto nodes from the 1D interval $[-1, 1]$ onto each coordinate axis\n",
    "  of every element."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  ![elements_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/9f486670-b579-4e42-8697-439540c8bbb4)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  The visualization of elements with nodes shown here includes spaces between elements, which do\n",
    "  not exist in reality. This spacing is included only for illustrative purposes to underscore the\n",
    "  separation of elements and the independent projection of nodes onto each element."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `init_interfaces(leaf_cell_ids, mesh, elements)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  At this point, the elements with nodes have been defined; however, they lack the necessary\n",
    "  communication functionality. This is crucial because the local solution polynomials on the\n",
    "  elements are not independent of each other. Furthermore, nodes on the boundary of adjacent\n",
    "  elements share the same spatial location, which requires a method to combine this into a\n",
    "  meaningful solution.\n",
    "  Here [Riemann solvers](https://en.wikipedia.org/wiki/Riemann_solver#Approximate_solvers)\n",
    "  come into play which can handle the principal ambiguity of a multi-valued solution at the\n",
    "  same spatial location."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  As demonstrated earlier, the elements can have varying sizes. Let us initially consider\n",
    "  neighbors with equal size. For these elements, the `init_interfaces` function generates\n",
    "  interfaces that store information about adjacent elements, their relative positions, and\n",
    "  allocate containers for sharing solution data between neighbors during the solution process."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  In our visualization, these interfaces would conceptually resemble tubes connecting the\n",
    "  corresponding elements."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  ![interfaces_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/bc3b6b02-afbc-4371-aaf7-c7bdc5a6c540)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `init_mortars(leaf_cell_ids, mesh, elements, dg.mortar)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  Returning to the consideration of different sizes among adjacent elements, within the\n",
    "  `TreeMesh`, adjacent leaf cells can vary in side length by a maximum factor of two. This\n",
    "  implies that a large element has one neighbor of\n",
    "  equal size with a connection through an interface, or two neighbors at half the size,\n",
    "  requiring a connection through so called \"mortars\". In 3D, a large element would have\n",
    "  four small neighbor elements."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  Mortars store information about the connected elements, their relative positions, and allocate\n",
    "  containers for storing the solutions along the boundaries between these elements."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  Due to the differing sizes of adjacent elements, it is not feasible to directly map boundary\n",
    "  nodes of adjacent elements. Therefore, the concept of mortars employs a mass-conserving\n",
    "  interpolation function to map boundary nodes from a larger element to a smaller one."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  In our visualization, mortars are represented as branched tubes."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  ![mortars_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/43a95a60-3a31-4b1f-8724-14049e7a0481)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `init_boundaries(leaf_cell_ids, mesh, elements)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  In order to apply boundary conditions, it is necessary to identify the locations of the\n",
    "  boundaries. Therefore, we initialize a \"boundaries\" object, which records the elements that\n",
    "  contain boundaries, specifies which side of an element is a boundary, stores the coordinates\n",
    "  of boundary nodes, and allocates containers for managing solutions at these boundaries."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  In our visualization, boundaries and their corresponding nodes are highlighted with green,\n",
    "  semi-transparent lines."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  ![boundaries_example](https://github.com/trixi-framework/Trixi.jl/assets/119304909/21996b20-4a22-4dfb-b16a-e2c22c2f29fe)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "All the structures mentioned earlier are collected as a cache of type `NamedTuple`. Subsequently,\n",
    "an object of type `SemidiscretizationHyperbolic` is initialized using this cache, initial and\n",
    "boundary conditions, equations, mesh and solver."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In conclusion, the primary purpose of a `SemidiscretizationHyperbolic` is to collect equations,\n",
    "the geometric representation of the domain, and approximation instructions, creating specialized\n",
    "structures to interconnect these components in a manner that enables their utilization for\n",
    "the numerical solution of partial differential equations (PDEs)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As evident from the earlier description of `SemidiscretizationHyperbolic`, it comprises numerous\n",
    "functions called subsequently. Without delving into details, the structure of the primary calls\n",
    "are illustrated as follows:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![SemidiscretizationHyperbolic_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/8bf59422-0537-4d7a-9f13-d9b2253c19d7)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Overview of the `semidiscretize` function"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "At this stage, we have defined the equations and configured the domain's discretization. The\n",
    "final step before solving is to select a suitable time span and apply the corresponding initial\n",
    "conditions, which are already stored in the initialized `SemidiscretizationHyperbolic` object."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The purpose of the `semidiscretize` function is to wrap the semidiscretization as an\n",
    "`ODEProblem` within the specified time interval. During this procedure the approximate solution\n",
    " is created at the given initial time via the specified `initial_condition` function from the\n",
    " `SemidiscretizationHyperbolic` object. This `ODEProblem` can be subsequently passed to the\n",
    "`solve` function from the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package\n",
    "or to `Trixi.solve`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "ode = semidiscretize(semi, (0.0, 1.0));"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The `semidiscretize` function involves a deep tree of subsequent calls, with the primary ones\n",
    "explained below."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `allocate_coefficients(mesh, equations, solver, cache)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  To apply initial conditions, a data structure (\"container\") needs to be generated to store the\n",
    "  initial values of the target variables for each node within each element."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  Since only one-dimensional `Array`s are `resize!`able in Julia, we use `Vector`s as an internal\n",
    "  storage for the target variables and `resize!` them whenever needed, e.g. to change the number\n",
    "  of elements. Then, during the solving process the same memory is reused by `unsafe_wrap`ping\n",
    "  multi-dimensional `Array`s around the internal storage."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `wrap_array(u_ode, semi)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  As previously noted, `u_ode` is constructed as a 1D vector to ensure compatibility with\n",
    "  OrdinaryDiffEq.jl. However, for internal use within Trixi.jl, identifying which part of the\n",
    "  vector relates to specific variables, elements, or nodes can be challenging."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  This is why the `u_ode` vector is wrapped by the `wrap_array` function using `unsafe_wrap`\n",
    "  to form a multidimensional array `u`. In this array, the first dimension corresponds to\n",
    "  variables, followed by N dimensions corresponding to nodes for each of N space dimensions.\n",
    "  The last dimension corresponds to the elements.\n",
    "  Consequently, navigation within this multidimensional array becomes noticeably easier."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  \"Wrapping\" in this context involves the creation of a reference to the same storage location\n",
    "  but with an alternative structural representation. This approach enables the use of both\n",
    "  instances `u` and `u_ode` as needed, so that changes are simultaneously reflected in both.\n",
    "  This is possible because, from a storage perspective, they share the same stored data, while\n",
    "  access to this data is provided in different ways."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "- `compute_coefficients!(u, initial_conditions, t, mesh::DG, equations, solver, cache)`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "  Now the variable `u`, intended to store solutions, has been allocated and wrapped, it is time\n",
    "  to apply the initial conditions. The `compute_coefficients!` function calculates the initial\n",
    "  conditions for each variable at every node within each element and properly stores them in the\n",
    "  `u` array."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "At this stage, the `semidiscretize` function has all the necessary components to initialize and\n",
    "return an `ODEProblem` object, which will be used by the `solve` function to compute the\n",
    "solution."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In summary, the internal workings of `semidiscretize` with brief descriptions can be presented\n",
    "as follows."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![semidiscretize_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/491eddc4-aadb-4e29-8c76-a7c821d0674e)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Functions `solve` and `rhs!`"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Once the `ODEProblem` object is initialized, the `solve` function and one of the ODE solvers from\n",
    "the OrdinaryDiffEq.jl package can be utilized to compute an approximated solution using the\n",
    "instructions contained in the `ODEProblem` object."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), dt = 0.01,\n",
    "            save_everystep = false);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Since the `solve` function and the ODE solver have no knowledge\n",
    "of a particular spatial discretization, it is necessary to define a\n",
    "\"right-hand-side function\", `rhs!`, within Trixi.jl."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Trixi.jl includes a set of `rhs!` functions designed to compute `du`, i.e.,\n",
    "$\\frac{\\partial u}{\\partial t}$ according to the structure\n",
    "of the setup. These `rhs!` functions calculate interface, mortars, and boundary fluxes, in\n",
    "addition to surface and volume integrals, in order to construct the `du` vector. This `du` vector\n",
    "is then used by the time integration method to obtain the solution at the subsequent time step.\n",
    "The `rhs!` function is called by time integration methods in each iteration of the solve loop\n",
    "within OrdinaryDiffEq.jl, with arguments `du`, `u`, `semidiscretization`, and the current time."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Trixi.jl uses a two-levels approach for `rhs!` functions. The first level is limited to a\n",
    "single function for each `semidiscretization` type, and its role is to redirect data to the\n",
    "target `rhs!` for specific solver and mesh types. This target `rhs!` function is responsible\n",
    "for calculating `du`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Path from the `solve` function call to the appropriate `rhs!` function call:"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![rhs_structure](https://github.com/trixi-framework/Trixi.jl/assets/119304909/dbea9a0e-25a4-4afa-855e-01f1ad619982)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Computed solution:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(sol)\n",
    "pd = PlotData2D(sol)\n",
    "plot!(getmesh(pd))"
   ],
   "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
}