Advection problem with artificial diffusion

[spossann]

We aim to discretize the following advection problem:

$$ \partial_t u + a(t, x) \cdot \nabla u = 0. $$

with $\nabla \cdot a(t, x) = 0$. For stabilization we will use artificial diffusion and change the problem to

$$ \partial_t u + a(t, x) \cdot \nabla u - \nabla \cdot (\sigma \nabla u) = 0, $$

where the matrix $\sigma$ has to be chosen carefully. For example, in SUPG (Streamline Upwind Petrov-Galerkin) methods:

$$ \sigma = \tau \ a a^\top, \qquad \tau \sim \frac{h}{2|a|}, $$

where $\tau$ is a stabilization parameter that can be tuned, and $h$ is the grid parameter in the direction of advection.

We formulate weakly: find $u \in H^1$ such that

$$ \int_\Omega \partial_t u v \ dx + \int_\Omega a(t, x) \cdot \nabla u \ v \ dx + \int_\Omega \nabla u \ \sigma \nabla v \ dx = 0 \qquad \forall v \in H^1. $$

Discretize in space with FEEC:

$$ u(t, x) \in H^1 \quad \approx \quad u_h(t, x) = \sum_{i=0}^{N_0-1} u_i(t) \Lambda^0_i(x) \quad \in V_0 \subset H^1 $$

ans similar for $v \in H^1$. The dimension of the space $V_0$ is $N_0$. For the derivative we have

$$ \nabla u_h(t, x) = \sum_{i=0}^{N_1-1} (\mathbb C \mathbf u(t))_i \Lambda^1_i (x) \quad \in V_1 \subset H(curl), $$

where $\mathbf u(t) = (u_i(t))_i \in \mathbb R^{N_0}$, $\mathbb C: \mathbb R^{N_0} \to \mathbb R^{N_1}$ is the gradient matrix (containing only 0 and 1), and $\Lambda^1(x)$ are the basis functions spanning the space $V_1$ of dimension $N_1$. For more details on the Derham sequence see https://struphy-hub.github.io/struphy/sections/subsections/numerics-geomFE.html

Additionally, we shall need the mass matrices

$$ \mathbb M^0_{ij} = \int_\Omega \Lambda^0_i(x) \Lambda^0_j(x) \ dx , \qquad \mathbb M^{1,\sigma}_{ij} = \int_\Omega \Lambda^1_i(x) \ \sigma \ \Lambda^1_j(x) \ dx $$

The space-discrete problem thus reads

$$ \int_\Omega \partial_t u_h v_h \ dx + \int_\Omega a(t, x) \cdot \nabla u_h \ v_h \ dx + \int_\Omega \nabla u_h \ \sigma \nabla v_h \ dx = 0 \qquad \forall v_h \in V_0. $$

The diffusion part of the equation is straight forward, by inserting the above definitions:

$$ \begin{aligned} \int_\Omega \partial_t u_h v_h \ dx \quad &\to \quad \dot{\mathbf u}^\top \mathbb M^0 \mathbf v \\ \int_\Omega \nabla u_h \ \sigma \nabla v_h \ dx \quad &\to \quad \mathbf u^\top \mathbb C^\top \mathbb M^{1,\sigma} \mathbb C \mathbf v \end{aligned} $$

We now have to decide on the projection used in the advection term. To clearly see the problem, let us rewrite it:

$$ \int_\Omega a(t, x) \cdot \nabla u_h \ v_h \ dx = \int_\Omega f \ v_h \ dx \qquad \forall v_h \in V_0, $$

where $f = a(t, x) \cdot \nabla u_h \in H^1$ is not in the discrete space $V_0$, and thus needs to be projected. The natural choice would be $L^2$-projection $\mathcal P_0: H^1 \to V_0$ defined by

$$ \int_\Omega \mathcal P_0(f) \ \Lambda^0_i \ dx = \int_\Omega f \ \Lambda^0_i \ dx \qquad \forall i. $$

In the advection term, this leads to

$$ \int_\Omega f \ v_h \ dx = \int_\Omega \mathcal P_0(f) \ v_h \ dx , $$

because we can substitute any $\Lambda^0_i$ for $v_h$. From the definition of the $L^2$-projection (two lines above) we obtain

$$ \mathcal P_0(f) = \sum_{i=0}^{N_0 - 1} f_i \Lambda^0_i(x) \qquad \textrm{with} \qquad f_i = (\mathbb{M}^0)^{-1}_{ij} L_j(u_h),\qquad \mathbf f = (\mathbb{M}^0)^{-1} \mathbf L(u_h), $$

where $\mathbf L:\mathbb R^{N_0} \to \mathbb R$ is a linear form that depends on $u_h$, defined by

$$ L_j(u_h) := \int_\Omega a(t, x) \cdot \nabla u_h \Lambda^0_j(x) \ dx , \qquad \mathbf L = (L_j)_j. $$

The advection term is thus discretized as follows:

$$ \int_\Omega a(t, x) \cdot \nabla u_h \ v_h \ dx \quad \to \quad \mathbf L^\top(u_h) (\mathbb{M}^0)^{-\top} \mathbb M^0 \mathbf v = \mathbf L^\top(u_h) \mathbf v. $$

In summary, the fully space-discrete problem reads

$$ \dot{\mathbf u}^\top \mathbb M^0 \mathbf v + \mathbf L^\top(u_h) \mathbf v + \mathbf u^\top \mathbb C^\top \mathbb M^{1,\sigma} \mathbb C \mathbf v = 0 \qquad \forall \mathbf v \in \mathbb R^{N_0}. $$

Since this holds for all $\mathbf v \in \mathbb R^{N_0}$, and using the symmetry of the mass matrices, we obtain

$$ \dot{\mathbf u} + (\mathbb{M}^0)^{-1}\mathbf L(u_h) + (\mathbb{M}^0)^{-1} \mathbb C^\top \mathbb M^{1,\sigma} \mathbb C \mathbf u = 0. $$