# Modal DG Formulation for the Burgers' Equation

## Notebook setup


```python
import numpy
from matplotlib import pyplot
%matplotlib inline
```

## Theoretical aspects and DG-FEM formulation

We consider the Burgers' inviscid equation (IB) described by the following conservation law,

$$
\frac{\partial u(x,t)}{\partial t} +  \frac{\partial f(u(x,t))}{\partial x} = 0, \ x \in [L,R] = \Omega
$$

with initial conditions

$$
u(x,0) = u_{0}(x) 
$$

and boundary conditions

$$
\begin{align}
u(0,t) &= g_0(t) \\
u(L,t) &= g_1(t)
\end{align}
$$

where $u(x,t)$ is a species concentraion, $x \in \mathbb{R}$ and $t \in \mathbb{R}_{+}$ are the space and time coordinates, and $\Omega$ is the problem domain, respectively. The typical quadratic flux function $f: \mathbb{R} \rightarrow \mathbb{R}$ used in this model can be written as

$$
f(u) = \frac{1}{2} u^2
$$

## The Weak Local Statement

In order to obtain an approximation for the solution of the previous equation, we will rewrite the previous initial-boundary value problem in an integral form. For this, we assume that our integral form will be valid for each individual element of a given mesh and the global approximate solution will be enforced by a mechanism called numerical flux which allows inter element communication through their boundaries. Using this premise, for each element $\Omega^{k}$ of the mesh, we form the local residual as follows

\begin{equation}
	x \in \Omega^{k}: \mathcal{R}_{h}(x,t) = \frac{\partial u^{k}_{h}}{\partial t} + \frac{\partial f^{k}_{h}}{\partial x}
\end{equation}

where $u^{k}_{h}$ and $f^{k}_{h}$ are numerical approximations to the solution and the flux, respectively.

As in the Galerkin finite element formulation, we require that the residual is orthogonal to all test functions $\phi_{h}$ belonging to a finite dimension vector space $V_{h}$, leading to

\begin{equation}
	\int_{\Omega^{k}} \mathcal{R}_{h}(x,t)\phi^{k}_{j}(x) \: dx = 0,
\end{equation}

for all the elemental test functions $\phi^{k}_{j}(x)$.

Substituting the residual expression into the orthogonal relation yields

\begin{equation}
	\int_{\Omega^{k}} \frac{\partial u^{k}_{h}}{\partial t}\phi^{k}_{j}(x) + \frac{\partial f^{k}_{h}}{\partial x}\phi^{k}_{j}(x) \: dx = 0.
\end{equation}

Integrating the second term of the latter equation by parts once gives us the weak local statement of the Burgers' equation

\begin{equation}
	\int_{\Omega^{k}} \left(\frac{\partial u^{k}_{h}}{\partial t}\phi^{k}_{j}(x) - f^{k}_{h}\frac{d\phi^{k}_{j}}{dx}\right) \: dx = -\left[f^{k}_{h}\phi^{k}_{j}\right]^{x^{k}_{r}}_{x^{k}_{l}}.
\end{equation}

### Discontinuous Galerkin-FEM Discretization

We approximate $\Omega$ by $K$ nonoverlapping elements, $x \in [x^{k}_{l},x^{k}_{r}] = \Omega^{k}$ that we call a _finite element mesh_. On each of these elements we will assume that the local solution, $u^{k}_{h}(x,t)$, __and__ the local flux, $f^{k}_{h}(u_{h}(x,t))$, can be expressed as a polynomial of order $N$

$$
\begin{aligned}
x \in \Omega^{k}: \ u^{k}_{h}(x,t) &= \sum^{N+1}_{j=1}\hat{u}^{k}_{j}(t)\psi^{k}_{j}(x), \\
f^{k}_{h}(u_{h}(x,t)) &= \sum^{N+1}_{j=1} \hat{f}^{k}_{j}(t)\psi^{k}_{j}(x).
\end{aligned}
$$

where $\psi^{k}_{j}(x)$ are the local basis functions which can be chosen as of the modal or nodal type. In this work we consider the Legendre polynomials, $P_{n}(\xi)$, defined on the reference interval $[-1,1]$, for the modal basis and the Lagrange polynomials, $\ell^{k}_{i}(x)$, defined through the $N_{p} = N+1$ Gauss-Legendre-Lobatto quadrature points, for the nodal case. These two basis functions representations are connected by the definition of the expansion coefficients, $\hat{u}_{i}$, as follows

$$
{\cal V}\hat{\mathbb{u}} = \mathbb{u}
$$

where ${\cal V}_{ij} = P_{j-1}(\xi_{i})$, $\hat{\mathbb{u}}_{i} = \hat{u}_{i}$, $\mathbb{u}_{i} = u(\xi_{i})$. The matrix ${\cal V}$ is recognized as the generalized Vandermond matrix.

We now return to the weak local statement and substitute the function $u$ and the flux $f$ by their discrete counterparts. We will also replace the arbitrary test function $\phi^{k}_{j}$ by the piecewise polynomial functions $\psi^{k}_{i}(x) \in V_{h}$. We also introduce the numerical flux $f^{*}$ in order to recover the global meaning of the discontinuous Galerkin (DG-FEM) approximate solution of the Burgers' equation. Now, the semi-discrete form of the DG-FEM formulation can be written as

$$
\int_{\Omega^{k}} \left(\frac{\partial u^{k}_{h}}{\partial t} \psi^{k}_{i}(x) - f^{k}_{h}(u^{k}_{h}) \frac{d\psi^{k}_{i}(x)}{dx}\right) dx = -\left[\psi^{k}_{i}(x)f^{*}\right]^{x^{k}_{r}}_{x^{k}_{l}}, \ 1 \leq i \leq N+1.
$$

Substituting the expressions for $u_h^k(x,y)$ and $f_h^k(u_h(x,t))$ into the last equation and using an affine transformation to perform integrations on the reference element $\hat{\Omega} = [-1,1]$, we can rewrite the semidiscrete DG-FEM weak formulation in its matrix form,

$$
\mathcal{M}^{k}\frac{d}{dt}\boldsymbol{\hat{u}}^{k} - \mathcal{S}^{T}\boldsymbol{\hat{f}}^{k} = -\boldsymbol{\hat{\psi}}(\xi_{r})f^{*}(x^{k}_{r},x^{k+1}_{l})+\boldsymbol{\hat{\psi}}(\xi_{l})f^{*}(x^{k}_{l},x^{k-1}_{r}),
$$

where, for $1 \leq i,j \leq N+1$,

$$
\begin{align}
\mathcal{M}^{k} &= \frac{h^{k}}{2}\int_{\hat{\Omega}} \hat{\psi}_{i}\hat{\psi}_{j} \: d\xi, \\
\mathcal{S}^{T} &= \int_{\hat{\Omega}} \hat{\psi}_{i}\frac{d\hat{\psi}_{j}}{d\xi} \: d\xi,
\end{align}
$$

are the local mass and stiffness matrices, respectively. The term $h^{k}/{2}$ is the jacobian associated to the affine transformation and its value is half of the size of the $k$-th element, $\Omega^{k}$. Furthermore, we have

$$
\begin{align}
\boldsymbol{\hat{u}}^{k} &= [\hat{u}^{k}_{1},\ldots,\hat{u}^{k}_{N+1}]^{T}, \\
\boldsymbol{\hat{f}}^{k} &= [\hat{f}^{k}_{1},\ldots,\hat{f}^{k}_{N+1}]^{T}, \\
\boldsymbol{\hat{\psi}} &= [\hat{\psi}_{1}(\xi),\ldots,\hat{\psi}_{N+1}(\xi)]^{T}.
\end{align}
$$

as the vectors of local coefficients, flux coefficients and local test functions, respectively.

Before proceeding with the numerical flux definition, we have to find a way to represent the flux coefficients $\hat{f}^{k}_{i}$ in terms of the approximate solution $u^{k}_{h}(x,t)$. To keep things simple, we propose a three step algorithm to approximate these coefficients:

1. Compute the solution, $u^{k}_{h}(x,t)$, at $N+1$ quadrature points

$$
u^{k}_{h}(x_{l},t) = \sum^{N+1}_{i=1}\hat{u}_{i}(t)\hat{\psi}_{i}(\xi_{l}), \ 1 \leq l \leq N+1.
$$

2. Compute the flux, $f^{k}_{h}(x,t)$, at $N+1$ quadrature points

$$
f^{k}_{h}(x_{l},t) = \frac{u^{k}_{h}(x_{l},t)^{2}}{2}, \ 1 \leq l \leq N+1.
$$

3. Project the nonlinear flux evaluated at the physical points $x_{l}$ into the polynomial space

$$
\begin{align*}
&f^{k}_{h}(x,t) = \sum^{N+1}_{j=1}\hat{f}^{k}_{j}(t)\psi_{j}(x) \\
&\sum^{N+1}_{j=1}\hat{f}^{k}_{j}(t)\int_{\Omega^{k}}\psi_{i}(x)\psi_{j}(x) \: dx = \int_{\Omega^{k}}f^{k}_{h}(x,t)\psi_{i}(x) \: dx \\
&\sum^{N+1}_{j=1}\hat{f}^{k}_{j}(t)\frac{h^{k}}{2}\int_{\hat{\Omega}}\hat{\psi}_{i}(\xi)\hat{\psi}_{j}(\xi) \: d\xi = \frac{h^{k}}{2}\int_{\hat{\Omega}}f^{k}_{h}(x(\xi),t)\hat{\psi}_{i}(\xi) \: d\xi \\
&\mathcal{M}\boldsymbol{\hat{f}}^{k} = \boldsymbol{b}^{k} \Rightarrow \boldsymbol{\hat{f}}^{k} = \mathcal{M}^{-1}\boldsymbol{b}^{k}.
\end{align*}
$$

where
$$
\begin{align*}
&\mathcal{M} = \int_{\hat{\Omega}} \hat{\psi}_{i}\hat{\psi}_{j} \: d\xi, \\
&\boldsymbol{b}^{k} = \int_{\hat{\Omega}}f^{k}_{h}(x(\xi),t)\hat{\psi}_{i}(\xi) \: d\xi,
\end{align*}
$$

are the local mass matrix and the local load vector defined on the reference element $\hat{\Omega}$, respectively.

We have used $N+1$ quadrature points to integrate the flux, $f^{k}_{h}(x,t)$. As a consequence of the quadratic non linearity of the flux, the order of the polynomial to be integrated in the local load vector corresponds to the sum of the polynomial orders of the flux, $2N$, and the basis function, $N$, resulting in a total polynomial order of $3N$. This requires a minimum of $(3N+1)/2$ quadrature points to exactly compute the vector $\boldsymbol{b}^{k}$. The solution evaluation at $N+1$ quadrature points in step one does not suffice to uniquely specify a $3N$-th-order polynomial and an aliasing error is introduced in the numerical approximation. To remedy this kind of error, we can evaluate the nonlinear flux at $(3N+1)/2$ quadrature points in order to consistently integrate the coefficients of the vector $\boldsymbol{b}^{k}$. This strategy results in more stability for the whole numerical scheme.

Now, we are ready to define the numerical flux, $f^{*}$, that appears in the equation above. We consider the monotone Lax-Friedrichs flux in order to stabilize the DG-FEM discretization

$$
f^{*}(u^{-}_{h},u^{+}_{h}) = \{\{f_{h}(u_{h})\}\}+\frac{C}{2}[[u_{h}]],
$$

where $\{\{.\}\}$ and $[[.]]$ are the average and jump operators usually defined as follows

$$
\{\{u\}\} = \frac{u^{-}+u^{+}}{2}
$$

where $u$ can be both a scalar or a vector quantity and

$$
\begin{align}
[[u]] &= \hat{\mathbf{n}}^{-}u^{-} + \hat{\mathbf{n}}^{+}u^{+}, \\
[[\mathbf{u}]] &= \hat{\mathbf{n}}^{-}\cdot\mathbf{u}^{-} + \hat{\mathbf{n}}^{+}\cdot\mathbf{u}^{+}
\end{align}
$$

where $\hat{\mathbf{n}}$ is defined as the normal vector along the left or right boundaries of the $k$-th element. We observe that the jump operator is defined differently depending on whether $u$ is a scalar or a vector, $\mathbf{u}$. The superscripts "$-$" and "$+$" stand for the interior and exterior information of the element, respectively. 

To write the final semidiscrete form of the proposed DG-FEM scheme, we substitute the definition of the numerical flux into the right hand side of the equation in matricial form, as follows

$$
-\boldsymbol{\hat{\psi}}(1)f^{*}(x^{k}_{r},x^{k+1}_{l})+\boldsymbol{\hat{\psi}}(-1)f^{*}(x^{k}_{l},x^{k-1}_{r}) =
$$

$$
\begin{align*}
&= \left[-\boldsymbol{\hat{\psi}}(1)\left(\{\{f^{k}_{h}\}\} + \frac{C}{2}[[u^{k}_{h}]]\right) + \boldsymbol{\hat{\psi}}(-1)\left(\{\{f^{k}_{h}\}\} + \frac{C}{2}[[u^{k}_{h}]]\right)\right] \\
&= \left[-\boldsymbol{\hat{\psi}}(1)\left(\frac{1}{2}(f^{-}_{h}+f^{+}_{h})+\frac{C}{2}(\hat{\mathbf{n}}^{-}_{r}u^{-}_{h} + \hat{\mathbf{n}}^{+}_{r}u^{+}_{h})\right)\right. \\
&+\left.\boldsymbol{\hat{\psi}}(-1)\left(\frac{1}{2}(f^{-}_{h}+f^{+}_{h})+\frac{C}{2}(\hat{\mathbf{n}}^{-}_{l}u^{-}_{h} + \hat{\mathbf{n}}^{+}_{l}u^{+}_{h})\right)\right]
\end{align*}
$$

where the subscripts $r$ or $l$ for the normal vector $\hat{\mathbf{n}}$ indicates that we are taking the normal vector at the right or left side of the $k$ element, respectively. We also notice that the local coordinates at the right and left boundaries of the reference element assume the values $\xi_{r} = 1$ and $\xi_{l} = -1$.

With a bit of elementary calculus, we can easily estimate the constant $C$ in the case of the Burgers's equation
\begin{align}
	&f_{u} = \frac{\partial f}{\partial u} = \frac{\partial}{\partial u}\left(\frac{u^{2}}{2}\right) = u, \nonumber \\
	&C = \max|f_{u}| = \max\left|u\right|
\end{align}

Substituting the expressions for the approximations $u^{k}_{h}$ and $f^{k}_{h}$ into the equation above and remembering that $\hat{\mathbf{n}}^{-} = -\hat{\mathbf{n}}^{+}$ yields

$$
-\boldsymbol{\hat{\psi}}(1)f^{*}(x^{k}_{r},x^{k+1}_{l})+\boldsymbol{\hat{\psi}}(-1)f^{*}(x^{k}_{l},x^{k-1}_{r}) =
$$

$$
\begin{align}
& \: \left[-\frac{1}{2}\boldsymbol{\hat{\psi}}(1)\left((\boldsymbol{\hat{f}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(1)+(\boldsymbol{\hat{f}}^{k+1}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)\right)\right. \\
&-\frac{C}{2}\hat{\mathbf{n}}^{-}_{r}\boldsymbol{\hat{\psi}}(1)\left((\boldsymbol{\hat{u}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(1)-(\boldsymbol{\hat{u}}^{k+1}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)\right) \\
&+\frac{1}{2}\boldsymbol{\hat{\psi}}(-1)\left((\boldsymbol{\hat{f}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)+(\boldsymbol{\hat{f}}^{k-1}_{h})^{T}\boldsymbol{\hat{\psi}}(1)\right) \\
&+\left.\frac{C}{2}\hat{\mathbf{n}}^{-}_{l}\boldsymbol{\hat{\psi}}(-1)\left((\boldsymbol{\hat{u}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(-1)-(\boldsymbol{\hat{u}}^{k-1}_{h})^{T}\boldsymbol{\hat{\psi}}(1)\right)\right] \\
\\
=& \: \left[-\frac{1}{2}\left([\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(1)]\: \boldsymbol{\hat{f}}^{k}_{h}+[\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(-1)]\: \boldsymbol{\hat{f}}^{k+1}_{h}\right)\right.\\
&-\frac{C}{2}\hat{\mathbf{n}}^{-}_{r}\left([\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\boldsymbol{\hat{u}}^{k}_{h}-[\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\boldsymbol{\hat{u}}^{k+1}_{h}\right) \\
&+\frac{1}{2}\left([\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\boldsymbol{\hat{f}}^{k}_{h}+[\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\boldsymbol{\hat{f}}^{k-1}_{h}\right) \\
&+\left.\frac{C}{2}\hat{\mathbf{n}}^{-}_{l}\left([\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(-1)]\:\boldsymbol{\hat{u}}^{k}_{h}-[\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(1)]\:\boldsymbol{\hat{u}}^{k-1}_{h}\right)\right]
\end{align}
$$

Tensor multiplication of vectors of elemental test functions gives rise to matrices which we will call flux matrices. They are defined as follows

$$
\begin{align*}
\mathcal{F}^{k}_{r} &= [\hat{\boldsymbol{\psi}}(1)\otimes\hat{\boldsymbol{\psi}}(1)], \\
\mathcal{F}^{k+1}_{l} &= [\hat{\boldsymbol{\psi}}(1)\otimes\hat{\boldsymbol{\psi}}(-1)], \\
\mathcal{F}^{k}_{l} &= [\hat{\boldsymbol{\psi}}(-1)\otimes\hat{\boldsymbol{\psi}}(-1)], \\
\mathcal{F}^{k-1}_{r} &= [\hat{\boldsymbol{\psi}}(-1)\otimes\hat{\boldsymbol{\psi}}(1)].
\end{align*}
$$

Now we use the flux matrices definitions and rearrange terms to write the final matrix form of the right hand side of equation

$$
-\boldsymbol{\hat{\psi}}(1)f^{*}(x^{k}_{r},x^{k+1}_{l})+\boldsymbol{\hat{\psi}}(-1)f^{*}(x^{k}_{l},x^{k-1}_{r}) = \:
$$

$$
\begin{align} 
\left[-\frac{1}{2}\mathcal{F}^{k+1}_{r}\left(\boldsymbol{\hat{f}}^{k+1}_{h}-(\hat{\mathbf{n}}^{-}_{r}\cdot C)\boldsymbol{\hat{u}}^{k+1}_{h}\right) \right. &- \frac{1}{2}\mathcal{F}^{k}_{r}\left(\boldsymbol{\hat{f}}^{k}_{h}+(\hat{\mathbf{n}}^{-}_{r}\cdot C)\boldsymbol{\hat{u}}^{k}_{h}\right) \nonumber \\
	+\frac{1}{2}\mathcal{F}^{k-1}_{l}\left(\boldsymbol{\hat{f}}^{k-1}_{h}-(\hat{\mathbf{n}}^{-}_{l}\cdot C)\boldsymbol{\hat{u}}^{k-1}_{h}\right) &+ \left.\frac{1}{2}\mathcal{F}^{k}_{l}\left(\boldsymbol{\hat{f}}^{k}_{h}+(\hat{\mathbf{n}}^{-}_{l}\cdot C)\boldsymbol{\hat{u}}^{k}_{h}\right)\right].
\end{align}
$$

Noting that the exterior normal vector, $\hat{\mathbf{n}}^{-}$, assumes only two values, $\hat{\mathbf{n}}^{-}_{r} = 1$ and $\hat{\mathbf{n}}^{-}_{l} = -1$, at the right and the left boundaries, respectively, of any $k$ element , the semidiscrete form of the proposed DG scheme for the (IB) equation can finally be written as

$$
\begin{align}
\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\boldsymbol{\hat{u}}^{k}_{h} - \mathcal{S}^{T}\boldsymbol{\hat{f}}^{k}_{h} = \: 
\left[-\frac{1}{2}\mathcal{F}^{k+1}_{r}\left(\boldsymbol{\hat{f}}^{k+1}_{h}-C\boldsymbol{\hat{u}}^{k+1}_{h}\right) \right. &- \frac{1}{2}\mathcal{F}^{k}_{r}\left(\boldsymbol{\hat{f}}^{k}_{h}+C\boldsymbol{\hat{u}}^{k}_{h}\right) \\
+\frac{1}{2}\mathcal{F}^{k-1}_{l}\left(\boldsymbol{\hat{f}}^{k-1}_{h}+C\boldsymbol{\hat{u}}^{k-1}_{h}\right) &+ \left.\frac{1}{2}\mathcal{F}^{k}_{l}\left(\boldsymbol{\hat{f}}^{k}_{h}-C\boldsymbol{\hat{u}}^{k}_{h}\right)\right].
\end{align}
$$

### Filter Stabilization

A popular strategy to stabilize discretizations when discontinuities and shocks develop in the solution which is the case when dealing with conservation laws is the use of filtering. This technique is an advantageous way of introducing the dissipation needed to preserve stability of the numerical solution. The filters are characterized by modifying the expansion coefficients of the solution leading to a global modification of the solution. In this work, filtered expansions are mathematically expressed as

$$
\mathcal{F}_{N}\boldsymbol{u}_{h}(x,t) = \sum^{N+1}_{i=1}\sigma_{E}\left(\frac{i-1}{N}\right)\hat{u}_{i}(t)P_{i-1}(x).
$$

where $\sigma_{E}(\eta)$ is the exponential filter function, $\hat{u}_{i}(t)$ are the coefficients of the modal representation of the approximated solution $u_{h}(x,t)$ and $P_{i-1}(x)$ are the hierarchical Legendre polynomials. 

The chosen exponential filter is defined as follows
$$
\sigma_{E}(\eta) = \left\{ 
\begin{array}{ll}
1, & 0 \leq \eta \leq \eta_{c} = \frac{N_{c}}{N} \\
\exp(-\alpha((\eta-\eta_{c})/(1-\eta_{c}))^{s}), & \eta_{c} \leq \eta \leq 1,
\end{array}
\right.
$$

where $\alpha$ is an arbitrary constant that sets the amount of introduced dissipation, $s$ is the order of the filter, $N$ is the polynomial order of the approximation and $N_{c}$ represents the cutoff below which the low order modes are left untouched. We usually choose $\alpha = -\log(\varepsilon_{M}) \approx 36$ with $\varepsilon_{M}$ being the machine precision in double precision. We also restrict $s$ to be an even value.

### DG-FEM Applied to Second Order PDEs
Although filter stabilization is a good cost-benefit technique to control Gibbs phenomenon in the vicinity of shock formation, it usually compromises the solution accuracy nearly the points of discontinuity. An alternative approach to control numerical wiggles around shocks which is commonly used in aerodynamic applications is the introduction of small quantities of artificial viscosity into the problem solution. 

In the context of conservation laws, this technique requires the addition of a second order diffusive term at the right hand side of the original equation (IB) as follows

$$
\frac{\partial u}{\partial t} +  \frac{\partial f(u)}{\partial x} = \frac{\partial}{\partial x}\varepsilon(x)\frac{\partial u}{\partial x},
\label{convec-diffusion-eq}
$$

where $\varepsilon(x)$ is the artificial viscosity introduced. This equation must be subjected to appropriate initial and boundary conditions.

In order to develop a suitable DG discretization for this second order PDE, we rewrite the last statement as a system of first order equations

$$
\begin{align}
&\frac{\partial u}{\partial t} +  \frac{\partial}{\partial x}(f(u)-\sqrt{\varepsilon}q) = 0 \\
&q = \sqrt{\varepsilon}\frac{\partial u}{\partial x}
\end{align}
$$

As before, we write approximations for the variables $u$, $f$ and $q$, and multiply both equations by arbitrary test functions $\phi^{k}_{j}$ and integrate them over each element $\Omega^{k}$ to get the local weak statement of the modified conservation law:

$$
\begin{align}
&\int_{\Omega^{k}} \frac{\partial u^{k}_{h}}{\partial t} \phi^{k}_{j} \: dx + \int_{\Omega^{k}} \frac{\partial f^{k}_{h}}{\partial x} \phi^{k}_{j} \: dx - \int_{\Omega^{k}}\frac{\partial}{\partial x}(\sqrt{\varepsilon}q^{k}_{h}) \phi^{k}_{j} \: dx = 0, \\
&\int_{\Omega^{k}} q^{k}_{h} \phi^{k}_{j} \: dx - \int_{\Omega^{k}} \sqrt{\varepsilon}\frac{\partial u^{k}_{h}}{\partial x} \phi^{k}_{j} \: dx = 0
\end{align}
$$

Integrating the second term in the equations above by parts once yields

$$
\begin{align}
&\int_{\Omega^{k}} \frac{\partial u^{k}_{h}}{\partial t} \phi^{k}_{j} \: dx - \int_{\Omega^{k}} f^{k}_{h}\frac{\partial \phi^{k}_{j}}{\partial x} \: dx + \int_{\Omega^{k}} (\sqrt{\varepsilon}q^{k}_{h})\frac{\partial \phi^{k}_{j}}{\partial x} \: dx = -\left[f^{k}_{h}\phi^{k}_{j}\right]^{x^{k}_{r}}_{x^{k}_{l}} + \left[(\sqrt{\varepsilon}q^{k}_{h})\phi^{k}_{j}\right]^{x^{k}_{r}}_{x^{k}_{l}}, \\
&\int_{\Omega^{k}} q^{k}_{h} \phi^{k}_{j} \: dx + \int_{\Omega^{k}} u^{k}_{h}\frac{\partial}{\partial x}(\sqrt{\varepsilon}\phi^{k}_{j}) \: dx = \left[(\sqrt{\varepsilon}u^{k}_{h})\phi^{k}_{j}\right]^{x^{k}_{r}}_{x^{k}_{l}}
\end{align}
$$

In order to discretize the later system of equations using the DG approach we assume

\begin{align}
q^{k}_{h}(x,t) &= \sum^{N+1}_{j=1}\hat{q}^{k}_{j}(t)\psi^{k}_{j}(x),
\end{align}

in exactly the same way we did for $u^{k}_{h}(x,t)$ and $f^{k}_{h}(u_{h}(x,t))$.

Putting all equations together, replacing the arbitrary test function $\phi^{k}_{j}$ by the piecewise polynomial functions $\psi^{k}_{i}(x)$ and introducing the numerical fluxes $f^{*}$, $(\sqrt{\varepsilon}q)^{*}$ and $(\sqrt{\varepsilon}u)^{*}$ at the element interfaces, we get the DG-FEM semi-discrete form of the modified traffic flow problem

\begin{align}
&\sum^{N+1}_{j=1}\int_{\Omega^{k}} \frac{d\hat{u}^{k}_{j}}{dt} \psi^{k}_{j}\psi^{k}_{i} \: dx - \sum^{N+1}_{j=1}\int_{\Omega^{k}} \hat{f}^{k}_{j} \psi^{k}_{j}\frac{d\psi^{k}_{i}}{dx} \: dx + \sum^{N+1}_{j=1}\int_{\Omega^{k}} \sqrt{\varepsilon} \: \hat{q}^{k}_{j} \psi^{k}_{j}\frac{d\psi^{k}_{i}}{dx} \: dx = -\left[\psi^{k}_{i}f^{*}\right]^{x^{k}_{r}}_{x^{k}_{l}} + \left[\psi^{k}_{i}(\sqrt{\varepsilon}q^{k}_{h})^{*}\right]^{x^{k}_{r}}_{x^{k}_{l}},& \\ 
&\sum^{N+1}_{j=1}\int_{\Omega^{k}} \hat{q}^{k}_{j} \psi^{k}_{j}\psi^{k}_{i} \: dx + \sum^{N+1}_{j=1}\int_{\Omega^{k}} \hat{u}^{k}_{j} \psi^{k}_{j}\frac{d}{dx}(\sqrt{\varepsilon}\psi^{k}_{i}) \: dx = \left[\psi^{k}_{i}(\sqrt{\varepsilon}u^{k}_{h})^{*}\right]^{x^{k}_{r}}_{x^{k}_{l}}, \ 1 \leq i \leq N+1.
\end{align}

To recover the matrix form of the above system, we use the affine mapping to perform integrations on the reference element $\hat{\Omega} = [-1,1]$ and rewrite the DG-FEM semi-discrete formulation as follows

\begin{align}
\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\boldsymbol{\hat{u}}^{k}_{h} - \mathcal{S}^{T}\boldsymbol{\hat{f}}^{k}_{h} + (\mathcal{S}^{\sqrt{\varepsilon}})^{T}\boldsymbol{\hat{q}}^{k}_{h} =& -\boldsymbol{\hat{\psi}}(\xi_{r})f^{*}(x^{k}_{r},x^{k+1}_{l})+\boldsymbol{\hat{\psi}}(\xi_{l})f^{*}(x^{k}_{l},x^{k-1}_{r}) \\
&+ \boldsymbol{\hat{\psi}}(\xi_{r})(\sqrt{\varepsilon}q^{k}_{h})^{*}(x^{k}_{r},x^{k+1}_{l})-\boldsymbol{\hat{\psi}}(\xi_{l})(\sqrt{\varepsilon}q^{k}_{h})^{*}(x^{k}_{l},x^{k-1}_{r}), \\
\frac{h^{k}}{2} \mathcal{M}\boldsymbol{\hat{q}}^{k}_{h} + (\mathcal{\tilde{S}}^{\sqrt{\varepsilon}})^{T}\boldsymbol{\hat{u}}^{k}_{h} =& \ \boldsymbol{\hat{\psi}}(\xi_{r})(\sqrt{\varepsilon}u^{k}_{h})^{*}(x^{k}_{r},x^{k+1}_{l})-\boldsymbol{\hat{\psi}}(\xi_{l})(\sqrt{\varepsilon}u^{k}_{h})^{*}(x^{k}_{l},x^{k-1}_{r}). 
\end{align}

where $\mathcal{M}$, $\mathcal{S}$ and $\boldsymbol{\hat{\psi}}$ have already been defined. We observe however that two new operators have to be introduced

\begin{align}
\mathcal{S}^{\sqrt{\varepsilon}}_{ij} &= \int_{\Omega^{k}} \sqrt{\varepsilon} \psi^{k}_{j}\frac{d\psi^{k}_{i}}{dx} \: dx, \\
\mathcal{\tilde{S}}^{\sqrt{\varepsilon}}_{ij} &= \int_{\Omega^{k}} \psi^{k}_{j}\frac{d}{dx}(\sqrt{\varepsilon}\psi^{k}_{i}) \: dx \ .
\end{align}

As before, we use the monotone Lax-Friedrichs flux for the convective numerical flux $f^{*}$. The developments follow exactly the same steps described for the convective part of the conservation law. For the numerical fluxes $(\sqrt{\varepsilon}q)^{*}$ and $(\sqrt{\varepsilon}u)^{*}$, corresponding to the second order dissipative operator, we choose the Local Discontinuous Galerkin (LDG) flux which yields optimal convergence rates. An usual way to define LDG scheme is as follows

$$
(\sqrt{\varepsilon}u_{h})^{*} = \{\{\sqrt{\varepsilon}\:\}\}u^{+}_{h}, \ (\sqrt{\varepsilon}q_{h})^{*} = \sqrt{\varepsilon^{-}}q^{-}_{h}.
$$

The LDG scheme for derivative operators of second order is a smart way to use the upwinding concept for problems that have no preferred direction of propagation (e.g. dissipation and diffusion problems). If fact, the upwinding in this case is done in a very careful way, ensuring that the upwinding for $u_{h}$ and $q_{h}$ is always done in opposite directions.

Now, using the LDG scheme, we can explicitly define the numerical fluxes $(\sqrt{\varepsilon}u_{h})^{*}$ and $(\sqrt{\varepsilon}q_{h})^{*}$ as follows

\begin{align}
&\boldsymbol{\hat{\psi}}(\xi_{r})(\sqrt{\varepsilon}u^{k}_{h})^{*}-\boldsymbol{\hat{\psi}}(\xi_{l})(\sqrt{\varepsilon}u^{k}_{h})^{*} = \boldsymbol{\hat{\psi}}(\xi_{r})\{\{\sqrt{\varepsilon}\:\}\}u^{k+1}_{h}(x^{k+1}_{l}) - \boldsymbol{\hat{\psi}}(\xi_{l})\{\{\sqrt{\varepsilon}\:\}\}u^{k}_{h}(x^{k}_{l}), \\
&\boldsymbol{\hat{\psi}}(\xi_{r})(\sqrt{\varepsilon}q^{k}_{h})^{*}-\boldsymbol{\hat{\psi}}(\xi_{l})(\sqrt{\varepsilon}q^{k}_{h})^{*} = \boldsymbol{\hat{\psi}}(\xi_{r})\sqrt{\varepsilon(x^{k}_{r})}q^{k}_{h}(x^{k}_{r}) - \boldsymbol{\hat{\psi}}(\xi_{l})\sqrt{\varepsilon(x^{k-1}_{r})}q^{k-1}_{h}(x^{k-1}_{r}).
\end{align}

In order to write the final matrix form of the equations above, we proceed expanding the approximation $u_{h}$ and the average operator, first and after we expand the approximation $q_{h}$.

\begin{align}
\boldsymbol{\hat{\psi}}(\xi_{r})(\sqrt{\varepsilon}u^{k}_{h})^{*}-\boldsymbol{\hat{\psi}}(\xi_{l})(\sqrt{\varepsilon}u^{k}_{h})^{*} &= \boldsymbol{\hat{\psi}}(1) \: \frac{1}{2}\left(\sqrt{\varepsilon(x^{k+1}_{l})}+\sqrt{\varepsilon(x^{k}_{r})}\right)(\boldsymbol{\hat{u}}^{k+1}_{h})^{T}\boldsymbol{\hat{\psi}}(-1) \\
&- \boldsymbol{\hat{\psi}}(-1) \: \frac{1}{2}\left(\sqrt{\varepsilon(x^{k-1}_{r})}+\sqrt{\varepsilon(x^{k}_{l})}\right)(\boldsymbol{\hat{u}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(-1), \\
&= \frac{1}{2}\left(\sqrt{\varepsilon(x^{k+1}_{l})}+\sqrt{\varepsilon(x^{k}_{r})}\right)[\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(-1)]\boldsymbol{\hat{u}}^{k+1}_{h} \\
&- \frac{1}{2}\left(\sqrt{\varepsilon(x^{k-1}_{r})}+\sqrt{\varepsilon(x^{k}_{l})}\right)[\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(-1)]\boldsymbol{\hat{u}}^{k}_{h}, \\
&= \frac{1}{2}\left(\sqrt{\varepsilon(x^{k+1}_{l})}+\sqrt{\varepsilon(x^{k}_{r})}\right)\mathcal{F}^{k+1}_{r}\boldsymbol{\hat{u}}^{k+1}_{h} \\
&- \frac{1}{2}\left(\sqrt{\varepsilon(x^{k-1}_{r})}+\sqrt{\varepsilon(x^{k}_{l})}\right)\mathcal{F}^{k}_{l}\boldsymbol{\hat{u}}^{k}_{h}.
\end{align}

\begin{align}
\boldsymbol{\hat{\psi}}(\xi_{r})(\sqrt{\varepsilon}q^{k}_{h})^{*}-\boldsymbol{\hat{\psi}}(\xi_{l})(\sqrt{\varepsilon}q^{k}_{h})^{*} &= \boldsymbol{\hat{\psi}}(1) \: \sqrt{\varepsilon(x^{k}_{r})} \: (\boldsymbol{\hat{q}}^{k}_{h})^{T}\boldsymbol{\hat{\psi}}(1) \\
&- \boldsymbol{\hat{\psi}}(-1) \: \sqrt{\varepsilon(x^{k-1}_{r})} \: (\boldsymbol{\hat{q}}^{k-1}_{h})^{T}\boldsymbol{\hat{\psi}}(1),  \\
&= \sqrt{\varepsilon(x^{k}_{r})} \: [\boldsymbol{\hat{\psi}}(1)\otimes\boldsymbol{\hat{\psi}}(1)] \: \boldsymbol{\hat{q}}^{k}_{h}  \\
&- \sqrt{\varepsilon(x^{k-1}_{r})} \: [\boldsymbol{\hat{\psi}}(-1)\otimes\boldsymbol{\hat{\psi}}(1)] \: \boldsymbol{\hat{q}}^{k-1}_{h}, \\
&= \sqrt{\varepsilon(x^{k}_{r})} \: \mathcal{F}^{k}_{r}\boldsymbol{\hat{q}}^{k}_{h} \\
&- \sqrt{\varepsilon(x^{k-1}_{r})} \: \mathcal{F}^{k-1}_{l}\boldsymbol{\hat{q}}^{k-1}_{h}.
\end{align}

Now, we are ready to write the final matrix form of the semi-discrete DG-FEM approximation of the modified Burgers' equation

\begin{align}
\frac{h^{k}}{2}\mathcal{M}\frac{d}{dt}\boldsymbol{\hat{u}}^{k}_{h} - \mathcal{S}^{T}\boldsymbol{\hat{f}}^{k}_{h} + (\mathcal{S}^{\sqrt{\varepsilon}})^{T}\boldsymbol{\hat{q}}^{k}_{h} = \: 
 &\left[- \: \frac{1}{2}\mathcal{F}^{k+1}_{r}\left(\boldsymbol{\hat{f}}^{k+1}_{h}-C\boldsymbol{\hat{u}}^{k+1}_{h}\right) - \frac{1}{2}\mathcal{F}^{k}_{r}\left(\boldsymbol{\hat{f}}^{k}_{h}+C\boldsymbol{\hat{u}}^{k}_{h}\right) \right. \\
& \ + \left. \frac{1}{2}\mathcal{F}^{k-1}_{l}\left(\boldsymbol{\hat{f}}^{k-1}_{h}+C\boldsymbol{\hat{u}}^{k-1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{l}\left(\boldsymbol{\hat{f}}^{k}_{h}-C\boldsymbol{\hat{u}}^{k}_{h}\right)\right. \\ 
& \ + \left. \sqrt{\varepsilon(x^{k}_{r})} \: \mathcal{F}^{k}_{r}\boldsymbol{\hat{q}}^{k}_{h} - \sqrt{\varepsilon(x^{k-1}_{r})} \: \mathcal{F}^{k-1}_{l}\boldsymbol{\hat{q}}^{k-1}_{h} \right], \\
\frac{h^{k}}{2} \mathcal{M}\boldsymbol{\hat{q}}^{k}_{h} + (\mathcal{\tilde{S}}^{\sqrt{\varepsilon}})^{T}\boldsymbol{\hat{u}}^{k}_{h} = & \ + \frac{1}{2}\left(\sqrt{\varepsilon(x^{k+1}_{l})}+\sqrt{\varepsilon(x^{k}_{r})}\right)\mathcal{F}^{k+1}_{r}\boldsymbol{\hat{u}}^{k+1}_{h} \\
& \ - \frac{1}{2}\left(\sqrt{\varepsilon(x^{k-1}_{r})}+\sqrt{\varepsilon(x^{k}_{l})}\right)\mathcal{F}^{k}_{l}\boldsymbol{\hat{u}}^{k}_{h}.
\end{align}

### Sub-cell Shock Capturing

The development of the previous alternative formulation with the second order dissipative term is the necessary step to apply the sub-cell shock capturing strategy. To avoid reducing the order of accuracy of the global solution near shocks, the current approach introduces dissipation locally and the amount of viscosity required for stability is determined by the resolution of the approximating space. A shock detection algorithm based on the rate of decay of the expansion coefficients of the approximated solution is used to track the troubled cells where artificial viscosity is need to ensure stability.

#### Shock Indicator

The method used in this work to detect oscillations near shocks is due to Persson \& Peraire. This method is based on the representation of the high-order approximated solution within each element by hierarchical orthogonal polynomials, namely Legendre polynomials. The expansion coefficients of the Legendre basis are treated as a sequence of Fourier modes. For smooth solutions the expansion coefficients are expected to decay very rapidly ($\sim 1/N^{2}$). On the other hand, in the presence of an actual discontinuity all frequency modes are present and the Fourier-like coefficients decrease much more slowly. In this context, solutions of polynomial order $N$ within each element can be represented as

$$
u(x) = \sum^{N+1}_{j=1}\hat{u}^{k}_{j}\psi^{k}_{j}(x),
$$

where $\psi^{k}_{j}(x)$ are the Legendre polynomials and $\hat{u}^{k}_{j}$ are the modal coefficients. For a truncated representation of the same solution at order $N-1$, we write

$$
\tilde{u}(x) = \sum^{N}_{j=1}\hat{u}^{k}_{j}\psi^{k}_{j}(x).
$$

With these two representations, the shock indicator is defined as

$$
F_{k} = \frac{\left\langle u-\tilde{u},u-\tilde{u} \right\rangle_{k}}{\left\langle u,u \right\rangle_{k}},
$$

where $\left\langle \cdot,\cdot \right\rangle_{k}$ represents the standard $L^{2}$ inner product over the element $k$.
The indicator can readily be computed if we recover the following expressions

\begin{align*}
u-\tilde{u} &= \hat{u}^{k}_{N+1}\psi^{k}_{N+1}, \\
\left\langle u-\tilde{u},u-\tilde{u} \right\rangle_{k} &= \int_{\Omega_{k}} (\hat{u}^{k}_{N+1}\psi^{k}_{N+1})^{2} \: dx, \\
\left\langle u,u \right\rangle_{k} &= \int_{\Omega_{k}} (\sum^{N+1}_{j=1}\hat{u}^{k}_{j}\psi^{k}_{j})^{2} \: dx = (\boldsymbol{\hat{u}}^{k}_{h})^{T} \mathcal{M}^{k} \boldsymbol{\hat{u}}^{k}_{h}. 
\end{align*}

#### Artificial Viscosity Quantification

Once the shock has been detected by the previous indicator, the amount of artificial viscosity is taken to be constant over each element and determined by the following function,

$$
	\varepsilon_{k} = \left\{ 
		\begin{array}{ll}
			0 & \text{if} \ s_{k} < s_{0}-\kappa \\
			\frac{\varepsilon_{0}}{2}\left(1+\sin\frac{\pi(s_{k}-s_{0})}{2\kappa}\right) & \text{if} \ s_{0}-\kappa \leq s_{k} \leq s_{0}+\kappa \\
			\varepsilon_{0} & \text{if} \ s_{k} > s_{0}+\kappa,
		\end{array}
	\right.
$$

where $\varepsilon_{0} = \frac{1}{2}max(u_{h})\:h^{k}/N$, $s_{0} = -4.0\log_{10}N$, $s_{k} = \log_{10}S_{k}$, $S_{k} = \min(C \cdot F_{k},1.0)$ and $C = 10^{3}(N^{4})$. The variable $\kappa$ is empirically chosen to obtain a sharp but smooth shock profile. In this work $\kappa = 1.0$.

### Time Marching

The DG-FEM spatial discretization of the Burgers' equation produces a system of first order ODEs.  

Before integrating the semi-discrete problem in time, we define the right hand side operator $\mathcal{L}_{h}$ for the variable $\boldsymbol{u}_{h}$ as

\begin{align}
\mathcal{L}_{h}(\boldsymbol{u}_{h},t) = \frac{2}{h^{k}} \mathcal{M}^{-1} \left\{\mathcal{S}^{T}\boldsymbol{\hat{f}}^{k}_{h} + \left[-\frac{1}{2}\right.\right.&\mathcal{F}^{k+1}_{r}\left(\boldsymbol{\hat{f}}^{k+1}_{h}-C\boldsymbol{\hat{u}}^{k+1}_{h}\right) - \frac{1}{2}\mathcal{F}^{k}_{r}\left(\boldsymbol{\hat{f}}^{k}_{h}+C\boldsymbol{\hat{u}}^{k}_{h}\right) \\
+\frac{1}{2} \: &\mathcal{F}^{k-1}_{l}\left(\boldsymbol{\hat{f}}^{k-1}_{h}+C\boldsymbol{\hat{u}}^{k-1}_{h}\right) + \left.\left.\frac{1}{2}\mathcal{F}^{k}_{l}\left(\boldsymbol{\hat{f}}^{k}_{h}-C\boldsymbol{\hat{u}}^{k}_{h}\right)\right]\right\}
\end{align}

or in the following alternative form when considering the add of the artificial viscosity operator

\begin{align}
\mathcal{L}_{h}(\boldsymbol{u}_{h},t) = \frac{2}{h^{k}} \mathcal{M}^{-1} \left\{ \mathcal{S}^{T}\boldsymbol{\hat{f}}^{k}_{h} - (\mathcal{S}^{\sqrt{\varepsilon}})^{T}\boldsymbol{\hat{q}}^{k}_{h} + 
  \left[- \: \frac{1}{2}\mathcal{F}^{k+1}_{r}\left(\boldsymbol{\hat{f}}^{k+1}_{h}\right.\right.\right. \! -& \left. C\boldsymbol{\hat{u}}^{k+1}_{h}\right) - \frac{1}{2}\mathcal{F}^{k}_{r}\left(\boldsymbol{\hat{f}}^{k}_{h}+C\boldsymbol{\hat{u}}^{k}_{h}\right) \\
+ \ \frac{1}{2}\mathcal{F}^{k-1}_{l}\left(\boldsymbol{\hat{f}}^{k-1}_{h} \right. +& \left. C\boldsymbol{\hat{u}}^{k-1}_{h}\right) + \frac{1}{2}\mathcal{F}^{k}_{l}\left(\boldsymbol{\hat{f}}^{k}_{h}-C\boldsymbol{\hat{u}}^{k}_{h}\right)  \\ 
+ \ \sqrt{\varepsilon(x^{k}_{r})} \: \mathcal{F}^{k}_{r}\boldsymbol{\hat{q}}^{k}_{h} \: -& \left.\left.\sqrt{\varepsilon(x^{k-1}_{r})} \: \mathcal{F}^{k-1}_{l}\boldsymbol{\hat{q}}^{k-1}_{h} \right]\right\}, \\
\boldsymbol{\hat{q}}^{k}_{h} = \frac{2}{h^{k}} \mathcal{M}^{-1} \left\{- (\mathcal{\tilde{S}}^{\sqrt{\varepsilon}})^{T}\boldsymbol{\hat{u}}^{k}_{h} + \frac{1}{2}\left(\sqrt{\varepsilon(x^{k+1}_{l})} \right.\right. +& \left.\sqrt{\varepsilon(x^{k}_{r})}\right)\mathcal{F}^{k+1}_{r}\boldsymbol{\hat{u}}^{k+1}_{h} \\
- \: \frac{1}{2}\left(\sqrt{\varepsilon(x^{k-1}_{r})} \right. +& \left.\left. \sqrt{\varepsilon(x^{k}_{l})}\right)\mathcal{F}^{k}_{l}\boldsymbol{\hat{u}}^{k}_{h}\right\}.
\end{align}

To march this system of equations in time, we can use the classical and reliable forth order four stage explicit Runge-Kutta method (ERK)

\begin{align}
&\boldsymbol{k}^{(1)} = \mathcal{L}_{h}(\boldsymbol{u}^{n}_{h},t^{n}),\\
&\boldsymbol{k}^{(2)} = \mathcal{L}_{h}\left(\boldsymbol{u}^{n}_{h}+\frac{1}{2}\mathit{\Delta} t\boldsymbol{k}^{(1)},t^{n}+\frac{1}{2}\mathit{\Delta} t\right), \\
&\boldsymbol{k}^{(3)} = \mathcal{L}_{h}\left(\boldsymbol{u}^{n}_{h}+\frac{1}{2}\mathit{\Delta} t\boldsymbol{k}^{(2)},t^{n}+\frac{1}{2}\mathit{\Delta} t\right), \\
&\boldsymbol{k}^{(4)} = \mathcal{L}_{h}\left(\boldsymbol{u}^{n}_{h}+\mathit{\Delta} t\boldsymbol{k}^{(3)},t^{n}+\mathit{\Delta} t\right), \\
&\boldsymbol{u}^{n+1}_{h} = \boldsymbol{u}^{n}_{h}+\frac{1}{6}\mathit{\Delta} t\left(\boldsymbol{k}^{(1)}+2\boldsymbol{k}^{(2)}+2\boldsymbol{k}^{(3)}+\boldsymbol{k}^{(4)}\right),
\end{align}

to advance from $\boldsymbol{u}^{n}_{h}$ to $\boldsymbol{u}^{n+1}_{h}$ which are separated by the timestep, $\mathit{\Delta} t$. It is important to note that the ERK method is conditionally stable and requires the timestep size $\mathit{\Delta} t$ be estimated using the $CFL$ (Courant-Friedrichs-Lewy) condition. We define the CFL as

$$
CFL := \max_{k}\frac{c^{k}\mathit{\Delta} t}{\mathit{\Delta} x^{k}},
$$

where $c^{k}$ is the local propagation wave speed on element $k$ and $\mathit{\Delta} x^{k} = h^{k}$ is the local spatial grid resolution. The specific value of the $CFL$ (constant) number depends on the method used for the space and time discretizations but, in general, its order of magnitude is about 1. For high order methods such as the DG scheme described in this text, we have to consider the minimum grid space as the minimum distance between the mapped quadrature points (GL or GLL) on the mesh elements.
	
When considering the traffic flow equation with the artificial viscosity operator, the scaling of the explicit time step is tipically given by

$$
\Delta t \sim \frac{1}{\left|V_{max}\right|\frac{N^{2}}{h}+\left\|\varepsilon\right\|_{L^{\infty}}\frac{N^{4}}{h^{2}}}
$$

where $V_{max}$ is the largest characteristic car speed, $\varepsilon$ is the magnitude of the viscosity, $h$ is the minimum local mesh size and $N$ is the polynomial order of approximation.

## Python implementation

### The finite element mesh

Given that we are on a 1D setting, a very simple mesh generator is easy to implement. Of course, on 2D or 3D we would resort to a third-party application.

The function is _MeshGen1D_, returning:

- `Nv`, number of vertices on the mesh
- `VX`, $x$ coordinate of every node on the mesh, and
- `EToV`, connectivity matrix, associating elements and nodes


```python
# Simple 1D mesh generator
def MeshGen1D():
    # function [Nv, VX, K, EToV] = MeshGen1D(xmin,xmax,K)
    # Purpose  : Generate simple equidistant grid with K elements

    Nv = K+1
    # Generate node coordinates
    VX = numpy.linspace(1,Nv,Nv)

    for i in range(1,Nv+1):
        VX[i-1] = (xmax-xmin)*(i-1)/(Nv-1) + xmin

    # read element to node connectivity
    EToV = numpy.zeros((K,2)) 
    for k in range(0,K): # for k = 1:K
        EToV[k,0] = k+1 #  EToV(k,1) = k; EToV(k,2) = k+1;
        EToV[k,1] = k+2

    return Nv, VX, EToV
```

### Gaussian quadrature points and weights

Having defined the mesh, now we proceed with the creation of the Gaussian quadrature points and integration weights. For this purpose we can resort to the _JacobiGaussQuad_ function, returning:

- `xi`, the integration points, and
- `wi`, the integration weights

Of course we need to fire up all Jacobi polynomial machinery to support the creation of points and weights.

### Jacobi polynomials


```python
#Jacobi polynomials
def JacobiP(r,m,alpha,beta):
# Pm, Pn -- mth-order Jacobi polynomial with weights alpha and beta evaluated at point r
# Pn0, Pn1, Pn2 -- Jacobi polynomials of order n-1, n and n+1
# a1n, a2n, a3n, a4n -- recursive coefficients to compute Jacobi polynomials
    
    #Jacobi polynomials evaluated at point r
    Pn0 = 1.
    Pn1 = 0.5*(alpha - beta + (alpha + beta + 2) * r);
    
    for n in range(1,m):
        a1n = 2.*(n+1)*(n + alpha + beta +1)*(2.*n + alpha + beta)
        a2n = (2.*n + alpha + beta + 1)*(alpha**2 - beta**2)
        a3n = (2.*n + alpha + beta)*(2.*n + alpha + beta + 1)*(2.*n + alpha + beta + 2)
        a4n = 2.*(n + alpha)*(n + beta)*(2.*n + alpha + beta + 2)
        
        Pn2 = (1./a1n)*((a2n + a3n * r)*Pn1 - a4n*Pn0)
        
        Pn0 = Pn1
        Pn1 = Pn2
    
    if m == 1:
        Pm = Pn1
    elif m > 1:
        Pm = Pn2
    else: 
        Pm = Pn0
    
    return Pm
```

### Jacobi polynomials derivatives


```python
#First derivative of Jacobi polynomials
def DJacobiP(r,m,alpha,beta):
# DPm, DPn1 -- Derivative of the mth-order Jacobi polynomial at point r
# Pn0, Pn1, Pn2 -- Jacobi polynomials of order n-1, n and n+1
# a1n, a2n, a3n, a4n -- recursive coefficients to compute Jacobi polynomials
# b1n, b2n, b3n -- recursive coefficients to compute first derivatives of Jacobi polynomials

    # Derivative of the Jacobi polynomial of order m evaluated at point r
    DPm = 0.
    Pn0 = 1.
    Pn1 = 0.5*(alpha - beta + (alpha + beta + 2) * r)
    Pn2 = 0.
    
    for n in range(1,m+1):
        a1n = 2.*(n+1)*(n + alpha + beta +1)*(2.*n + alpha + beta)
        a2n = (2.*n + alpha + beta + 1)*(alpha**2 - beta**2)
        a3n = (2.*n + alpha + beta)*(2.*n + alpha + beta + 1)*(2.*n + alpha + beta + 2)
        a4n = 2.*(n + alpha)*(n + beta)*(2.*n + alpha + beta + 2)

        Pn2 = (1./a1n)*((a2n + a3n * r)*Pn1 - a4n*Pn0)
        
        if r == -1.:
            DPn1 = (-1.0)**(m-1)*0.5*(alpha+beta+m+1)*math.factorial(m+beta)/(math.factorial(beta+1.)*math.factorial(m-1))
        elif r == 1.:
            DPn1 = 0.5*(alpha+beta+m+1.)*math.factorial(m+alpha)/(math.factorial(alpha+1)*math.factorial(m-1))
        else:
            b1n = (2.*n + alpha + beta)*(1 - r**2)
            b2n = (n)*(alpha - beta - (2.*n + alpha + beta) * r)
            b3n = 2.*(n + alpha)*(n + beta)

            DPn1 = (1./b1n)*(b2n*Pn1 + b3n*Pn0)

        Pn0 = Pn1
        Pn1 = Pn2

        DPm = DPn1
        
    return DPm
```

### Roots of Jacobi polynomials through Newton's method with deflation


```python
#Newton-Raphson algorithm with polynomial deflation to compute polynomial roots
def JacobiRoots(alpha,beta,m):
    
    #Allocates vector for quadrature points
    x = numpy.zeros(m)
    #Stop criterium for Newton-Raphson algorithm
    epsilon = 1.e-16

    #Newton-Raphson algorithm with polynomial deflation
    for k in range(0,m):
        #Initial guess: roots of the mth-order Chebyshev polynomial
        r = -numpy.cos(numpy.pi * (2.*k + 1)/(2.*m))
        if k > 0:
            r = (r + x[k-1])/2.
            
        delta = 1.
        while abs(delta) > epsilon:
            s = 0.
            for i in range(0,k):
                s = s + 1./(r - x[i])
                
            Pm = JacobiP(r,m,alpha,beta)
            DPm = DJacobiP(r,m,alpha,beta)
            delta = -Pm/(DPm - Pm * s)
            r = r + delta
            
        x[k] = r
        
    return x
```

### Gaussian quadrature weights


```python
#Quadrature weights
def JacobiGaussWeights(x,alpha,beta,m):
    
    #Allocates vector for quadrature weights
    w = numpy.zeros(m)
    
    if quad_pts == 'GL':
        #case'GL'
        c = 2.**(alpha + beta + 1)

        g1 = alpha + m + 1
        g2 = beta + m + 1
        g3 = alpha + beta + m + 1

        #Computes factorials g1, g2, g3 e M (WARNING: For n integer Gamma[n] = (n-1)!)
        gama1 = math.factorial(g1-1)
        gama2 = math.factorial(g2-1)
        gama3 = math.factorial(g3-1)
        M = math.factorial(m)

        #Computes quadrature weights
        C1 = c*gama1*gama2
        for i in range (0,m):
            r = x[i]
            DPm = DJacobiP(r,m,alpha,beta)
            C2 = M*gama3*(1-r**2)
            w[i] = (C1/C2)*(DPm**(-2))

    elif quad_pts == 'GLL':
        #case 'GLL'
        c = 2.**(alpha + beta + 1)

        g1 = alpha + m
        g2 = beta + m
        g3 = alpha + beta + m + 1

        #Computes factorials g1, g2, g3 e M (WARNING: For n integer Gamma[n] = (n-1)!)
        gama1 = math.factorial(g1-1)
        gama2 = math.factorial(g2-1)
        gama3 = math.factorial(g3-1)
        M = math.factorial(m-1)

        #Computes quadrature weights
        C1 = c*gama1*gama2
        C2 = (m - 1.)*M*gama3
        C3 = beta + 1.
        C4 = alpha + 1.

        r = x[0]
        Pm = JacobiP(r,m-1,alpha,beta)
        w[0] = C3*(C1/C2)*(Pm**(-2))

        r = x[m-1]
        Pm = JacobiP(r,m-1,alpha,beta)
        w[m-1] = C4*(C1/C2)*(Pm**(-2))

        for i in range (1,(m-1)):
            r = x[i]
            Pm = JacobiP(r,m-1,alpha,beta)
            w[i] = (C1/C2)*(Pm**(-2))
            
    else:
        print('Invalid integration type!')
        
    return w
```

### Putting all together


```python
#Quadrature points and weights
import math

def JacobiGaussQuad(N):
    
    #Gets number of quadrature points based on the order of the polynomial approximation
    if quad_pts == 'GL':
        n_GL = (N+1.)/2.
        n = math.ceil(n_GL)
        
    elif quad_pts == 'GLL': 
        n_GLL = (N+3.)/2.
        n = math.ceil(n_GLL) #math.floor(n_GLL) keeps the nodal spectral approach n = N+1
        
    else:
        print('Invalid integration type!')
        
    n=int(n) #JacobiRoots and JacobiGaussWeights dont accept n as float!
    
    #Sets parameters to obtain quadrature points from Legendre polynomials
    #Legendre == Jacobi(alpha,beta) with alpha == 0 and beta == 0
    alpha = 0
    beta = 0
    
    #Sobroutines that evaluates the quadrature points and weights using modified 
    #Newton-Raphson algorithm with polynomial deflation
    if quad_pts == 'GL':
        xi = JacobiRoots(alpha,beta,n)
        wi = JacobiGaussWeights(xi,alpha,beta,n)
        
    elif quad_pts == 'GLL':
        #Allocating vector for quadrature points
        xi = numpy.zeros(n)
        x = JacobiRoots(alpha+1,beta+1,n-2)
        xi[1:-1] = x.copy()
        xi[0] = -1.0
        xi[n-1] = 1.0
        wi = JacobiGaussWeights(xi,alpha,beta,n)
        
    else:
        print('Invalid integration type!')
    
    return xi, wi
```

### Basis functions

For this specific problem we will resort to the Legendre orthogonal/orthonormal basis, implemented respectively in _LegendreBasis_ and _OrthonormalLegendreBasis_ functions. We also need the basis derivative, so we present also the collocation matrix used for derivation.

### Orthogonal Legendre basis


```python
#Orthogonal polynomial basis
def LegendreBasis(N,xi):
    
    # Number of basis functions
    Q = N+1
    # Number of evaluation points
    npts = len(xi) #npts = size(xi,1)
    # Allocates matrix to store Legendre basis evaluated at npts points
    Lj = numpy.zeros((npts,Q))
    
    #Legendre basis
    for j in range(1,Q+1):
        for i in range(1,npts+1):
            Lk = JacobiP(xi[i-1],j-1,0,0) 
            Lj[i-1,j-1] = Lk
            #gamma = 2/(2*(j-1)+1); orthonormalization process
            #Lj(i,j) = sqrt(gamma)*Lk;
            
    return Lj
```

### Orthonormal Legendre basis


```python
#Orthonormal polynomial basis
def OrthonormalLegendreBasis(N,xi):
    
    # Number of basis functions
    Q = N+1
    # Number of evaluation points
    npts = len(xi) #npts = size(xi,1)
    # Allocates matrix to store Legendre basis evaluated at npts points
    Lj = numpy.zeros((npts,Q))
    
    #Legendre basis
    for j in range(0,Q):
        for i in range(0,npts):
            Lk = JacobiP(xi[i],j,0,0) 
            gamma = 2.0/(2.0*j+1.0) # orthonormalization process
            Lj[i,j] = Lk/(numpy.sqrt(gamma))
            
    return Lj
```

### Auxiliar derivatives to assemble the collocation matrix


```python
#First derivative of the auxiliar polynomial PQ
def dPQ(ri,Q,alpha,beta):
    
    #case 'GL'
    if quad_pts == 'GL':
        dPq = DJacobiP(ri,Q,alpha,beta)
        
    #case 'GLL'
    elif quad_pts == 'GLL':
        
        if ri == -1.0:
            c = 2.*(-1)**Q
            gama1 = math.factorial(Q + beta - 1)
            gama2 = math.factorial(Q - 2)
            gama3 = math.factorial(beta + 1)

            dPq = (c*gama1)/(gama2*gama3)
        elif ri == 1.0:
            c = -2.
            gama1 = math.factorial(Q + alpha - 1)
            gama2 = math.factorial(Q - 2)
            gama3 = math.factorial(alpha + 1)

            dPq = (c*gama1)/(gama2*gama3)
        else:
            PJacobi = JacobiP(ri,Q-1,alpha,beta)

            dPq = -2.*(Q - 1)*PJacobi
            
    else:
        print('Invalid integration type!')       
        
    return dPq
```


```python
#Second derivative of the auxiliar polynomial PQ
def d2PQ(ri,Q,alpha,beta):
    
    #case 'GL'
    if quad_pts == 'GL':
        DPJacobi = DJacobiP(ri,Q,alpha,beta)
        ddPq = ( (alpha - beta + (alpha + beta + 2)*ri)/(1 - ri**2) )*DPJacobi
        
    #case 'GLL'
    elif quad_pts == 'GLL':
    
        if ri == -1.0:
            c = (((-1)**Q)*2*(alpha - (Q - 1)*(Q + alpha + beta)))/(beta + 2)
            gama1 = math.factorial(Q + beta - 1)
            gama2 = math.factorial(Q - 2)
            gama3 = math.factorial(beta + 1)
            
            ddPq = (c*gama1)/(gama2*gama3)
        elif ri == 1.0:
            c = (2*(beta - (Q - 1)*(Q + alpha + beta)))/(alpha + 2)
            gama1 = math.factorial(Q + alpha - 1)
            gama2 = math.factorial(Q - 2)
            gama3 = math.factorial(alpha + 1)
            
            ddPq = (c*gama1)/(gama2*gama3)
        else:
            PJacobi = JacobiP(ri,Q-1,alpha,beta)
            DPJacobi = (-2*(Q - 1))/((1 - ri)*(1 + ri))*PJacobi
            
            ddPq = (alpha - beta + (alpha + beta)*ri)*DPJacobi
                         
    else:
        print('Invalid integration type!') 
            
    return ddPq
```

### The collocation differentiation matrix


```python
#Collocation differeniation matrix
def DMatrix1D(xi,nip):
    
    # Sets parameters to obtain quadrature points from Legendre polynomials
    # Legendre == Jacobi(alpha,beta) with alpha == 0 and beta == 0
    if nip == 1:
        D = 0
        
    else:  
        alpha = 0
        beta = 0
        D = numpy.zeros((nip,nip))

        for i in range(0,nip):
            ri = xi[i]

            dPi = dPQ(ri,nip,alpha,beta)
            d2Pi = d2PQ(ri,nip,alpha,beta)

            for j in range(0,nip):
                rj = xi[j]
                dPj = dPQ(rj,nip,alpha,beta)
                if i == j:
                    D[i,i] = d2Pi/(2.*dPi)
                else:
                    D[i,j] = (dPi/dPj)*(1./(ri - rj))
                
    return D
```

### Jacobian

As we have seen in the theoretical part, there is a Jacobian involved in the numerical integration due to the transformation of all integrals from the real elements to the standard one. The Jacobian is evaluated, unsurprisingly, by the _Jacobian_ function. It returns:

- _J_, the Jacobian matrix, and
- _xcoord_, a mapping of the integration points (remember, they are all in [-1,1]) to the real elements


```python
def Jacobian():
    
    J = numpy.zeros(K)
    B = numpy.zeros(K)
    xcoord = numpy.zeros((nip,K))
    
    J = 0.5*(VX[1:] - VX[:-1])
    B = 0.5*(VX[1:] + VX[:-1])
                                  
    # [xcoord] = [xi]*[J]^T + [ones]*[B]^T assuming [xi], [J], [ones] and [B] as column vectors 
    xcoord = numpy.dot(xi.reshape((nip,1)),J.reshape((1,K))) + numpy.dot(numpy.ones(nip).reshape((nip,1)),B.reshape((1,K)))
    
    return xcoord, J
```

### Calculating $\Delta t$ and the number of time steps


```python
def Timestep():
    
    # smalest grid space
    xmin = min(abs(xc[0,:]-xc[1,:]))
    uic = 0.5 + numpy.sin(2.0*numpy.pi*xc[:,:])
    umax = max(abs(uic.reshape((nip*K))))

    #CFL=0.25 
    #dt = CFL*xmin/Vm
    h = 2*J[0] 
    ni_max = 0.005
    dt = (umax*N**2/h + ni_max*N**4/h**2)**(-1)
    Nsteps = math.ceil(tmax/dt) 
    dt = tmax/Nsteps
    
    #dt = 0.0005
    #dt = 0.00025    
    #Nsteps = 5000 # reduced number of time steps to evaluate overall behaviour
    
    return dt, int(Nsteps), umax
```

### Projecting the initial condition onto the polynomial space

As initial condition we are using

$$
u(x,0) = \frac{1}{2} + \sin(2\pi\,x).
$$

The function _IC_ is responsible for projecting $u(x,0)$ onto the polynomial space, as well as, evaluating the mass matrix. Note that due to the orthogonal nature of the basis, the mass matrix is diagonal.


```python
def InviscidBurgers():
    
    fxi = 0.5 + numpy.sin(2.0*numpy.pi*xc[:,:])
    
    return fxi

def IC():
    
    #Allocates vector for the initial condition 
    uhat = numpy.zeros((Nldof,K))
    # IC
    fxi = InviscidBurgers()
    
    #Allocates vector for the inverse of the diagonal Mass matrix
    invdiag = numpy.zeros(Nldof)
    
    for i in range(0,Nldof):
        #d[i] = sum(wi*psi[:,i]**2)
        invdiag[i] = 1./(sum(wi*psi[:,i]**2))
        
    # Defining the diagonal Mass matrix
    #m = numpy.dot(wi*psi.T,psi)   
    #m = numpy.diag(d)
    
    # RHS of the projection operation
    b = numpy.dot(wi*psi.T,fxi)
    
    # Projecting IC into the polynomial space
    #fproj = LA.solve(m,b)
    #phi = fproj.T.reshape((Nldof*K,1))  

    fproj = invdiag*b.T
    uhat = fproj.T.copy()
    
    return invdiag, uhat
```

### Stiffness matrix $\mathcal{S}$


```python
def StiffMatrix():

    #Allocates local stifness matrix
    S = numpy.zeros((Nldof,Nldof))

    #Computes local stiffness matrix
    S = numpy.dot(wi*psi.T,Dpsi)
    
    return S
```

### Lift matrices $\mathcal{F}_r^{k,k+1}$ and $\mathcal{F}_l^{k-1,k}$


```python
def LiftMatrix():
    
    #Allocates lift matrices
    Frk = numpy.ones((Nldof,Nldof))
    Flk = numpy.ones((Nldof,Nldof))
    Flkp1 = numpy.ones((Nldof,Nldof))
    Frkm1 = numpy.ones((Nldof,Nldof))
    
    #Computes lift matrices
    Flkp1[:,1:Nldof:2] = (-1.)*Flkp1[:,1:Nldof:2]
    Frkm1[1:Nldof:2,:] = (-1.)*Frkm1[1:Nldof:2,:]
    
    Flk[:Nldof:2,1:Nldof:2] = (-1.)*Flk[:Nldof:2,1:Nldof:2]
    Flk[1:Nldof:2,:Nldof:2] = (-1.)*Flk[1:Nldof:2,:Nldof:2]
    
    return Flk, Frk, Flkp1, Frkm1
```

### Projection of the numerical flux onto the polynomial space


```python
def FluxProjection(ft,ut):
        
    uhat = ut[:,1:-1,0].copy()
    
    #Recovers solution at (nonlinear) quadrature points for over integration
    uh = numpy.dot(nlpsi,uhat)
    
    #Reconstructs nonlinear flux at quadrature points for over integrataion
    fh = (uh[:,:]**2)/2.0
    
    # RHS of the projection operation
    b = numpy.dot(nlwi*nlpsi.T,fh)
    
    #Physical flux projection onto the polynomial space
    Fluxproj = InvM*b.T
    
    ft[:,1:-1] = Fluxproj.T.copy()
    
    return ft
```

### Persson and Peraire shock indicator


```python
def ShockIndicator(uti,i,sp):
    
    #Gets element size
    hk = 2*J[i-1]
    
    #Sensor parameters
    #s0 = -(4.0 + 4.25*math.log10(N))
    s0 = -4.*math.log10(N)
    kappa = 1.0
    #ni0 = 0.5*umax*hk/N
    ni0 = 0.5*sp*hk/N
    #ni0 = hk/N
    
    
    #Recovering Mass matrix
    M = 1./InvM
    Mk = (hk/2.)*M
    
    #Smoothness indicator
    InnerProdUdiff = uti[N]**2*(hk/2.)*numpy.dot(wi*psi[:,N],psi[:,N])
    InnerProdU = numpy.dot(uti[:],Mk*uti[:])
    
    Fk = 0.
    #nik = 0.
    
    #Avoids divisoin by zero
    if InnerProdU < 1.e-15:
        nik = 0.
    
    #Computes artificial viscosity amount based on the smoothness indicator factors Fk and sk
    else:
    
        Fk = InnerProdUdiff/InnerProdU

        #C = 1000.*(5.0*N**4+1.0) #Exponential IC example
        #C = 0.0001*(N**4)      #Inviscid Burgers' example
        C = 100.0*(N**4)
        #C = 10000.*(5.0*N**4)    #Red light example
        Sk = numpy.minimum(C*Fk,1.0)

        sk = math.log10(Sk)
        #sk = math.log10(Fk)

        #Sets element artificial viscosity amout
        if sk < (s0-kappa):
            nik = 0.
        elif (s0-kappa) <= sk <= (s0+kappa):
            nik = (ni0/2.0)*(1.0 + numpy.sin(numpy.pi*(sk-s0)/(2.0*kappa)))
        elif sk > (s0+kappa):
            nik = ni0 
    
    return nik, Fk
```

### $\mathcal{L}_h$ operator

Note that there are several commented out options for viscous fluxes and boundary conditions.


```python
def Lh(phi,t,epst):
    
    uhat = phi.copy()
    
    #Recovers physical space solution at quadrature points
    uh = numpy.dot(psi,uhat)
    
    #C = abs(umax)
    C = max(abs(uh.reshape((nip*K))))
    
    #Auxiliary arrays
    ut = numpy.zeros((Nldof,K+2,2))
    qt = numpy.zeros((Nldof,K+2))
    ft = numpy.zeros((Nldof,K+2))
    #epst = numpy.zeros((nip,K+2))
    #Fk = numpy.zeros(K)
    
    #Solution u at time step t
    ut[:,1:-1,0] = phi.copy()
    
    #Shock Indicator and computation of the Elemental Artificial Viscosity
    #for i in range(1,K+1):
    #    (epst[:,i],Fk[i-1]) = ShockIndicator(ut[:,i,0],i,C)
    
    #Apply BCs ----------------------------
    #epst[:,0] = epst[:,1]
    #epst[:,-1] = epst[:,-2]
    
    #Recovers the physical solution uh at the left boundary of the first element and at
    #the right boundary of the last element (using the first row of matrices Flkp1 and Frk)
    uhL = numpy.dot(Flkp1[0,:],ut[:,1,0])
    uhR = numpy.dot(Frk[0,:],ut[:,-2,0])
    
    #Non homogeneous Dirichlet at left boundary
    #ut[0,0,0] = -uhL + 0.5*(1.0 + numpy.sin(t))    #entrance
    #ut[0,-1,0] = -uhR + 0.5*(1.0 + numpy.sin(t))    #exit
    #ut[0,-1,0] = ut[0,0,0]

    #Inviscid Burgers' example
    ut[0,0,0] = -uhL + 2.0*(0.5 + numpy.sin(-2.0*numpy.pi*uhL*t)) #entrance
    ut[0,-1,0] = ut[0,0,0]  #exit
    
    #Homogeneous Neumann at right boundary
    #ut[:,-1,0] = ut[:,-2,0].copy()
    #ut[0,-1,0] = uhR
    
    #Recovers qt from ut to apply BCs
    #Non homogeneous Dirichlet (Left)
    
    #LDG scheme - Persson-Peraire approach
    #qhL = InvM*( J[0]**(-1)*(-epst[0,1]*numpy.dot(S.T,ut[:,1,0]) \
    #                        + 0.5*(epst[0,2] + epst[0,1])*numpy.dot(Flkp1,ut[:,2,0]) \
    #                        - 0.5*(epst[0,0] + epst[0,1])*numpy.dot(Flk,ut[:,1,0]) ) ).T
    
    #LDG scheme - Hesthaven-Warburton approach
    qhL = InvM*( J[0]**(-1)*(-numpy.sqrt(epst[0,1])*numpy.dot(S.T,ut[:,1,0]) \
                       + 0.5*(numpy.sqrt(epst[0,2]) + numpy.sqrt(epst[0,1]))*numpy.dot(Flkp1,ut[:,2,0]) \
                       - 0.5*(numpy.sqrt(epst[0,0]) + numpy.sqrt(epst[0,1]))*numpy.dot(Flk,ut[:,1,0]) ) ).T
    
    #BR1 scheme (central flux) - Hesthaven-Warburton approach
    #qhL = InvM*( J[0]**(-1)*(-numpy.sqrt(epst[0,1])*numpy.dot(S.T,ut[:,1,0]) \
    #                   + 0.5*(numpy.sqrt(epst[0,2])*numpy.dot(Flkp1,ut[:,2,0]) + numpy.sqrt(epst[0,1])*numpy.dot(Frk,ut[:,1,0])) \
    #                   - 0.5*(numpy.sqrt(epst[0,0])*numpy.dot(Frkm1,ut[:,0,0]) + numpy.sqrt(epst[0,1])*numpy.dot(Flk,ut[:,1,0])) ) ).T
    
    qt[0,0] = numpy.dot(Flkp1[0,:],qhL)
    #qt[:,0] = qhL.copy()
    
    ##Homogeneous Neumann (Right)
    #Non-Homogeneous Dirichlet (Right)
    
    #LDG scheme - Persson-Peraire approach
    #qhR = InvM*( J[-1]**(-1)*(-epst[0,-2]*numpy.dot(S.T,ut[:,-2,0]) \
    #                         + 0.5*(epst[0,-1] + epst[0,-2])*numpy.dot(Flkp1,ut[:,-1,0]) \
    #                         - 0.5*(epst[0,-3] + epst[0,-2])*numpy.dot(Flk,ut[:,-2,0]) ) ).T
    
    #LDG scheme - Hesthaven-Warburton approach
    qhR = InvM*( J[-1]**(-1)*(-numpy.sqrt(epst[0,-2])*numpy.dot(S.T,ut[:,-2,0]) \
                         + 0.5*(numpy.sqrt(epst[0,-1]) + numpy.sqrt(epst[0,-2]))*numpy.dot(Flkp1,ut[:,-1,0]) \
                         - 0.5*(numpy.sqrt(epst[0,-3]) + numpy.sqrt(epst[0,-2]))*numpy.dot(Flk,ut[:,-2,0]) ) ).T
    
    #BR1 scheme (central flux) - Hesthaven-Warburton approach
    #qhR = InvM*( J[-1]**(-1)*(-numpy.sqrt(epst[0,-2])*numpy.dot(S.T,ut[:,-2,0]) \
    #                   + 0.5*(numpy.sqrt(epst[0,-1])*numpy.dot(Flkp1,ut[:,-1,0]) + numpy.sqrt(epst[0,-2])*numpy.dot(Frk,ut[:,-2,0])) \
    #                   - 0.5*(numpy.sqrt(epst[0,-3])*numpy.dot(Frkm1,ut[:,-3,0]) + numpy.sqrt(epst[0,-2])*numpy.dot(Flk,ut[:,-2,0])) ) ).T
    
    qt[0,-1] = numpy.dot(Frk[0,:],qhR) # minus sign before numpy.dot() for Neumann conditions
    #qt[:,-1] = -qhR.copy()
    
    #Physical flux BC
    #ft[0,0] = Vm*(ut[0,0,0] - (ut[0,0,0]**2)/um)
    #ft[0,-1] = Vm*(uhR - (uhR**2)/um)
    #ft[0,-1] = Vm*(ut[0,-1,0] - (ut[0,-1,0]**2)/um)
    ft[0,0] = (ut[0,0,0]**2)/2.0
    ft[0,-1] = ft[0,0]
    
    #----------------------------
    
    #Physical non linear flux projection onto the polynomial space using over integration
    ft = FluxProjection(ft,ut)
    
    #Updates solution coefficients form time t to time t+1
    
    #Dissipative operator discretized by LDG scheme and Persson-Peraire approach
    #qt[:,1:-1] = (InvM*( J[:]**(-1)*(- epst[0,1:-1]*numpy.dot(S.T,ut[:,1:-1,0]) \
    #                                 + 0.5*(epst[0,2:] + epst[0,1:-1])*numpy.dot(Flkp1,ut[:,2:,0]) \
    #                                 - 0.5*(epst[0,:-2] + epst[0,1:-1])*numpy.dot(Flk,ut[:,1:-1,0]) ) ).T).T.copy()
    
    #ut[:,1:-1,1] = (InvM*( J[:]**(-1)*( numpy.dot(S.T,ft[:,1:-1]) - epst[0,1:-1]*numpy.dot(S.T,qt[:,1:-1]) \
    #                                  - 0.5*numpy.dot(Flkp1,(ft[:,2:] - C*ut[:,2:,0])) \
    #                                  - 0.5*numpy.dot(Frk,(ft[:,1:-1] + C*ut[:,1:-1,0])) \
    #                                  + 0.5*numpy.dot(Frkm1,(ft[:,:-2] + C*ut[:,:-2,0])) \
    #                                  + 0.5*numpy.dot(Flk,(ft[:,1:-1] - C*ut[:,1:-1,0])) \
    #                                  + numpy.dot(Frk,qt[:,1:-1])  \
    #                                  - numpy.dot(Frkm1,qt[:,:-2]) ) ).T).T.copy()
    
    #Dissipative operator discretized by LDG scheme and Hesthaven-Warburton approach
    qt[:,1:-1] = (InvM*( J[:]**(-1)*(- numpy.sqrt(epst[0,1:-1])*numpy.dot(S.T,ut[:,1:-1,0]) \
                                    + 0.5*(numpy.sqrt(epst[0,2:]) + numpy.sqrt(epst[0,1:-1]))*numpy.dot(Flkp1,ut[:,2:,0]) \
                                    - 0.5*(numpy.sqrt(epst[0,:-2]) + numpy.sqrt(epst[0,1:-1]))*numpy.dot(Flk,ut[:,1:-1,0]) ) ).T).T.copy()
    
    ut[:,1:-1,1] = (InvM*( J[:]**(-1)*( numpy.dot(S.T,ft[:,1:-1]) - numpy.sqrt(epst[0,1:-1])*numpy.dot(S.T,qt[:,1:-1]) \
                                      - 0.5*numpy.dot(Flkp1,(ft[:,2:] - C*ut[:,2:,0])) \
                                      - 0.5*numpy.dot(Frk,(ft[:,1:-1] + C*ut[:,1:-1,0])) \
                                      + 0.5*numpy.dot(Frkm1,(ft[:,:-2] + C*ut[:,:-2,0])) \
                                      + 0.5*numpy.dot(Flk,(ft[:,1:-1] - C*ut[:,1:-1,0])) \
                                      + numpy.sqrt(epst[0,1:-1])*numpy.dot(Frk,qt[:,1:-1])  \
                                      - numpy.sqrt(epst[0,:-2])*numpy.dot(Frkm1,qt[:,:-2]) ) ).T).T.copy()
    
    #Dissipative operator discretized by BR1 scheme (central flux) and Hesthaven-Warburton approach
    #qt[:,1:-1] = (InvM*( J[:]**(-1)*(- numpy.sqrt(epst[0,1:-1])*numpy.dot(S.T,ut[:,1:-1,0]) \
    #                                + 0.5*(numpy.sqrt(epst[0,2:])*numpy.dot(Flkp1,ut[:,2:,0]) + \
    #                                       numpy.sqrt(epst[0,1:-1])*numpy.dot(Frk,ut[:,1:-1,0])) \
    #                                - 0.5*(numpy.sqrt(epst[0,:-2])*numpy.dot(Frkm1,ut[:,:-2,0]) + \
    #                                       numpy.sqrt(epst[0,1:-1])*numpy.dot(Flk,ut[:,1:-1,0])) ) ).T).T.copy()
    
    #ut[:,1:-1,1] = (InvM*( J[:]**(-1)*( numpy.dot(S.T,ft[:,1:-1]) - numpy.sqrt(epst[0,1:-1])*numpy.dot(S.T,qt[:,1:-1]) \
    #                                  - 0.5*numpy.dot(Flkp1,(ft[:,2:] - C*ut[:,2:,0])) \
    #                                  - 0.5*numpy.dot(Frk,(ft[:,1:-1] + C*ut[:,1:-1,0])) \
    #                                  + 0.5*numpy.dot(Frkm1,(ft[:,:-2] + C*ut[:,:-2,0])) \
    #                                  + 0.5*numpy.dot(Flk,(ft[:,1:-1] - C*ut[:,1:-1,0])) \
    #                                  + 0.5*(numpy.sqrt(epst[0,2:])*numpy.dot(Flkp1,qt[:,2:]) + \
    #                                         numpy.sqrt(epst[0,1:-1])*numpy.dot(Frk,qt[:,1:-1]))  \
    #                                  - 0.5*(numpy.sqrt(epst[0,:-2])*numpy.dot(Frkm1,qt[:,:-2]) + \
    #                                         numpy.sqrt(epst[0,1:-1])*numpy.dot(Flk,qt[:,1:-1])) ) ).T).T.copy()
    
    ki = ut[:,1:-1,1].copy()
    
    return ki
```

### Four-step Runge-Kutta time integrator


```python
def RK4(phi,n):
    
    global u, visc, ind, wavespeed
    
    uhat = phi.copy()
    
    #Recovers physical space solution at quadrature points
    uht = numpy.dot(psi,uhat)
    
    #C = abs(umax)
    C = max(abs(uht.reshape((nip*K))))
    
    #Vector containing the artificial viscosity computed for each element
    eps = numpy.zeros((nip,K+2))
    Fk = numpy.zeros(K)
    
    #Shock Indicator and computation of the Elemental Artificial Viscosity
    for i in range(1,K+1):
        (eps[:,i],Fk[i-1]) = ShockIndicator(phi[:,i-1],i,C)
        
    eps[:,0] = eps[:,1]
    eps[:,-1] = eps[:,-2]
    
    #Runge-Kutta 4-steps
    k1 = Lh(phi,t,eps)
    k2 = Lh(phi + dt/2.*k1,t + 0.5*dt,eps)
    k3 = Lh(phi + dt/2.*k2,t + 0.5*dt,eps)
    k4 = Lh(phi + dt*k3,t + dt,eps)
    phi = phi + (dt/6.)*(k1 + 2.*k2 + 2.*k3 + k4)
    
    uh = numpy.dot(psi,phi).T.reshape((nip*K))
    
    u[:,n] = uh.copy()
    
    #Copying artificial viscosity for plotting purposes
    nit = eps[:,1:-1].copy()
    visc[:,n] = nit.T.reshape((nip*K)).copy()
    ind[:,n] = Fk.copy()
    wavespeed[n] = C
    
    return phi
```

### Putting all together: the main program


```python
# Inviscid Burgers' equation: du/dt + d/dx(f(u)) = 0, x in [0,2*pi]

# Nonlinear flux
# f(u) = (u^2)/2;

# Monotone Lax-Friedrich flux: 
# f*(u-,u+) = {{f(uh)}} + C/2*[[uh]]
# C = max|df/du| = max|u|

# IC: u(x,0) = 0.5 + sin(2.0*pi*x)
# BC: u(0,t) = 0.5 + sin(-2.0*pi*uL*t)) - Dirichlet
#     u(L,t) = 0.5 + sin(-2.0*pi*uR*t)) - Dirichlet

# Main

K = 20          #Total number of elements
N = 10          #Polynomial order of approximation
P = 2*N         #Polynomial order of the integrands associated to the Mass matrix coefficients
nlP = 3*N       #Polynomial order of the integrands associated to the non-linear terms
Nldof = N+1     #Number of local degrees of freedom per element

# Artificial dissipation coefficient (test with fixed global artificial viscosity)
#eps = 0.05 # GLL
#eps = 0.025 # GL
#eps = 0.03

# Spatial bounds of the domain [km]
#xb = [0.,2.*numpy.pi] #Exponential IC example
xb = [0.0,1.0] #Green light and Red light examples
xmin = xb[0]
xmax = xb[1]
#Domain size [km]
L = xb[1]-xb[0]

#Temporal bounds [s]
tb = [0.,0.5]
tmin = tb[0]
tmax = tb[1]

#Point distribution for the quadrature rule
#quad_pts = 'GLL' # Gauss-Legendre-Lobatto
quad_pts = 'GL' # Gauss-Legendre

#Quadrature points and weights
(xi,wi) = JacobiGaussQuad(P)
nip = len(xi)

#Quadrature points and weights for the non linear term
(nlxi,nlwi) = JacobiGaussQuad(nlP)
nlnip = len(nlxi)

#Basis functions
(psi) = LegendreBasis(N,xi)
(nlpsi) = LegendreBasis(N,nlxi)
#(V) = OrthonormalLegendreBasis(N,xi)

#Derivative of basis functions
Dhj = DMatrix1D(xi,nip) #WARNING! P=2 at least. (P = 2 => N = 1)
Dpsi = numpy.dot(Dhj,psi)

#1D mesh generation
(Nv, VX, EToV) = MeshGen1D()

#Global coordinates of quadrature points and Jacobian of the local to global mapping
(xc,J) = Jacobian()

# Time step size estimation
(dt,Nsteps,umax) = Timestep()

Nsteps = 10000
# Main vectors that evolve in time: car density (u), artificial viscosity (visc) and auxiliar matrix (phi)
u = numpy.zeros((nip*K,Nsteps+1))
visc = numpy.zeros((nip*K,Nsteps+1))
ind = numpy.zeros((K,Nsteps+1))
wavespeed = numpy.zeros(Nsteps+1)
#wsp = numpy.zeros(4)

#Initial conditions
(InvM,uhat) = IC()
u[:,0] = numpy.dot(psi,uhat).T.reshape((nip*K))

#Plot IC
pyplot.figure(1)
pyplot.plot(xc.T.reshape((nip*K)), u[:,0], color='#003366', ls='-', lw=3)
pyplot.ylim(-0.8,1.8)

#Stiffness matrix
S = StiffMatrix()

#Element boundary matrices (lift operator)
(Flk,Frk,Flkp1,Frkm1) = LiftMatrix()

#Filter parameters
#Nc = 0.
#s = 32.

phi = uhat.copy()

#Runge Kutta time marching
t = tmin
dt = 0.00005
for n in range(1,Nsteps+1):
    phi = RK4(phi,n)
    t = t + dt
```

### Viewing the solution

If you have installed the JSAnimation library you can generate a small movie showing the evolution of the solution along the time. If you get a message __`ImportError: No module named JSAnimation.IPython_display`__ do not dispair! We can still resort to matplotlib to plot a few steps.

#### Creating a movie


```python
from matplotlib import animation
from JSAnimation.IPython_display import display_animation
```


```python
# Lets make some movies!
# From https://jakevdp.github.io/blog/2012/08/18/matplotlib-animation-tutorial/

# First set up the figure, the axis, and the plot element we want to animate
fig = pyplot.figure(2)
ax = pyplot.axes(xlim=(xb[0], xb[1]), ylim=(numpy.min(u), numpy.max(u)))
ax.grid(True)
line, = ax.plot([], [], lw=3)

# Initialization function: plot the background of each frame
def init():
    line.set_data([], [])
    return line,

# Animation function.  This is called sequentially
def animate(data):
    x = xc.T.flat
    y = data
    line.set_data(x,y)
    return line,

# Call the animator.  blit=True means only re-draw the parts that have changed.

# Choose a portion of the solution to speed things up
v = u[:,4000:5000]

anim = animation.FuncAnimation(fig, animate, init_func=init, frames=v.T, interval=50)
display_animation(anim, default_mode='once')
```

#### Plotting some steps of the solution

The solution _u_ has _Nsteps_ columns, precisely one for each time step. In order to view some is simple, as shown in the code below.


```python
# Plotting a few samples of the solution along the time
pyplot.figure(2)
pyplot.clf()
x = xc.T.flat
pyplot.plot(x, u[:,1])    # tstep = 1
pyplot.plot(x, u[:,5020]) # tstep = 5020
pyplot.plot(x, u[:,10000])# tstep = 10000
pyplot.show()
```

## Viscous Burgers'

The viscous Burgers' (VB) equation is composed by its inviscid counterpart (IB) with an extra viscous term

$$
\frac{\partial u(x,t)}{\partial t} +  \frac{\partial f(u(x,t))}{\partial x} = \frac{\partial}{\partial x}\biggr(\varepsilon(x) \frac{\partial u}{\partial x}\biggl), \ x \in [L,R] = \Omega
$$

with initial conditions

$$
u(x,0) = u_{0}(x) 
$$

and boundary conditions

$$
\begin{align}
u(0,t) &= g_0(t) \\
u(L,t) &= g_1(t)
\end{align}
$$
with $u$ and $f$ as before.

If we consider $\varepsilon(x) = \epsilon =$ constant, we obtain

$$
\frac{\partial u(x,t)}{\partial t} +  \frac{\partial f(u(x,t))}{\partial x} = \epsilon \frac{\partial^2 u}{\partial x^2}, \ x \in [L,R] = \Omega
$$

The obtained equation is a _simpler form_ of that presented in the section "_DG-FEM applied to second-order PDEs_", in the sense that we are considering a dissipation coefficient $\varepsilon$ constant in the whole domain, as well as in time.

## Challenge

> Modify the code above in order to solve the viscous Burgers' equation.


```python

```

db5fa4a13bdbbf0dc11f11f13efa05b024793560
ipynb
Jupyter Notebook
lessons/08_dg/08_02_Burgers-DG.v3.ipynb
albertonogueira/numerical-mooc
dd95e650310502b5cdfe6e405ed7ab7e1496d233
CC-BY-3.0
2017-02-10T12:09:09.000Z
2017-02-10T12:09:09.000Z
lessons/08_dg/08_02_Burgers-DG.v3.ipynb
albertonogueira/numerical-mooc
dd95e650310502b5cdfe6e405ed7ab7e1496d233
CC-BY-3.0
lessons/08_dg/08_02_Burgers-DG.v3.ipynb
albertonogueira/numerical-mooc
dd95e650310502b5cdfe6e405ed7ab7e1496d233
CC-BY-3.0
Qwen/Qwen-72B
1. YES
2. YES