# Adaptive Refinement in Advection-Diffusion Problems by Anomaly Detection: A Numerical Study
## 2 Problem Formulation  
We consider scalar advection-diffusion-reaction problems that, in general, can be modeled as:  
\[\left\{\begin{array}{ll}-div(\varepsilon\
abla u)+\mathbf{b}\cdot\
abla u+ \alpha u=f&\text{in }\Omega\\ u=u_{D}&\text{on }\Gamma_{D}\\ \varepsilon\frac{\partial u}{\partial n}=g&\text{on }\Gamma_{N},\end{array}\right. \tag{1}\]  
where \(\Omega\) in \(\mathbb{R}^{2}\) is a Lipschitz domain with a polygonal boundary \(\Gamma\). In particular, \(\Gamma\) consists of two non-overlapping components \(\Gamma_{D}\) e \(\Gamma_{N}\), such that \(\Gamma=\bar{\Gamma}_{D}\cup\bar{\Gamma}_{N}\) and \(\Gamma_{D}\cap\Gamma_{N}=\emptyset\). In Equation (1), the function \(u\) is the unknown, while \(\varepsilon\) is the diffusivity coefficient, \(\mathbf{b}=[b_{1},b_{2}]\) is the advective field with \(\
abla\cdot\mathbf{b}=0\), \(\alpha\geq 0\) a.e. in \(\Omega\) is the reaction coefficient, \(f\) is the source and \(\mathbf{n}\) is the unit outward normal vector. We will also assume the following:1. \(0<\varepsilon\ll 1\).
2. \(\mathbf{b}\in W^{1,\infty}(\Omega)^{2}\), \(\alpha\in L^{\infty}(\Omega)\), \(f\in L^{2}(\Omega)\) and \(g\in L^{2}(\Gamma_{N})\).
3. There exist two constants, \(\beta\geq 0\) and \(c_{b}\geq 0\), independent of \(\varepsilon\), such that \[-\frac{1}{2}div\mathbf{b}+\alpha\geq\beta\text{ and }\|\ \alpha\ \|_{L^{\infty}}\leq c_{b}\beta.\]  
As already mentioned in the Introduction, scalar advection-diffusion equations describe the _transport_ of scalar quantities, such as temperature or concentration. This transport term is composed of a diffusive part and an advective part.  
From the mathematical proprieties of these two components, it is possible to distinguish regular or exponential layers, parabolic layers and internal or boundary layers. The _accurate_ computation of these singularity regions is crucial to obtain a reliable approximation of the numerical solution of advection-diffusion problems. Standard finite element methods may lead to numerical solutions with nonphysical oscillations, due to the lack of _stability_ of the considered method. (Here, _stability_ is not meant in the sense of the continuous dependence of the solution on the initial data, but rather in the sense of reducing the oscillations occurring in the numerical solution).  
Setting  
\[V:=H_{D}^{1}(\Omega)=\{\varphi\in H^{1}(\Omega):\varphi=0\text{ on }\Gamma_{D}\}, \tag{2}\]  
the weak formulation of problem (1) is:  
\[\text{Find }u\in V:a(u,v)=\ell(v)\quad\forall v\in V, \tag{3}\]  
where we introduced the bilinear form \(a:V\times V\to\mathbb{R}\)  
\[a(u,v)=\int_{\Omega}\varepsilon\
abla u\cdot\
abla v\,d\Omega+\int_{\Omega}v \mathbf{b}\cdot\
abla u\,d\Omega+\int_{\Omega}\alpha uv\,d\Omega \tag{4}\]  
and the functional \(\ell:V\to\mathbb{R}\):  
\[\ell(v)=\int_{\Omega}f\,v\,d\Omega+\int_{\Gamma_{N}}g\,v\,d\gamma. \tag{5}\]  
Let \(\{\mathcal{T}_{h}\}_{h}\) be a family of conforming triangulations of \(\Omega\) into triangles \(T\) with diameter \(h_{T}\leq h\), with \(h=\max\limits_{T\in\mathcal{T}_{h}}h_{T}\). Let \(\mathbb{P}_{m}\) be the space of polynomials of degree at most \(m\), and we set  
\[V_{h}:=\big{\{}\phi\in C(\bar{\Omega})\ |\ \phi|_{T}\in\mathbb{P}_{m}\,\forall T \in\mathcal{T}_{h},\ \phi=0\text{ on }\Gamma_{D}\big{\}}.\]  
By projecting in the space \(V_{h}\subset V\) of continuous piecewise linear finite elements and stabilizing with the SUPG method [1], we can write the generalized formulation:  
\[\text{Find }u_{h}\in V_{h}:a_{h}(u_{h},v_{h})=\ell_{h}(v_{h})\quad\forall v_{h} \in V_{h}, \tag{6}\]  
where  
\[a_{h}(u_{h},v_{h})=a(u_{h},v_{h})+b_{h}(u_{h},v_{h})\quad\text{and }\ell_{h}(v_{h})=\ell(v_{h})+m_{h}(v_{h}).\]  
The terms \(b_{h}(u_{h},v_{h})\) and \(m_{h}(v_{h})\) are the _stabilizing terms_ and they are necessary to reduce the possible numerical oscillations due to the classical [PERSON] formulation.