Advection Diffusion Reaction Equations

$\newcommand{\bm}[1]{\boldsymbol{\mathbf{#1}}} \newcommand{\bvec}[1]{\mathbf{#1}} \newcommand{\bx}{\mathbf{x}} \newcommand{\bu}{\mathbf{u}} \newcommand{\bb}{\mathbf{b}} \newcommand{\bmv}{\mathbf{v}} \newcommand{\bmw}{\mathbf{w}} \newcommand{\Nh}{N_h} \newcommand{\bU}{\mathbf{U}} \newcommand{\bp}{\bvec{p}} \newcommand{\VV}{\mathbf{V}} \newcommand{\WW}{\mathbf{W}} \newcommand{\QQ}{\mathbf{Q}} \renewcommand{\AA}{\mathbf{A}} \newcommand{\BB}{\mathbf{B}} \newcommand{\FF}{\mathbf{F}} \newcommand{\GG}{\mathbf{G}} \newcommand{\CC}{\mathbf{C}} \newcommand{\MM}{\mathbf{M}} \newcommand{\XX}{\mathbf{X}} \newcommand{\KK}{\mathbf{K}} \newcommand{\UU}{\mathbf{U}} \newcommand{\II}{\mathbf{I}} \renewcommand{\SS}{\mathbf{S}} \newcommand{\HH}{\mathbf{H}} \newcommand{\DD}{\mathbf{D}} \newcommand{\NN}{\mathbf{N}} \newcommand{\bmg}{\mathbf{g}} \newcommand{\R}{\mathbb{R}} \newcommand{\JJ}{\mathbf{J}} \newcommand{\RR}{\mathbf{R}}$

Steady problems
Unsteady problems
References

Steady problems

We introduce steady advection-diffusion-reaction equations and their finite element approximation as implemented in redbKIT.

Strong formulation

Let $\Omega \subset \mathbb{R}^{d}$ ($d=2,3$) be an open bounded domain with piecewise smooth boundary $\Gamma = \partial \Omega$. The latter is decomposed into Dirichlet ($\Gamma_D$), Neumann ($\Gamma_N$) and Robin ($\Gamma_R$) components such that $\Gamma = \Gamma_D \cup \Gamma_N \cup \Gamma_R$. We consider the following advection-diffusion-reaction PDE:

$\begin{equation} \left\{ \begin{aligned} - \text{div}( \mu \nabla u) + \bb \cdot \nabla u + c u \ &= s & \qquad \quad \text{in } & \Omega \\ u &= g & \text{on } &\Gamma_{D}\\ \mu \nabla u \bm{n} &= r & \text{on } &\Gamma_{N} \\ \mu \nabla u \bm{n} + \alpha u &= \gamma & \text{on } &\Gamma_{R} \\ \end{aligned} \right. \label{eq:ADRstrong} \end{equation}$

where $\mu : \mathbb{R}^d \rightarrow \mathbb{R}$ is the diffusion coefficient, $\bm{b} : \mathbb{R}^d \rightarrow \mathbb{R}^d$ is a given velocity (convective) field, $c : \mathbb{R}^d \rightarrow \mathbb{R}$ is a reaction coefficient, $s: \mathbb{R}^d \rightarrow \mathbb{R}^d$ is a distributed source term, $g(\bx)$, $r(\bx)$, $\alpha(\bx)$, and $\gamma(\bx)$ are given scalar-valued functions; here $\bm{n}$ denotes the outward unit normal on $\partial \Omega$.

Weak formulation

Let us introduce the following functional spaces:

$V_D = \{ v \in H^1(\Omega): v|_{\Gamma_D} = g\},$ $V = \{ v \in H^1(\Omega): v|_{\Gamma_D} = 0\}.$

The weak formulation of \eqref{eq:ADRstrong} reads:

find $u \in V_D $ such that \begin{equation} \label{eq:ADRweak} a(u,v) = F(v) \qquad \forall v \in V, \end{equation} where $$ a(u,v) = \int_{\Omega} \left( \mu \nabla u \cdot \nabla v + \bb \cdot \nabla u v + cu v \right) \,d\Omega + \int_{\Gamma_R} \alpha u v \, d\Gamma $$ and $$ F(v) = \int_{\Omega} s v \, d\Omega + \int_{\Gamma_N} r v \, d\Gamma + \int_{\Gamma_R} \gamma v \, d\Gamma. $$

Finite element approximation

Let us introduce a FE partition $\mathcal{T}_h$ of the domain $\Omega$ from which we construct a conforming finite element space $X_h \subset H^1(\Omega)$. We also define $V_h = X_h \cap V$ and $V_{Dh} = X_h \cap V_D$. The FE approximation of \eqref{eq:ADRweak} reads: find $u_h \in V_{Dh}$ such that

$\begin{equation} \label{eq:ADRweakFE} a(u_h,v_h) = F(v_h) \qquad \forall v \in V_h. \end{equation}$

We denote by $\{\varphi_i\}_{i=1}^{\Nh}$ a Lagrangian FE basis for $X_h$; $\bu \in \mathbb{R}^{\Nh}$ is the vector of coefficients in the expansion of $u_h$ with respect to the FE basis functions, i.e.

$u_h = \sum_{j=1}^{N_h} \bu_j \varphi_j.$

The algebraic formulation of \eqref{eq:ADRweakFE} reads:

find $\bu \in \R^{N_h}$ such that $$ \begin{equation} \label{eq:ADRweakFEalg} \AA \bu = \mathbf{f}. \end{equation} $$ where $$ \AA_{ij} = a(\varphi_j, \varphi_i), \qquad \mathbf{f}_i = F(\varphi_i), \qquad 1 \leq i,j \leq N_h. $$

The assembly of the matrix $\AA$ and the right-hand side vector $\mathbf{f}$ is performed by the ADR_Assembler function. Once the system has been assembled, Dirichlet boundary conditions should be imposed; a detailed description of this step can be found here.

The resulting linear system can be solved using the LinearSolver class, which provides a flexible layer for both direct and (preconditioned) iterative solvers.

Unsteady problems

We now consider the time-dependent version of problem \eqref{eq:ADRstrong}.

Strong and weak formulations

The strong form of unsteady advection-diffusion-reaction equations can be stated as follows: find $u : \R^{d} \times [0,T] \rightarrow \R$ such that

$\begin{equation} \left\{ \begin{aligned} \frac{\partial u}{\partial t}- \text{div}( \mu \nabla u) + \bb \cdot \nabla u + c u \ &= s & \qquad \quad \text{in } & \Omega \times (0,T) \\ u &= g & \text{on } &\Gamma_{D} \times (0,T)\\ \mu \nabla u \bm{n} &= r & \text{on } &\Gamma_{N} \times (0,T)\\ \mu \nabla u \bm{n} + \alpha u &= \gamma & \text{on } &\Gamma_{R} \times (0,T) \\ u(\cdot , 0) &= u_0 & \text{in } &\Omega, \end{aligned} \right. \label{eq:ADRtstrong} \end{equation}$

Note that all the problem coeffients ($\mu$, $\bb$, $c$, $s$, $g$, $r$, $\alpha$, $\gamma$) might depend on time $t$. The weak formulation of \eqref{eq:ADRtstrong} reads:

for all $t \in (0,T]$, find $u = u(t) \in V_D$ such that \begin{equation} \label{eq:ADRtweak} \left(\frac{\partial u}{\partial t} , v\right) + a(u,v) = F(v) \qquad \forall v \in V, \end{equation} with $u(0) = u_0$.

Finite element approximation

At the algebraic level, the finite element approximation of problem \eqref{eq:ADRtweak} leads to the following first-order linear dynamical system: given $\bu_0 \in \R^{N_h}$, for all $t \in (0,T]$, find $\bu = \bu(t) \in \R^{N_h}$ such that

\begin{equation} \label{eq:ADRtweakAlg} \MM \frac{d \bu}{d t} + \AA(t) \bu = \mathbf{f}(t), \end{equation}

where the mass matrix $\MM$ is defined as

$\MM_{ij} = \left( \varphi_j, \varphi_i \right), \qquad 1 \leq i,j \leq N_h.$

Time discretization

For the time discretization of \eqref{eq:ADRtweakAlg} we consider the Backward Differentiation Formula (BDF) scheme. To begin with, we partition the time interval $[0,T]$ into $N_t$ subintervals of equal size $\Delta t = T / N_t$ and denote by $t_n = n \Delta t$, for $n = 0 , \dots, N_t$ the discrete time instances. Moreover, we denote by $\bu_n$ the approximation of $\bu$ at time $t_n$. We approximate the time derivative of $u$ as

$\frac{\partial \bu }{\partial t}\bigg|_{t=t_n} \approx \frac{\alpha \bu^{n+1} - \bu^{n,\sigma}}{\Delta t}$

where (limiting ourselves to BDF schemes of order $\sigma = 1,2,3$)

$\bu^{n,\sigma} = \begin{cases} \bu^n, & \text{ if } n \geq 0, \quad \text{for } \sigma =1 \text{ (BDF1)}, \\ 2\bu^n - \tfrac{1}{2}\bu^{n-1}, & \text{ if } n \geq 1, \quad \text{for } \sigma =2 \text{ (BDF2)}, \\ 3\bu^n - \tfrac{3}{2}\bu^{n-1} + \tfrac{1}{3}\bu^{n-1}, & \text{ if } n \geq 2, \quad \text{for } \sigma =3 \text{ (BDF3)}, \end{cases}$

and

$\alpha = \begin{cases} 1, & \text{ for } \sigma =1 \text{ (BDF1)}, \\ \tfrac{3}{2}, & \text{ for } \sigma =2 \text{ (BDF2)}, \\ \tfrac{11}{6}, & \text{ for } \sigma =3 \text{ (BDF3)}. \end{cases}$

We end up with following fully-dicrete scheme: given $\bu^{n}, \dots, \bu^{n+1-\sigma}$, for $n \geq \sigma-1$ find $\bu^{n+1}$ such that

$\begin{equation} \label{eq:ADRtweakAlgT} \left( \frac{\alpha}{\Delta t} \MM + \AA(t_{n+1}) \right) \bu^{n+1} = \frac{1}{\Delta t} \MM \bu^{n, \sigma} + \mathbf{f}(t_{n+1}). \end{equation}$

References

Navier-Stokes equations:

K. Eleda Brenan, S.L. Campbell, and L. R. Petzold. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations. SIAM, Philadelphia, 1996.
A. Quarteroni and A. Valli. Numerical Approximation of Partial Differential Equations. Springer-Verlag, Berlin-Heidelberg, 1994.
A. Quarteroni, R. Sacco, and F. Saleri. Numerical Mathematics. Springer, New York, second edition, 2007.