Numerical experiment of the guaranteed local error estimation.

In this project, we give numerical results shown in JCAM article(https://www.sciencedirect.com/science/article/pii/S0377042723000055).

Theoretical discussion

Through this note, the following problem will be considered. $$ (D):\quad -\Delta u = f \mbox{ in } \Omega, \quad \frac{\partial u}{ \partial \mathbf{n}} = g_N \mbox{ on } \Gamma_N, \quad u = g_D \mbox{ on } \Gamma_D. $$ Here, $\Gamma_N$ and $\Gamma_D$ are disjoint subsets of $\partial \Omega$ satisfying $\Gamma_N \cup \Gamma_D = \partial \Omega$; $\mathbf{n}$ is the unit outer normal direction on the boundary and $\frac{\partial }{\partial \mathbf{n}}$ is the directional derivative along $\mathbf{n}$ on $\partial \Omega$.

Notation: - $\mathcal{T}^h$ : the regular triagulation of $\Omega$.

• $h_K, h$ : Diameter of the trangle $K$ and the maximum of $h_K$ for $K \in \mathcal{T}^h$, respectively.

• $C_0, C(h)$ :$L^2$ projection error constant into the piecewise constant space and error constant of the Galerkin projection(computed by the hypercircle), respectively.

• $S, \alpha$ : Intersested subdomain and the cutoff function which support $\Omega'$ is including $S$ , respectively.

• $u_h, \mathbf{p}_h$ : The CG and RT FEM solutions of $(D)$, respectively.

Theorem (Theorem 3 in the article) The local and global error estimators $\widehat{E}_L, \widehat{E}_G$ are given as follows: $$ \widehat{E}_L = \left \{ E_1^2 + E_2^2 \right \}^{\frac{1}{2}} + \mbox{Osc}(f), \quad \widehat{E}_G = |\nabla u_h -\mathbf{p}_h| + C_0 h |f-\pi_0 f |. $$

Here, $$ \kappa_h := \max_{f_h \in X_h} ~~ \min_{ \substack{v_h \in V_{h,0} ,~ {\mathbf{q}}_h \in RT_{h,0}, \ \text{div } {\mathbf{q}}_h + f_h =0} } \frac{| \nabla v_h -{\mathbf{q}}_h|}{|f_h|}, \:C(h) = \sqrt{\kappa_h^2 +C_0^2 h^2}, \ \mbox{Osc}(f) := 0.261 \left \{ \sum_{K \in \mathcal{T}^h } h_K^2 |f-\pi_0f|_{K}^2 \right \}^{\frac{1}{2}}, \ E_1 := |\nabla u_h - \mathbf{p}_h|_{\alpha} +\mbox{Osc}(f) , \quad E_2 := 2\sqrt{2} C(h)\: |\nabla\alpha|_{L^{\infty}(\Omega)} \: |\nabla u_h - \mathbf{p}_h|_{\Omega}. $$

Demonstration

How to run

We summarized codes and implementation in Numerical_Example.ipynb. Press Shift + Enter then you can run codes in the cell.

Codes

The main notebook of this project refers the following python source codes. If you want to know how to use these codes, please see below of line if __name__ == '__main__':.

1. PoissonSolvers.py: Fundametal solvers for Poisson's BVP.
2. ErrorEstimators.py: Global error estimators (including $\kappa_h$ computation, Child-class of PoissonSolvers.py).
3. LocalErrorEstimator.py: This module provides proposed local error estimators (Child-class of ErrorEstimators.py).
4. WeightFunction2D.py : This module provides weight function (α in the article).
5. SingularRHS.py : Singular right hand side data appearing in the example over L-shape domain.
6. MyLocalRefinement.py It refines the part $S=[0.375, 0.625]^2$ of the square mesh.