Project 1: Discontinuous Galerkin methods for advection-diffusion-reaction problems

Project 1: Discontinuous Galerkin methods for advection-diffusion-reaction problems#

In this project you are asked to formulate, analyze, implement and test a discontinuous Galerkin based solver for the diffusion-advection problem

(58)#\[\begin{split} \begin{alignedat}{3} -\epsilon \Delta u_{\epsilon} + \boldsymbol{b} \cdot \nabla u_{\epsilon} + c u_{\epsilon}&= f &&\quad \text{in } \Omega \\ u_{\epsilon} &= u_D &&\quad \text{on } \partial \Omega. \end{alignedat} \end{split}\]

where \(\epsilon > 0\) is the diffusivity parameter, \(\boldsymbol{b} \in \mathbb{R}^2\) is the advection field, \(c \in \mathbb{R}\) is the reaction term, and \(f\) is a source term. The boundary condition \(u_D\) is a given Dirichlet boundary condition on the boundary \(\partial \Omega\) of the domain \(\Omega\).

Part 1: The symmetric interior penalty method#

For this part, set \(\boldsymbol{b} = 0\) in (58).

Task 1: Local conservation properties of the symmetric interior penalty method#

1) Show that the solution \(u\) to the diffusion problem (58) with \(\epsilon = 1\), \(\boldsymbol{b}=0\) and \(c=0\) satisfies the conservation property

(59)#\[ - \int_{\partial U} \partial_n u = \int_U f \]

for any open domain \(U \subset \Omega\). On each mesh facet \(F \in \mathcal{F}\), define the exact flux \(\Phi_F(u)\) by

\[ \Phi_F(u) = \partial_{n_F} u = -\nabla u \cdot \boldsymbol{n}_F \]

and introduce the “facet normal orientation function”

\[ \varepsilon_{T,F}(x) = \boldsymbol{n}_{\partial T}(x) \cdot \boldsymbol{n}_F(x) \]

so that

\[ \varepsilon_{T,F}(x) \boldsymbol{n}_F(x) = \boldsymbol{n}_{\partial T}(x) \]

holds. In particular, (59) means that \(u\) satisfies

\[ \sum_{F\in \mathcal{F}(T)} \varepsilon_{T,F} \Phi_F(u) = \int_T f. \]

where \(\mathcal{F}(T) = \{F \in \mathcal{F}_h \mid F \subset \partial T \}.\)

2) Now, consider the corresponding symmetric interior penalty and try to find a discrete flux \(\phi_F(u_h)\) which is uniquely defined for each facet such

\[ \sum_{F\in \mathcal{F}(T)} \varepsilon_{T,F} \phi_F(u_h) = \int_T f. \]

holds.

Hints

It will turn out that the “natural” discrete \(\phi_F(u_h)\) is not equal \(\Phi_F(u_h) = -\partial_{n_F} u_h\)!
Find the correct flux function by testing the SIP formulation with \(v_h = \chi_T \in \mathbb{P}^k_{\mathrm{dc}}\). Which terms in the bilinear form \(a_h\) will not vanish?

Task 2: Implementation and testing of the symmetric interior penalty method#

1) Implement the DG formulation in either Gridap or ngsolve which solves (58) on the domain \(\Omega = [0,1]^2\) for \(\mathbf{b} = 0 \) but arbitrary positive constants \(\epsilon\) and \(c\). Set \(c=1\).

Use the technique of manufactured solutions to verify your implementation as follows: Generate a series of meshes \(\{\mathcal{T}_i\}_{k=0}^{N}\) with \(N \geqslant 3\) and maximum mesh sizes \(h_{k} = 0.2 \cdot 2^{-k}\). On each mesh you solve the diffusion-reaction problem numerically and compute the norm of the error \(e_k := u-u_k\) in the following norms

\( \epsilon^{1/2}\|\nabla \cdot\|_{\Omega}\)
\( c^{1/2} \|\cdot\|_{\Omega} \)
\( \epsilon^{1/2}\|h^{-1/2} [\cdot] \|_{\mathcal{F}_h^i}\)
\( \epsilon^{1/2}\|h^{-1/2} (u_k - u_D) \|_{\mathcal{F}_h^b}\)

Here, \(u_k\) refers to the discrete solution \(u_h\) on mesh \(\mathcal{T}_k\).

For each computed error norm \(E_k\), calculate the experimental order of convergence defined by

\[ EOC(k) = \dfrac{\log(E_{k-1}/E_k)}{\log(h_{k-1}/h_k)} \]

by either tabulating or plotting the error norm contribution \(E_k\) as function of the mesh size in a \(\log\)-\(\log\) plot. Do such a convergence study for \(\epsilon \in \{ 1, 10^{-1}, 10^{-3}, 10^{-5} \}\).

Hint:
Both ngsolve and gridap provide examples for various DG formulations, see the DG tutorial for ngsolve and the DG-tutorial for Gridap.

2) Next we use our new solver to investigate what happens if we increase the Damköhler number

(60)#\[ \mathrm{Da} = \dfrac{\|c\|_{L^{\infty}(\Omega)} C_P^2}{\epsilon} \]

a) Set \(f = 0.5\) and \(u_D = 0\) and generate a unit square with maxh = 0.1.

Now solve the problem numerically using your SIP implementation for \(\epsilon \in \{ 1, 10^{-1}, 10^{-3}, 10^{-5} \}\). What do you observe, in particular with respect to the fulfillment of the boundary condition?

b) Now repeat the experiment using a standard FEM solver diffusion-reaction where the boundary conditions are built into the function spaces. Make your live easy, and just swap the function space in SIP solver and use

V = H1(mesh, order=order, dgjumps=True, dirichlet=".*")

This is of course not the most efficient implementation as the code will assembles a lot of zeros and reserve a larger sparsity pattern for the system matrix, but effectively, you will obtain a classical FEM solver.

What do you observe now? Can you explain it?

c) Now simply remove the argument dirichlet=".*" in your implementation from b) What do you observe now? Can you explain it?

Part 2: A discontinuous Galerkin method for the advection-diffusion-reaction problem#

Task 1: Derivation and analysis#

Formulate a discontinuous Galerkin formulation for the advection-diffusion-reaction problem (58) which combines the symmetric interior penalty (SIP) method for elliptic problems with the upwind-flux/stabilized discontinuous Galerkin we considered for the pure advection-reaction problem. Give a short justification of your final formulation.

Based on the a priori error estimates we derived individually for the SIP method and upwind flux formulation, what kind of a priori estimate do you expect for your formulation?

Prove the formulated a priori error estimate from Task 1 by combining the approaches we used to prove a priori error estimates for the SIP and the upwind flux formulation.

Hints: Be clever and don’t repeat step by step the entire theoretical analysis we developed for the DG methods previously, but rather try combine them:

To prove that the resulting bilinear form for the advection-diffusion-reaction is coercive with respect to a certain norm, use the discrete coercivity results we derived for SIP and for upwind flux formulation separately to 1) find the right norms form the combined bilinear form and 2) to give a quick proof of its coercivity.
To establish an a priori error estimate, use the “orthogonal subscale” argument to handle the advection-reaction part in the total bilinear form.

Task 2: Implementation and testing#

Implement the DG formulation in either Gridap or ngsolve which solves (58) on the domain \(\Omega = [0,1]^2\). In your implementation allow for the possibility to switch between the upwind-flux and central-flux formulations in the advective part.

Use the manufactured solution

\[\begin{align*} u_{\epsilon} = \cos(\pi x)\dfrac{1-e^{(y-1)/\epsilon}}{1-e^{-2/\epsilon}} + 0.5 \cos(\pi x)\sin(\pi y) \quad \text{in } \Omega \end{align*}\]

with the velocity field \(\boldsymbol{b} = (0,1)\) and \(c = 0.1\) to conduct the following convergence studies for both the upwind-flux and the central flux formulation:

Generate a series of meshes \(\{\mathcal{T}_i\}_{k=0}^{N}\) with \(N \geqslant 3\) and maximum mesh sizes \(h_{k} = 0.2 \cdot 2^{-k}\). On each mesh you solve the advection-diffusion-reaction problem numerically and compute the norm of the error \(e_k := u-u_k\) in the following norms

\( \epsilon^{1/2}\|\nabla \cdot\|_{\Omega}\)
\( c_0^{1/2} \|\cdot\|_{\Omega} \)
\( \epsilon^{1/2}\|h^{-1/2} [\cdot] \|_{\mathcal{F}_h^i}\)
\( \| |\boldsymbol{b}\cdot\boldsymbol{n}_F|^{1/2} [\cdot] \|_{\mathcal{F}_h}\)
\(\| (\tfrac{h}{b_c})^{1/2} \boldsymbol{b} \cdot \nabla (\cdot) \|_{\Omega}\)

Here, \(u_k\) refers to the discrete solution \(u_h\) on mesh \(\mathcal{T}_k\). Note, that we did not prove any error estimate for scaled stream-line derivative norm \(\| h^{1/2}/b_c \boldsymbol{b} \cdot \nabla (\cdot) \|_{\Omega}\), but one can prove that the solution \(u_h \in \mathbb{P}^k_{\mathrm{dc}}(\mathcal{T}_h)\) to the upwind-flux formulation also satisfies

\[ \| (\tfrac{h}{b_c})^{1/2} \boldsymbol{b} \cdot \nabla (u-u_h) \|_{\Omega} \lesssim h^{k+1/2} |u|_{k+1, \Omega} \]