# SUPG Stabilization

For streamline upwind/Petrov-Galerkin (SUPG) (see for example
{cite:t}`john:knobloch:2006,clevenger:heister:2019`), we add to the weak form
$a(\cdot,\cdot)$ the cell-wise defined weak form
```{math}
a_{\text{SUPG}} (T, \varphi) =
 \sum_{K \in \mathcal{T}_h}
  \delta_K \left( \rho C_p \frac{\partial T}{\partial t} - k \triangle T + \mathbf{\beta} \cdot \nabla T - F, \mathbf{\beta} \cdot \nabla \varphi \right)_K,
```
where $K \in \mathcal{T}_h$ are the cells in the computation,
$\delta_K \geq 0$ is a stabilization coefficient defined on each cell,
$\mathbf{\beta} = \rho C_p \mathbf{u}$ is the effective advection velocity.
The standard literature about SUPG does not contain $\rho C_p$, so it makes
sense to include this in the velocity. The first argument in the inner product
is the strong form of the residual of PDE, which is tested with the expression
$\mathbf{\beta} \cdot \nabla \varphi$ representing the solution in streamline
direction. We have to assume $k$ to be constant per cell, as we can not
compute the spatial derivatives easily.

For the implementation, $\frac{\partial T}{\partial t}$ is replaced by the
BDF2 approximation, and its terms from older timesteps and $-F$, are moved to
the right-hand side of the PDE.

We use the parameter design presented in {cite:t}`john:knobloch:2006` for
$\delta_K$:
```{math}
\delta_K = \frac{h}{2d\|\mathbf{\beta}\|_{\infty,K}} \left( \coth(Pe)-\frac{1}{Pe} \right)
```
where the Peclet number is given by
```{math}
Pe = \frac{ h \| \mathbf{\beta} \|_{\infty,K}}{2 d k_\text{max}},
```
$d$ is the polynomial degree of the temperature or composition element (typically 2),
$\coth(x) = (1+\exp(-2x)) / (1-\exp(-2x)),$ and
$k_\text{max}=\| k \|_{\infty, K}$ is the maximum conductivity in the cell $K$.

If $Pe<1$, the equation is diffusion-dominated and no stabilization is needed,
so we set $\delta_K=0$. Care needs to be taken in the definition if
$\| \beta \|$ or $k$ become zero:

1.  If $k$ is zero, then $Pe=\infty$ and the right part of the product in the
    definition of $\delta_K$ is equal to one.

2.  If $\| \beta \|$ is zero, $Pe < 1$, so we set $\delta_K=0$.

3.  If both are zero, no stabilization is needed (the field remains constant).