Linear operator for a method-of-lines discretization of the unsteady advection-diffusion equation with SUPG (streamline upwind Petrov-Galerkin) stabilization. More...

#include <unsteady_advection_diffusion_supg.hpp>

Public Types
using	SrcVectorType = linalg::VectorQ1Scalar< ScalarT >

using	DstVectorType = linalg::VectorQ1Scalar< ScalarT >

using	ScalarType = ScalarT

Public Member Functions
	UnsteadyAdvectionDiffusionSUPG (const grid::shell::DistributedDomain &domain, const grid::Grid3DDataVec< ScalarT, 3 > &grid, const grid::Grid2DDataScalar< ScalarT > &radii, const grid::Grid4DDataScalar< grid::shell::ShellBoundaryFlag > &boundary_mask, const linalg::VectorQ1Vec< ScalarT, VelocityVecDim > &velocity, const ScalarT diffusivity, const ScalarT dt, bool treat_boundary, bool diagonal=false, ScalarT mass_scaling=1.0, bool lumped_mass=false, linalg::OperatorApplyMode operator_apply_mode=linalg::OperatorApplyMode::Replace, linalg::OperatorCommunicationMode operator_communication_mode=linalg::OperatorCommunicationMode::CommunicateAdditively)

ScalarT &	dt ()

const ScalarT &	dt () const

void	apply_impl (const SrcVectorType &src, DstVectorType &dst)

void	operator() (const int local_subdomain_id, const int x_cell, const int y_cell, const int r_cell) const

Detailed Description

template<typename ScalarT, int VelocityVecDim = 3>
class terra::fe::wedge::operators::shell::UnsteadyAdvectionDiffusionSUPG< ScalarT, VelocityVecDim >

Linear operator for a method-of-lines discretization of the unsteady advection-diffusion equation with SUPG (streamline upwind Petrov-Galerkin) stabilization.

Continuous problem

The unsteady advection-diffusion equation is given by

\[ \frac{\partial}{\partial t}T + \mathbf{u} \cdot \nabla T - \nabla \cdot (\kappa \nabla T) = f \]

where \( T \) is the scalar temperature solution, \( \mathbf{u} \) a given velocity field, \( f \) a given source term, and \( \kappa \) a given diffusivity function.

Space discretization

Note: We assume here that we have \( \kappa|_E = \mathrm{const} \) on each element \(E\) which simplifies the implementation of the SUPG stabilization for linear finite elements as certain terms drop. Currently, \(\kappa = \mathrm{const}\) globally, but once that is changed, we will need to average on every element it in the kernel or pass it in as an elementwise constant (FE) function.

We first discretize in space, then in time (method of lines). After discretization in space, we get the system of ODEs in time

\[ M \frac{d}{dt}T + (C + K + G)T = F + F_\mathrm{SUPG} \]

where

Term	Description	Bilinear form
\(M_{ij}\)	mass	\( \int \phi_i \phi_j \)
\(C_{ij}\)	advection	\( \int \phi_i (\mathbf{u} \cdot \nabla \phi_j) \)
\(K_{ij}\)	diffusion	\( \int \nabla \phi_i \cdot (\kappa \nabla \phi_j) \)
\(G_{ij}\)	SUPG adv-adv	\( \sum_E \int_E \tau_E (\mathbf{u} \cdot \nabla \phi_i) (\mathbf{u} \cdot \nabla \phi_j) \)
\(F_i\)	forcing	\( \int \phi_i f \)
\((F_\mathrm{SUPG})_{i}\)	SUPG forcing	\( \sum_E \int_E \tau_E (\mathbf{u} \cdot \nabla \phi_i) f \).

Note: After brief testing, it seems that in general the term \(F_\mathrm{SUPG}\) does not always improve the computed solutions. It even appears to slightly increase the error sometimes. So setting \(F_\mathrm{SUPG} = 0\) can really be just fine in certain settings.

Time discretization

For the time discretization, we employ implicit BDF schemes. A general formula is

\[ (\alpha_0 M + \Delta t A) T^{n+1} = - \sum_{j=1}^k \alpha_j M T^{n+1-j} + \Delta t R^{n+1} \]

where \(A = (C + K + G)\) and \(R = F + F_{\mathrm{SUPG}}\).

We recover the common BDF1 (backward or implicit Euler) and BDF2 schemes by choosing:

Scheme	\(k\)	\(\alpha\)	full equation
BDF1 (backward/implicit Euler)	1	\([1, -1]\)	\((M + \Delta t A) T^{n+1} = M T^{n} + \Delta t R^{n+1}\)
BDF2	2	\([\frac{3}{2}, -2, \frac{1}{2}]\)	\((\frac{3}{2} M + \Delta t A) T^{n+1} = 2 M T^{n} - \frac{1}{2} M T^{n-1} + \Delta t R^{n+1}\)

The RHS term must be built with appropriate other classes/function. This class is only concerned with the matrix-free evaluation of the LHS system matrix

\[ \alpha_0 M + \Delta t A. \]

The parameters \(\alpha_0\) and \(\Delta t\) are set through the constructor via mass_scaling and dt respectively.

SUPG stabilization

Several choices for the stabilization parameter \( \tau_E \) are possible. As it is commonly selected in the literature (see, e.g., John and Knobloch (2006)), we set on element \(E\)

\[ \tau_E = \frac{h_E}{2 \|\mathbf{u}\|} \left(\coth(\mathrm{Pe_E}) - \frac{1}{\mathrm{Pe_E}}\right) \]

where

\[ \mathrm{Pe_E} = \frac{\|\mathbf{u}\|h_E}{2 \kappa} \]

with some precautionary measures to avoid edge cases that would result in division by zero etc.