## 3 Theoretical Baseline
### Discontinuous Galerkin Method  
This method was originally proposed in 1973 to solve challenging hyperbolic transport equations in nuclear physics [25]. In spirit, it can be held as a synthesis of finite element and finite volume schemes and poses a generalized variant of both. To derive such a scheme, we begin by stating the strong form of a given PDE. The most straightforward example in this case would be a first-order linear advection equation with a homogeneous von Neumann boundary condition:  
\[\partial_{t}\alpha+\underline{u}\cdot\
abla\alpha=0,\qquad\frac{\partial \alpha}{\partial n}=0\quad\forall x\in\partial\Omega, \tag{1}\]  
where \(\partial_{t}\) signifies the temporal derivative, \(x\) is the set of spatial coordinates, \(\partial\Omega\) denotes the boundary of the computational domain \(\Omega\), and \(n\) is the unit normal with respect to \(\partial\Omega\). In this case, and henceforth in this article, we assume for reasons of simplicity that the velocity field is divergence-free, i.e., \(\
abla\cdot\underline{u}=0\). We now state the weak form of Equation (1), that is, we multiply with a test function \(v\), integrate over the entirety of the domain \(\Omega\), and apply partial integration to the second term on the left-hand side that contains the nabla operator. For a divergence-free velocity field, one may set \(\underline{u}\cdot\
abla\alpha=\
abla\cdot(\underline{u}\alpha)\), which results in the following formulation: Find \(\alpha\in V\) such that  
\[\int_{\Omega}v\cdot\partial_{t}\alpha\,dx+\int_{\partial\Omega}v\alpha( \underline{u}\cdot\underline{u})\,ds-\int_{\Omega}\
abla v\cdot\underline{u }\alpha\,dx=0,\qquad\forall v\in W, \tag{2}\]  
where we need to make an appropriate choice for the solution space \(V\) and the test space \(W\) which may, but do not need to, differ from each other. Due to partial integration, we now encounter an additional term that has to be integrated over the domain boundary \(\partial\Omega\), where \((\underline{u}\cdot\underline{n})\) denotes the velocity component normal to the boundary. To make such a problem solvable by a computer, one must additionally choose the discretization of the solution space \(V\), denoted \(V_{h}\). A particularly popular choice of space is the set of Lagrange polynomials. In addition, the physical space must be discretized in the form of a triangulation. The DG scheme then consists of assembling the finite-dimensional, linear system on the element level. This enables high locality of the solution process, which leads to efficient computation on parallel architectures, as less data transfer is required. One resulting key feature of the DG scheme is that the elements now do not overlap anymore in terms of their degrees of freedom. Thus, the global problem is broken up into individual problems. This, in general, leads to large systems that are, however, sparse and, in the case of the mass matrix, even block diagonal. The remaining term, often denoted the numerical flux, is the only term within the physical domain that ensures coupling across elements. Through evaluation of this surface integral, adjacent degrees of freedom are coupled and, thus, global conservation of quantities can be assured. As the polynomial space of DG schemes only belongs to the \(L^{2}\) space of functions but not \(H^{1}\), the basis functions are discontinuous; thus, the derivative at boundaries is not well defined. Solving PDEs involving second derivatives is, thus, not possible as is. As aconsequence, there have been many successful extensions of this method to circumvent that problem. At this point, we name the most prevalent schemes, namely the symmetric interior penalty [26], hybridizable [27] and local DG scheme [28]. These methods, despite having different approaches, have been extensively studied and compared to each other [29]. As it turns out, all methods work well and have individual advantages and disadvantages. For this work, we will take the hybridizable DG scheme as a general framework.