Skip to content

Commit ad3df66

Browse files
sophy1029sophy1029
sophy1029
authored and
sophy1029
committed
Add step-61.
This tutorial program adds a code that solves the Poisson equation using the weak Galerkin formulation.
1 parent d68d619 commit ad3df66

15 files changed

+1324
-0
lines changed

‎examples/step-61/CMakeLists.txt

+39
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,39 @@
1+
##
2+
# CMake script for the step-61 tutorial program:
3+
##
4+
5+
# Set the name of the project and target:
6+
SET(TARGET "step-61")
7+
8+
# Declare all source files the target consists of. Here, this is only
9+
# the one step-X.cc file, but as you expand your project you may wish
10+
# to add other source files as well. If your project becomes much larger,
11+
# you may want to either replace the following statement by something like
12+
# FILE(GLOB_RECURSE TARGET_SRC "source/*.cc")
13+
# FILE(GLOB_RECURSE TARGET_INC "include/*.h")
14+
# SET(TARGET_SRC ${TARGET_SRC} ${TARGET_INC})
15+
# or switch altogether to the large project CMakeLists.txt file discussed
16+
# in the "CMake in user projects" page accessible from the "User info"
17+
# page of the documentation.
18+
SET(TARGET_SRC
19+
${TARGET}.cc
20+
)
21+
22+
# Usually, you will not need to modify anything beyond this point...
23+
24+
CMAKE_MINIMUM_REQUIRED(VERSION 2.8.12)
25+
26+
FIND_PACKAGE(deal.II 9.0.0 QUIET
27+
HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR}
28+
)
29+
IF(NOT ${deal.II_FOUND})
30+
MESSAGE(FATAL_ERROR "\n"
31+
"*** Could not locate a (sufficiently recent) version of deal.II. ***\n\n"
32+
"You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n"
33+
"or set an environment variable \"DEAL_II_DIR\" that contains this path."
34+
)
35+
ENDIF()
36+
37+
DEAL_II_INITIALIZE_CACHED_VARIABLES()
38+
PROJECT(${TARGET})
39+
DEAL_II_INVOKE_AUTOPILOT()

‎examples/step-61/doc/builds-on

+1
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
step-51 step-7 step-20

‎examples/step-61/doc/intro.dox

+279
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,279 @@
1+
<br>
2+
3+
<i>
4+
This program was contributed by Zhuoran Wang.
5+
</i>
6+
7+
<a name="Intro"></a>
8+
<h1>Introduction</h1>
9+
10+
This tutorial program presents an implementation of the "weak Galerkin"
11+
finite element method for the Poisson equation. In some sense, the motivation for
12+
considering this method starts from the same point as in step-51: We would like to
13+
consider discontinuous shape functions, but then need to address the fact that
14+
the resulting problem has a large number of degrees of freedom (because, for
15+
example, each vertex carries as many degrees of freedom as there are adjacent cells).
16+
We also have to address the fact that, in general, every degree of freedom
17+
on one cell couples with all of the degrees of freedom on each of its face neighbor
18+
cells. Consequently, the matrix one gets from the "traditional" discontinuous
19+
Galerkin methods are both large and relatively dense.
20+
21+
Both the hybridized discontinuous Galerkin method (HDG) in step-51 and the weak
22+
Galerkin (WG) method in this tutorial address the issue of coupling by introducing
23+
additional degrees of freedom whose shape functions only live on a face between
24+
cells (i.e., on the "skeleton" of the mesh), and which therefore "insulate" the
25+
degrees of freedom on the adjacent cells from each other: cell degrees of freedom
26+
only couple with other cell degrees of freedom on the same cell, as well as face
27+
degrees of freedom, but not with cell degrees of freedom on neighboring cells.
28+
Consequently, the shape functions for these cell degrees of freedom are
29+
discontinuous and only "live" on exactly one cell.
30+
31+
For a given equation, say the second order Poisson equation,
32+
the difference between the HDG and the WG method is how exactly one formulates
33+
the problem that connects all of these different shape functions. The HDG does
34+
things by reformulating second order problems in terms of a system of first
35+
order equations and then conceptually considers the face degrees of freedom
36+
to be "fluxes" of this first order system. In contrast, the WG method keeps things
37+
in second order form and considers the face degrees of freedom as of the same
38+
type as the primary solution variable, just restricted to the lower-dimensional
39+
faces. For the purposes of the equation, one then needs to somehow "extend"
40+
these shape functions into the interior of the cell when defining what it means
41+
to apply a differential operator to them. Compared to the HDG, the method
42+
has the advantage that it does not lead to a proliferation of unknowns due
43+
to rewriting the equation as a first-order system, but it is also not quite
44+
as easy to implement. However, as we will see in the following, this
45+
additional effort is not prohibitive.
46+
47+
48+
<h3> Weak Galerkin finite element methods</h3>
49+
50+
Weak Galerkin Finite Element Methods (WGFEMs) use discrete weak functions
51+
to approximate scalar unknowns and discrete weak gradients to
52+
approximate classical gradients.
53+
It was introduced by Junping Wang and Xiu Ye
54+
in the paper: <i>A weak Galerkin finite element method for second order elliptic problems</i>,
55+
J. Comput. Appl. Math., 2013, 103-115.
56+
Compared to the continuous Galerkin method,
57+
the weak Galerkin method satisfies important physical properties, namely
58+
local mass conservation and bulk normal flux continuity.
59+
It results in a SPD linear system, and expected convergence rates can
60+
be obtained with mesh refinement.
61+
62+
63+
<h3> WGFEM applied to the Poisson equation </h3>
64+
This program solves the Poisson equation
65+
using the weak Galerkin finite element method:
66+
@f{eqnarray*}
67+
\nabla \cdot \left( -\mathbf{K} \nabla p \right)
68+
= f,
69+
\quad \mathbf{x} \in \Omega, \\
70+
p = p_D,\quad \mathbf{x} \in \Gamma^D, \\
71+
\mathbf{u} \cdot \mathbf{n} = u_N,
72+
\quad \mathbf{x} \in \Gamma^N,
73+
@f}
74+
where $\Omega \subset \mathbb{R}^n (n=2,3)$ is a bounded domain.
75+
In the context of the flow of a fluid through a porous medium,
76+
$p$ is the pressure, $\mathbf{K}$ is a permeability tensor,
77+
$ f $ is the source term, and
78+
$ p_D, u_N $ represent Dirichlet and Neumann boundary conditions.
79+
We can introduce a flux, $\mathbf{u} = -\mathbf{K} \nabla p$, that corresponds
80+
to the Darcy velocity (in the way we did in step-20) and this variable will
81+
be important in the considerations below.
82+
83+
In this program, we will consider a test case where the exact pressure
84+
is $p = \sin \left( \pi x)\sin(\pi y \right)$ on the unit square domain,
85+
with homogenous Dirichelet boundary conditions and identity matrix $\mathbf{K}$.
86+
Then we will calculate $L_2$ errors of pressure, velocity and flux.
87+
88+
89+
<h3> Weak Galerkin scheme </h3>
90+
91+
Via integration by parts, the weak Galerkin scheme for the Poisson equation is
92+
@f{equation*}
93+
\mathcal{A}_h\left(p_h,q \right) = \mathcal{F} \left(q \right),
94+
@f}
95+
where
96+
@f{equation*}
97+
\mathcal{A}_h\left(p_h,q\right)
98+
:= \sum_{T \in \mathcal{T}_h}
99+
\int_T \mathbf{K} \nabla_{w,d} p_h \cdot \nabla_{w,d} q \mathrm{d}x,
100+
@f}
101+
and
102+
@f{equation*}
103+
\mathcal{F}\left(q\right)
104+
:= \sum_{T \in \mathcal{T}_h} \int_T f \, q^\circ \mathrm{d}x
105+
- \sum_{\gamma \in \Gamma_h^N} \int_\gamma u_N q^\partial \mathrm{d}x,
106+
@f}
107+
$q^\circ$ is the shape function of the polynomial space in the interior,
108+
$q^\partial$ is the shape function of the polynomial space on faces, $ \nabla_{w,d} $ means the discrete weak gradient used to approximate the classical gradient.
109+
We use FE_DGQ as the interior polynomial space,
110+
FE_FaceQ as the face polynomial space, and Raviart-Thomas elements for the velocity
111+
$\mathbf{u} = -{\mathbf{K}} \nabla_{w,d} p$.
112+
113+
<h3> Assembling the linear system </h3>
114+
115+
First, we solve for the pressure.
116+
We collect two local spaces together in one FESystem,
117+
the first component in this finite element system denotes
118+
the space for interior pressure, and the second denotes
119+
the space for face pressure.
120+
For the interior component, we use the polynomial space FE_DGQ.
121+
For the face component, we use FE_FaceQ.
122+
123+
We use shape functions defined on spaces FE_DGQ and FE_FaceQ to
124+
approximate pressures, i.e., $p_h = \sum a_i \phi_i,$
125+
where $\phi_i$ are shape functions of FESystem.
126+
We construct the local system by using discrete weak gradients of
127+
shape functions of FE_DGQ and FE_FaceQ.
128+
The discrete weak gradients of shape functions $\nabla_{w,d} \phi$ are defined as
129+
$\nabla_{w,d} \phi = \sum_{i=1}^m c_i \mathbf{w}_i,$
130+
where $\mathbf{w}_i$ is the basis function of $RT(k)$.
131+
132+
Using integration by parts, we have a small linear system
133+
on each element $T$,
134+
@f{equation*}
135+
\int_{T} \left(\nabla_{w,d} \phi \right) \cdot \mathbf{w} \mathrm{d}x=
136+
\int_{T^\partial} \phi^{\partial} \left(\mathbf{w} \cdot \mathbf{n}\right) \mathrm{d}x-
137+
\int_{T^\circ} \phi^{\circ} \left(\nabla \cdot \mathbf{w}\right) \mathrm{d}x,
138+
\quad \forall \mathbf{w} \in RT_{[k]}(E),
139+
@f}
140+
141+
@f{equation*}
142+
\sum_{i=1}^m c_i \int_T \mathbf{w}_i \cdot \mathbf{w}_j \mathrm{d}x =
143+
\int_{T^{\partial}} \phi_i^{\partial}
144+
\left(\mathbf{w}_j \cdot \mathbf{n} \right) \mathrm{d}x -
145+
\int_{T^{\circ}} \phi_i^{\circ} \left (\nabla \cdot \mathbf{w}_j \right)\mathrm{d}x,
146+
@f}
147+
which can be simplified to be
148+
@f{equation*}
149+
\mathbf{C}_{E}\mathbf{M}_{E} = \mathbf{F}_{E},
150+
@f}
151+
where $\mathbf{C}_E$ is the matrix with unknown coefficients $c$,
152+
$\mathbf{M}_E$ is the Gram matrix
153+
$\left[ \int_T \mathbf{w}_i \cdot \mathbf{w}_j \right] \mathrm{d}x$,
154+
$\mathbf{F}_E$ is the matrix of right hand side,
155+
$\mathbf{w}$ and $\phi_i^{\circ}$ are in FEValues,
156+
$\phi_i^{\partial}$ is in FEFaceValues.
157+
Then we solve for $\mathbf{C}_E = \mathbf{F}_E \mathbf{M}_E^{-1}$.
158+
Now, discrete weak gradients of shape functions are written as
159+
linear combinations of basis functions of the $RT$ space.
160+
In our code, we name $\mathbf{C}_E$ as <code>cell_matrix_C</code>,
161+
$\mathbf{M}_E$ as <code>cell_matrix_rt</code>,
162+
$\mathbf{F}_E$ as <code>cell_matrix_F</code>.
163+
164+
The components of the local cell matrices $\mathbf{A}$ are
165+
@f{equation*}
166+
\mathbf{A}_{ij} =
167+
\int_{T} \mathbf{K} \nabla_{w,d} \phi_i \cdot \nabla_{w,d} \phi_j \mathrm{d}x.
168+
@f}
169+
From previous steps, we know $\nabla_{w,d} \phi_i = \sum_{k=1}^m c_{ik} \mathbf{w}_k,$
170+
and $\nabla_{w,d} \phi_j = \sum_{l=1}^m c_{jl} \mathbf{w}_l.$
171+
Then combining the coefficients we have calculated, components of $\mathbf{A}$ are calculated as
172+
@f{equation*}
173+
\int_T \sum_{k,l = 1}^{m}c_{ik} c_{jl} \left(\mathbf{K} \mathbf{w}_i \cdot \mathbf{w}_j\right) \mathrm{d}x
174+
= \sum_{k,l = 1}^{m}c_{ik} c_{jl} \int_{T} \mathbf{K} \mathbf{w}_i \cdot \mathbf{w}_j \mathrm{d}x.
175+
@f}
176+
177+
Next, we use ConstraintMatrix::distribute_local_to_global to
178+
distribute contributions from local matrices $\mathbf{A}$ to the system matrix.
179+
180+
In the scheme
181+
$\mathcal{A}_h\left(p_h,q \right) = \mathcal{F} \left( q \right),$
182+
we have system matrix and system right hand side,
183+
we can solve for the coefficients of the system matrix.
184+
The solution vector of the scheme represents the pressure values in interiors and on faces.
185+
186+
<h3> Post-processing and $L_2$-errors </h3>
187+
188+
After we have calculated the numerical pressure $p$,
189+
we use discrete weak gradients of $p$ to calculate the velocity on each element.
190+
191+
On each element the gradient of the numerical pressure $\nabla p$ can be
192+
approximated by discrete weak gradients $ \nabla_{w,d}\phi_i$, so
193+
@f{equation*}
194+
\nabla_{w,d} p_h = \sum_{i} a_i \nabla_{w,d}\phi_i.
195+
@f}
196+
197+
The numerical velocity $ \mathbf{u}_h = -\mathbf{K} \nabla_{w,d}p_h$ can be written as
198+
@f{equation*}
199+
\mathbf{u}_h = -\mathbf{K} \nabla_{w,d} p =
200+
-\sum_{i} \sum_{j} a_ic_{ij}\mathbf{K}\mathbf{w}_j,
201+
@f}
202+
where $c_{ij}$ is the coefficient of Gram matrix,
203+
$\mathbf{w}_j$ is the basis function of the $RT$ space.
204+
$\mathbf{K} \mathbf{w}_j$ may not be in the $RT$ space.
205+
So we need $L_2$-projection to project it back to the $RT$ space.
206+
207+
We define the projection as
208+
$ \mathbf{Q}_h \left( \mathbf{K}\mathbf{w}_j \right) =
209+
\sum_{k} d_{jk}\mathbf{w}_k$.
210+
For any $j$,
211+
$\left( \mathbf{Q}_h \left( \mathbf{Kw}_j \right),\mathbf{w}_k \right)_E =
212+
\left( \mathbf{Kw}_j,\mathbf{w}_k \right)_E.$
213+
So the numerical velocity becomes
214+
@f{equation*}
215+
\mathbf{u}_h = \mathbf{Q}_h \left( -\mathbf{K}\nabla_{w,d}p_h \right) =
216+
-\sum_{i=0}^{4} \sum_{j=1}^{4}a_ib_{ij}\mathbf{Q}_h \left( \mathbf{K}\mathbf{w}_j \right),
217+
@f}
218+
and we have the following system to solve for the coefficients $d_{jk}$,
219+
@f{equation*}
220+
\left[
221+
\begin{matrix}
222+
\left(\mathbf{w}_i,\mathbf{w}_j \right)
223+
\end{matrix}
224+
\right]
225+
\left[
226+
\begin{matrix}
227+
d_{jk}
228+
\end{matrix}
229+
\right]
230+
=
231+
\left[
232+
\begin{matrix}
233+
\left( \mathbf{Kw}_j,\mathbf{w}_k \right)
234+
\end{matrix}
235+
\right].
236+
@f}
237+
$
238+
\left[
239+
\begin{matrix}
240+
d_{jk}
241+
\end{matrix}
242+
\right]
243+
$
244+
is named <code>cell_matrix_D</code>,
245+
$
246+
\left[
247+
\begin{matrix}
248+
\left( \mathbf{Kw}_j,\mathbf{w}_k \right)
249+
\end{matrix}
250+
\right]
251+
$
252+
is named <code>cell_matrix_E</code>.
253+
254+
Then the elementwise velocity is
255+
@f{equation*}
256+
\mathbf{u}_h = -\sum_{i} \sum_{j}a_ic_{ij}\sum_{k}d_{jk}\mathbf{w}_k =
257+
\sum_{k}- \left(\sum_{j} \sum_{i} a_ic_{ij}d_{jk} \right)\mathbf{w}_k,
258+
@f}
259+
where $-\sum_{j} \sum_{i} a_ic_{ij}d_{jk}$ is named
260+
<code>beta</code> in the code.
261+
262+
We calculate the $L_2$-errors of pressure, velocity and flux
263+
by the following formulas,
264+
@f{eqnarray*}
265+
\|p-p_h^\circ\|^2
266+
= \sum_{T \in \mathcal{T}_h} \|p-p_h^\circ\|_{L^2(E)}^2, \\
267+
\|\mathbf{u}-\mathbf{u}_h\|^2
268+
= \sum_{T \in \mathcal{T}_h} \|\mathbf{u}-\mathbf{u}_h\|_{L^2(E)^2}^2,\\
269+
\|(\mathbf{u}-\mathbf{u}_h) \cdot \mathbf{n}\|^2
270+
= \sum_{T \in \mathcal{T}_h} \sum_{\gamma \subset T^\partial}
271+
\frac{|T|}{|\gamma|} \|\mathbf{u} \cdot \mathbf{n} - \mathbf{u}_h \cdot \mathbf{n}\|_{L^2(\gamma)}^2,
272+
@f}
273+
where $| T |$ is the area of the element,
274+
$\gamma$ are faces of the element,
275+
$\mathbf{n}$ are unit normal vectors of each face.
276+
277+
We will extract interior pressure solutions of each cell
278+
from the global solution and calculate the $L_2$ error
279+
by using function VectorTools::integrate_difference.

‎examples/step-61/doc/kind

+1
Original file line numberDiff line numberDiff line change
@@ -0,0 +1 @@
1+
techniques

0 commit comments

Comments
 (0)