SUPG Stabilization#

For streamline upwind/Petrov-Galerkin (SUPG) (see for example Clevenger and Heister [2019], John and Knobloch [2006]), we add to the weak form \(a(\cdot,\cdot)\) the cell-wise defined weak form

\[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 John and Knobloch [2006] for \(\delta_K\):

\[\delta_K = \frac{h}{2d\|\mathbf{\beta}\|_{\infty,K}} \left( \coth(Pe)-\frac{1}{Pe} \right)\]

where the Peclet number is given by

\[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).