Overview</A>
</H3>
This example is devoted to the <EM>discontinuous Galerkin method</EM>, or
-in short: DG method. It includes following topics.
+in short: DG method. It includes the following topics.
<UL>
<LI>Discretization of the linear transport equation with the DG method
-<LI>Two different assembling routines of the system matrix based on face terms given
+<LI>Two different assembling routines for the system matrix based on
+ face terms given as a sum of integrals that
<DL COMPACT>
<DT>1.
-<DD>as a sum of integrals over all cell boundaries, and
+<DD>loops over all cell and all their faces, or that
<DT>2.
-<DD>as a sum of integrals over all faces.
+<DD>loops over all faces, whereas each face is treated only once.
</DL>
<LI>Time comparison of the two assembling routines.
</UL>
</DIV>
<BR CLEAR="ALL">
<P></P>
-the inflow boundary of the domain, whereas <I>n</I> denotes the unit outer
-normal to boundary
+the inflow part of the boundary of the domain, whereas <I>n</I> denotes the
+unit outer normal to the boundary
<!-- MATH: $\partial\Omega$ -->
<IMG
WIDTH="27" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
SRC="step-12.data/img6.gif"
- ALT="$\partial\Omega$">
-of the domain. This equation,
-(<A HREF="step-12.html#transport-equation">1</A>), is the conservative version of the
-transport equation already considered in <A HREF="step-9.html#SECTION00000010000000000000">step 9</A> of this tutorial.
+ ALT="$\partial\Omega$">.
+This
+equation, (<A HREF="step-12.html#transport-equation">1</A>), is the conservative version of
+the transport equation already considered in <A HREF="step-9.html#SECTION00000010000000000000">step 9</A> of this tutorial.
+
+<P>
+In particular will consider problem (<A HREF="step-12.html#transport-equation">1</A>) on
+
+<!-- MATH: $\Omega=[0,1]^2$ -->
+<IMG
+ WIDTH="83" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/img7.gif"
+ ALT="$\Omega=[0,1]^2$">
+with
+<!-- MATH: $\beta=\frac{1}{|x|}(-x_2, x_1)$ -->
+<IMG
+ WIDTH="127" HEIGHT="38" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/img8.gif"
+ ALT="$\beta=\frac{1}{\vert x\vert}(-x_2, x_1)$">
+representing a
+circular counterclockwise flow field, <I>f</I>=0 and <I>g</I>=1 on
+
+<!-- MATH: $x\in\partial\Omega_-^1:=[0,0.5]\times\{0\}$ -->
+<IMG
+ WIDTH="194" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/img9.gif"
+ ALT="$x\in\partial\Omega_-^1:=[0,0.5]\times\{0\}$">
+and <I>g</I>=0 on
+<!-- MATH: $x\in
+\partial\Omega_-\setminus \partial\Omega_-^1$ -->
+<IMG
+ WIDTH="121" HEIGHT="37" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/img10.gif"
+ ALT="$x\in
+\partial\Omega_-\setminus \partial\Omega_-^1$">.
<H3><A NAME="SECTION00003000000000000000">
Discretization</A>
</H3>
-To follow the general theory of deriving DG discretizations for purely
+To follow the general paradigm of deriving DG discretizations for purely
hyperbolic equations we first consider the general hyperbolic problem
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG
WIDTH="170" HEIGHT="31"
- SRC="step-12.data/img7.gif"
+ SRC="step-12.data/img11.gif"
ALT="\begin{displaymath}\nabla\cdot F(u)=f \qquad\mbox{in }\Omega,
\end{displaymath}">
</DIV>
<BR CLEAR="ALL">
<P></P>
subject to appropriate boundary conditions. Here <I>F</I> denotes the flux
-function of the equation under consideration that for in our case, see equation,
-(<A HREF="step-12.html#transport-equation">1</A>), is represented by
+function of the equation under consideration that in our case, see
+
equation (<A HREF="step-12.html#transport-equation">1</A>), is represented by
<!-- MATH: $F(u)=\beta u$ -->
<IMG
WIDTH="86" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img8.gif"
+ SRC="step-12.data/img12.gif"
ALT="$F(u)=\beta u$">.
-For
-deriving the DG discretization we start with a variational,
+For deriving the DG discretization we start with a variational,
mesh-dependent formulation of the problem,
<BR><P></P>
<DIV ALIGN="CENTER">
<IMG
WIDTH="363" HEIGHT="50"
- SRC="step-12.data/img9.gif"
+ SRC="step-12.data/img13.gif"
ALT="\begin{displaymath}\sum_K\left\{-(F(u),\nabla v)_K+(F(u)\cdot n, v)_{\partial K}\right\}=(f,v)_\Omega,
\end{displaymath}">
</DIV>
<!-- MATH: $(\cdot, \cdot)_K$ -->
<IMG
WIDTH="49" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img10.gif"
+ SRC="step-12.data/img14.gif"
ALT="$(\cdot, \cdot)_K$">
and
<!-- MATH: $(\cdot, \cdot)_{\partial
K}$ -->
<IMG
WIDTH="57" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img11.gif"
+ SRC="step-12.data/img15.gif"
ALT="$(\cdot, \cdot)_{\partial
K}$">
simply denote the integrals over the cell and the boundary of
-the cell, respectively. To discretize the problem, the functions <I>u</I>and <I>v</I> are replaced by discrete functions <I>u</I><SUB><I>h</I></SUB> and <I>v</I><SUB><I>h</I></SUB> that in the
+the cell, respectively. To discretize the problem, the functions <I>u</I>
+and <I>v</I> are replaced by discrete functions <I>u</I><SUB><I>h</I></SUB> and <I>v</I><SUB><I>h</I></SUB> that in the
case of discontinuous Galerkin methods belong to the space <I>V</I><SUB><I>h</I></SUB> of
discontinuous piecewise polynomial functions of some degree <I>p</I>. Due
to the discontinuity of the discrete function <I>u</I><SUB><I>h</I></SUB> on interelement
<!-- MATH: $F(u)\cdot n$ -->
<IMG
WIDTH="65" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img12.gif"
+ SRC="step-12.data/img16.gif"
ALT="$F(u)\cdot n$">
must be replaced by a <EM>numerical
flux</EM> function
<!-- MATH: $H(u_h, \hat u_h, n)$ -->
<IMG
WIDTH="97" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img13.gif"
+ SRC="step-12.data/img17.gif"
ALT="$H(u_h, \hat u_h, n)$">,
where
<!-- MATH: $u_h|_{\partial K}$ -->
<IMG
WIDTH="49" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img14.gif"
+ SRC="step-12.data/img18.gif"
ALT="$u_h\vert _{\partial K}$">denotes the inner trace (w.r.t. the cell <I>K</I>) of <I>u</I><SUB><I>h</I></SUB> and
<!-- MATH: $\hat
u_h|_{\partial K}$ -->
<IMG
WIDTH="49" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img15.gif"
+ SRC="step-12.data/img19.gif"
ALT="$\hat
u_h\vert _{\partial K}$">
the outer trace, i.e. the value of <I>u</I><SUB><I>h</I></SUB> on the
neighboring cell. Furthermore the numerical flux function <I>H</I>, among
-other things, is consistent, i.e.
+other things, must be consistent, i.e.
<BR><P></P>
<DIV ALIGN="CENTER">
<!-- MATH: \begin{displaymath}
<IMG
WIDTH="162" HEIGHT="31"
- SRC="step-12.data/img16.gif"
+ SRC="step-12.data/img20.gif"
ALT="\begin{displaymath}H(u,u,n)=F(u)\cdot n,
\end{displaymath}">
</DIV>
<!-- MATH: $u_h\in V_h$ -->
<IMG
WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img17.gif"
+ SRC="step-12.data/img21.gif"
ALT="$u_h\in V_h$">
such that
<BR><P></P>
<TABLE WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-scheme"> </A><IMG
WIDTH="513" HEIGHT="50"
- SRC="step-12.data/img18.gif"
+ SRC="step-12.data/img22.gif"
ALT="\begin{displaymath}
\sum_K\left\{-(F(u_h),\nabla v_h)_K+(H(u_h,\hat u_h,n), v_h)_{\partial K}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h.
\end{displaymath}"></TD>
<BR CLEAR="ALL"><P></P>
Boundary conditions are realized by replacing <IMG
WIDTH="23" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img19.gif"
+ SRC="step-12.data/img23.gif"
ALT="$\hat u_h$">
on the inflow boundary
<!-- MATH: $\partial\Omega_-$ -->
<IMG
WIDTH="39" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img20.gif"
+ SRC="step-12.data/img24.gif"
ALT="$\partial\Omega_-$">
by the boundary function <I>g</I>.
In the special case of the transport equation
<!-- MATH: \begin{equation}
H(u_h,\hat u_h,n)(x)=\left\{\begin{array}{ll}
(\beta\cdot n\, u_h)(x),&\mbox{for } \beta(x)\cdot n(x)<0,\\
- (\beta\cdot n\, \hat u_h)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0.
+ (\beta\cdot n\, \hat u_h)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0,
\end{array}
\right.
\end{equation} -->
<TABLE WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="flux-transport-equation"> </A><IMG
WIDTH="417" HEIGHT="50"
- SRC="step-12.data/img21.gif"
+ SRC="step-12.data/img25.gif"
ALT="\begin{displaymath}
H(u_h,\hat u_h,n)(x)=\left\{\begin{array}{ll}
(\beta\cdot ...
-...)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0.
+...)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0,
\end{array}\right.
\end{displaymath}"></TD>
<TD WIDTH=10 ALIGN="RIGHT">
</TABLE>
</DIV>
<BR CLEAR="ALL"><P></P>
-Hence, the discontinuous Galerkin scheme for the transport equation
-(<A HREF="step-12.html#transport-equation">1</A>) is given by: find
+where on the inflow part of the cell the value is taken from the
+neighboring cell, <IMG
+ WIDTH="23" HEIGHT="32" ALIGN="MIDDLE" BORDER="0"
+ SRC="step-12.data/img23.gif"
+ ALT="$\hat u_h$">,
+and on the outflow part the value is
+taken from the current cell, <I>u</I><SUB><I>h</I></SUB>. Hence, the discontinuous Galerkin
+scheme for the transport equation (<A HREF="step-12.html#transport-equation">1</A>) is given
+by: find
<!-- MATH: $u_h\in V_h$ -->
<IMG
WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img17.gif"
+ SRC="step-12.data/img21.gif"
ALT="$u_h\in V_h$">
such that for all
<!-- MATH: $v_h\in V_h$ -->
<IMG
WIDTH="62" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img22.gif"
+ SRC="step-12.data/img26.gif"
ALT="$v_h\in V_h$">
-following equation holds:
+following
+equation holds:
<BR><P></P>
<DIV ALIGN="CENTER">
<TABLE WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-transport"> </A><IMG
WIDTH="718" HEIGHT="50"
- SRC="step-12.data/img23.gif"
+ SRC="step-12.data/img27.gif"
ALT="\begin{displaymath}
\sum_K\left\{-(u_h,\beta\cdot\nabla v_h)_K+(\beta\cdot n\, ...
...ght\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
<!-- MATH: $\partial K_-:=\{x\in\partial K, \beta(x)\cdot n(x)<0\}$ -->
<IMG
WIDTH="263" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img24.gif"
+ SRC="step-12.data/img28.gif"
ALT="$\partial K_-:=\{x\in\partial K, \beta(x)\cdot n(x)<0\}$">
denotes
the inflow boundary and
K_-$ -->
<IMG
WIDTH="143" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img25.gif"
+ SRC="step-12.data/img29.gif"
ALT="$\partial K_+=\partial K\setminus \partial
K_-$">
-the outflow part of cell <I>K</I>. This equation will be
-refered to as <EM>first version</EM> of the DG method in the following. We
-note that after integrating by parts again, we end up with: find
-
+the outflow part of cell <I>K</I>. Below, this equation will be
+referred to as <EM>first version</EM> of the DG method. We note that
+after integrating by parts again, we end up with: find
<!-- MATH: $u_h\in V_h$ -->
<IMG
WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img17.gif"
- ALT="$u_h\in V_h$">
-such that
+ SRC="step-12.data/img21.gif"
+ ALT="$u_h\in V_h$">such that
<BR><P></P>
<DIV ALIGN="CENTER">
<!-- MATH: \begin{displaymath}
<IMG
WIDTH="509" HEIGHT="50"
- SRC="step-12.data/img26.gif"
+ SRC="step-12.data/img30.gif"
ALT="\begin{displaymath}\sum_K\left\{(\nabla\cdot\{\beta u_h\},\nabla v_h)_K-(\beta\c...
..._{\partial K_-}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h,
\end{displaymath}">
<!-- MATH: $[u_h]=u_h-\hat u_h$ -->
<IMG
WIDTH="114" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img27.gif"
+ SRC="step-12.data/img31.gif"
ALT="$[u_h]=u_h-\hat u_h$">
denotes the jump of the discrete function
between two neighboring cells, that on the boundary of the domain is defined to be
K\cap\partial \hat K$ -->
<IMG
WIDTH="111" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img28.gif"
+ SRC="step-12.data/img32.gif"
ALT="$\Gamma=\partial
K\cap\partial \hat K$">
between two neighboring cells <I>K</I> and <IMG
WIDTH="21" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img29.gif"
+ SRC="step-12.data/img33.gif"
ALT="$\hat K$">occurs twice: Find
<!-- MATH: $u_h\in V_h$ -->
<IMG
WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img17.gif"
+ SRC="step-12.data/img21.gif"
ALT="$u_h\in V_h$">
such that
<BR><P></P>
<TABLE WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-scheme-faces-long"> </A><IMG
WIDTH="723" HEIGHT="50"
- SRC="step-12.data/img30.gif"
+ SRC="step-12.data/img34.gif"
ALT="\begin{displaymath}
-\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma\left\{(H(u_h,\hat ...
...s\partial\Omega}\right\}=(f,v)_\Omega \quad\forall v_h\in V_h,
<!-- MATH: $u_h\in V_h$ -->
<IMG
WIDTH="63" HEIGHT="31" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img17.gif"
+ SRC="step-12.data/img21.gif"
ALT="$u_h\in V_h$">
such that
<BR><P></P>
<DIV ALIGN="CENTER">
<!-- MATH: \begin{equation}
--\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n), [v_h])_{\Gamma\setminus\partial\Omega}+\sum_\Gamma(H(u_h,\hat u_h,n), v_h)_{\Gamma\cap\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
+-\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n), [v_h])_{\Gamma\setminus\partial\Omega}+(H(u_h,g,n), v_h)_{\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
\end{equation} -->
<TABLE WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-scheme-faces"> </A><IMG
- WIDTH="748" HEIGHT="50"
- SRC="step-12.data/img31.gif"
+ WIDTH="695" HEIGHT="50"
+ SRC="step-12.data/img35.gif"
ALT="\begin{displaymath}
-\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n),...
-...Gamma\cap\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
+...), v_h)_{\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
\end{displaymath}"></TD>
<TD WIDTH=10 ALIGN="RIGHT">
(7)</TD></TR>
</TABLE>
</DIV>
<BR CLEAR="ALL"><P></P>
-Whereas the outer unit normal
-
+Whereas the outer unit normal
<!-- MATH: $n|_{\partial K}$ -->
<IMG
WIDTH="41" HEIGHT="34" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img32.gif"
+ SRC="step-12.data/img36.gif"
ALT="$n\vert _{\partial K}$">
-is uniquely defined this is not so for <IMG
+is uniquely defined
+this is not so for <IMG
WIDTH="24" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img33.gif"
- ALT="$n_\Gamma$">.
-<IMG
- WIDTH="24" HEIGHT="30" ALIGN="MIDDLE" BORDER="0"
- SRC="step-12.data/img33.gif"
+ SRC="step-12.data/img37.gif"
ALT="$n_\Gamma$">
-might be the normal from
-either side of the face. Hence, we need to fix the normal <I>n</I> on the face to
-be one of the two normals and denote the other normal by
--<I>n</I>. This way we get -<I>n</I> in the second face term in
+as the latter might be the normal from
+either side of the face. Hence, we need to fix the normal <I>n</I> on the
+face to be one of the two normals and denote the other normal by -<I>n</I>.
+This way we get -<I>n</I> in the second face term in
(<A HREF="step-12.html#dg-scheme-faces-long">6</A>) that finally produces the minus sign in
the jump [<I>v</I><SUB><I>h</I></SUB>] in equation (<A HREF="step-12.html#dg-scheme-faces">7</A>).
<TABLE WIDTH="100%" ALIGN="CENTER">
<TR VALIGN="MIDDLE"><TD ALIGN="CENTER" NOWRAP><A NAME="dg-transport-gamma"> </A><IMG
WIDTH="744" HEIGHT="50"
- SRC="step-12.data/img34.gif"
+ SRC="step-12.data/img38.gif"
ALT="\begin{displaymath}
-\sum_K(u_h,\beta\cdot\nabla v_h)_K+\sum_\Gamma\left\{(\bet...
...ght\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
</TABLE>
</DIV>
<BR CLEAR="ALL"><P></P>
-what will be refered to as <EM>second version</EM> of the DG method.
+which will be refered to as <EM>second version</EM> of the DG method.
<H3><A NAME="SECTION00004000000000000000">
Implementation</A>
</H3>
As already mentioned at the beginning of this example we will
-implement the assembling of the system matrix in two different ways.
+implement assembling the system matrix in two different ways.
The first one will be based on the first version (<A HREF="step-12.html#dg-transport">5</A>)
of the DG method that includes a sum of integrals over all cell
boundaries
<!-- MATH: $\partial K$ -->
<IMG
WIDTH="31" HEIGHT="16" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img35.gif"
+ SRC="step-12.data/img39.gif"
ALT="$\partial K$">.
This is realized by a loop over all cells and
a nested loop over all faces of each cell. Thereby each inner face
<!-- MATH: $\Gamma=\partial K\cap\partial \hat K$ -->
<IMG
WIDTH="111" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img28.gif"
+ SRC="step-12.data/img32.gif"
ALT="$\Gamma=\partial
K\cap\partial \hat K$">
is treated twice, the first
time when the outer loop treats cell <I>K</I> and the second time when it
treats cell <IMG
WIDTH="21" HEIGHT="21" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img29.gif"
+ SRC="step-12.data/img33.gif"
ALT="$\hat K$">.
This way some values like the shape function
values at quadrature points on faces need to be computed twice.
<P>
-To overcome this overhead and for comparison, we implement the
-assembling of a matrix also in a second and different way. This will
+To overcome this overhead and for comparison, we implement
+assembling of matrix also in a second and different way. This will
be based on the second version (<A HREF="step-12.html#dg-transport-gamma">8</A>) that
includes a sum of integrals over all faces <IMG
WIDTH="15" HEIGHT="15" ALIGN="BOTTOM" BORDER="0"
- SRC="step-12.data/img36.gif"
+ SRC="step-12.data/img40.gif"
ALT="$\Gamma$">.
Here, several
-difficulties occur.
+difficulties occurs.
<UL>
-<LI>As degrees of freedom are associated to cells (and not to faces)
- and as a normal is only defined w.r.t a cell adjacent to the face we
+<LI>As degrees of freedom are associated with cells (and not to faces)
+ and as a normal is only defined w.r.t. a cell adjacent to the face we
cannot simply run over all faces of the triangulation but need to
perform the nested loop over all cells and all faces of each cell
like in the first implementation. This, because in <EM>deal.II</EM>
faces are accessible from cells but not visa versa.
<LI>Due to the nested loop we arrive twice at each face. In order to
- assemble the face terms only once we either need to track which
+ assemble face terms only once we either need to track which
faces we have treated before, or we introduce a simple rule that decides
which of the two adjacent cells the face should be accessed and
treated from. Here, we employ the second approach and define the
following rule:
<DL COMPACT>
<DT>a)
-<DD>If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index.
+<DD>If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index on this level.
<DT>b)
<DD>If the two cells are of different refinement levels we access
and treat the face from the coarser cell.
<DD>neighboring cell is of the same refinement level, and
<DT>4.
<DD>neighboring cell is coarser
-</DL>where `neighboring cell' denotes one of the two cells adjacent to the
-face which is not the current cell. In all these cases the assembling
-of the face terms are either the same or similar to each other. Hence,
-we can implement the assembling of the face terms either by `copy and
-paste' (the lazy way, that revenges not later than the scheme or the
-equation might want to be changed afterwards) or by calling a separate
-function that covers all four cases. To be kind of educational within
-this tutorial we perform the latter approach, of course. We go even
-further and encapsulate this function and everything that is needed
-for assembling the specific equation under consideration within a
-class called <EM>DGTransportEquation</EM>. This class includes objects
-of all equation-specific functions, as there are the <EM>RHS</EM> and
-the <EM>BoundaryValues</EM> class, both derived from the <EM>Function</EM>
+</DL>where the `neighboring cell' and the current cell have the mentioned
+faces in common. In last three cases the assembling of the face terms
+are almost the same. Hence, we can implement the assembling of the
+face terms either by `copy and paste' (the lazy way, whose
+disadvantages come up when the scheme or the equation might want to be
+changed afterwards) or by calling a separate function that covers all
+three cases. To be kind of educational within this tutorial we perform
+the latter approach, of course. We go even further and encapsulate
+this function and everything that is needed for assembling the
+specific equation under consideration within a class called
+<EM>DGTransportEquation</EM>. This class includes objects of all
+equation-specific functions, the <EM>RHS</EM> and the
+<EM>BoundaryValues</EM> class, both derived from the <EM>Function</EM>
class, and the <EM>Beta</EM> class representing the vector field.
Furthermore, the <EM>DGTransportEquation</EM> class comprises member
functions <EM>assemble_face_terms1</EM> and
<EM>assemble_system2</EM> functions of the <EM>DGMethod</EM> class will
be general enough to serve for different DG methods, different
equations, even for systems of equations (!) and, under small
-modifications, for nonlinear problems.
+modifications, for nonlinear problems. Finally we note that the
+program is dimension independent, i.e. after replacing
+<EM>DGMethod<2></EM> by <EM>DGMethod<3></EM> the code runs in 3d.
\begin{document}
\subsection{Overview}
This example is devoted to the \emph{discontinuous Galerkin method}, or
-in short: DG method. It includes following topics.
+in short: DG method. It includes the following topics.
\begin{itemize}
\item Discretization of the linear transport equation with the DG method
-\item Two different assembling routines of the system matrix based on face terms given
+\item Two different assembling routines for the system matrix based on
+ face terms given as a sum of integrals that
\begin{enumerate}
-\item as a sum of integrals over all cell boundaries, and
-\item as a sum of integrals over all faces.
+\item loops over all cell and all their faces, or that
+\item loops over all faces, whereas each face is treated only once.
\end{enumerate}
\item Time comparison of the two assembling routines.
\end{itemize}
\[
\partial\Omega_-:=\{x\in\partial\Omega, \beta(x)\cdot n(x)<0\}
\]
-the inflow boundary of the domain, whereas $n$ denotes the unit outer
-normal to boundary $\partial\Omega$ of the domain. This equation,
-\eqref{transport-equation}, is the conservative version of the
-transport equation already considered in step 9 of this tutorial.
+the inflow part of the boundary of the domain, whereas $n$ denotes the
+unit outer normal to the boundary $\partial\Omega$. This
+equation, \eqref{transport-equation}, is the conservative version of
+the transport equation already considered in step 9 of this tutorial.
+
+In particular will consider problem \eqref{transport-equation} on
+$\Omega=[0,1]^2$ with $\beta=\frac{1}{|x|}(-x_2, x_1)$ representing a
+circular counterclockwise flow field, $f=0$ and $g=1$ on
+$x\in\partial\Omega_-^1:=[0,0.5]\times\{0\}$ and $g=0$ on $x\in
+\partial\Omega_-\setminus \partial\Omega_-^1$.
\subsection{Discretization}
-To follow the general theory of deriving DG discretizations for purely
+To follow the general paradigm of deriving DG discretizations for purely
hyperbolic equations we first consider the general hyperbolic problem
\[
\nabla\cdot F(u)=f \qquad\mbox{in }\Omega,
\]
subject to appropriate boundary conditions. Here $F$ denotes the flux
-function of the equation under consideration that for in our case, see equation,
-\eqref{transport-equation}, is represented by $F(u)=\beta u$. For
-deriving the DG discretization we start with a variational,
+function of the equation under consideration that in our case, see
+equation \eqref{transport-equation}, is represented by $F(u)=\beta u$.
+For deriving the DG discretization we start with a variational,
mesh--dependent formulation of the problem,
\[
\sum_K\left\{-(F(u),\nabla v)_K+(F(u)\cdot n, v)_{\partial K}\right\}=(f,v)_\Omega,
a test function and integration by parts on each cell $K$ of the
triangulation. Here $(\cdot, \cdot)_K$ and $(\cdot, \cdot)_{\partial
K}$ simply denote the integrals over the cell and the boundary of
-the cell, respectively. To discretize the problem, the functions $u$
+the cell, respectively. To discretize the problem, the functions $u$
and $v$ are replaced by discrete functions $u_h$ and $v_h$ that in the
case of discontinuous Galerkin methods belong to the space $V_h$ of
discontinuous piecewise polynomial functions of some degree $p$. Due
denotes the inner trace (w.r.t. the cell $K$) of $u_h$ and $\hat
u_h|_{\partial K}$ the outer trace, i.e. the value of $u_h$ on the
neighboring cell. Furthermore the numerical flux function $H$, among
-other things, is consistent, i.e.
+other things, must be consistent, i.e.
\[
H(u,u,n)=F(u)\cdot n,
\]
\begin{equation}\label{flux-transport-equation}
H(u_h,\hat u_h,n)(x)=\left\{\begin{array}{ll}
(\beta\cdot n\, u_h)(x),&\mbox{for } \beta(x)\cdot n(x)<0,\\
- (\beta\cdot n\, \hat u_h)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0.
+ (\beta\cdot n\, \hat u_h)(x),&\mbox{for } \beta(x)\cdot n(x)\geq 0,
\end{array}
\right.
\end{equation}
-Hence, the discontinuous Galerkin scheme for the transport equation
-\eqref{transport-equation} is given by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following equation holds:
+where on the inflow part of the cell the value is taken from the
+neighboring cell, $\hat u_h$, and on the outflow part the value is
+taken from the current cell, $u_h$. Hence, the discontinuous Galerkin
+scheme for the transport equation \eqref{transport-equation} is given
+by: find $u_h\in V_h$ such that for all $v_h\in V_h$ following
+equation holds:
\begin{equation}\label{dg-transport}
\sum_K\left\{-(u_h,\beta\cdot\nabla v_h)_K+(\beta\cdot n\, u_h, v_h)_{\partial K_+\setminus\partial\Omega}+(\beta\cdot n\, \hat u_h, v_h)_{\partial K_-\setminus\partial\Omega}\right\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
\end{equation}
where $\partial K_-:=\{x\in\partial K, \beta(x)\cdot n(x)<0\}$ denotes
the inflow boundary and $\partial K_+=\partial K\setminus \partial
-K_-$ the outflow part of cell $K$. This equation will be
-refered to as \emph{first version} of the DG method in the following. We
-note that after integrating by parts again, we end up with: find
-$u_h\in V_h$ such that
+K_-$ the outflow part of cell $K$. Below, this equation will be
+referred to as \emph{first version} of the DG method. We note that
+after integrating by parts again, we end up with: find $u_h\in V_h$
+such that
\[
\sum_K\left\{(\nabla\cdot\{\beta u_h\},\nabla v_h)_K-(\beta\cdot n [u_h], v_h)_{\partial K_-}\right\}=(f,v)_\Omega, \quad\forall v_h\in V_h,
\]
By employing conservativity \eqref{conservative} of the numerical flux
this equation simplifies to: find $u_h\in V_h$ such that
\begin{equation}\label{dg-scheme-faces}
- -\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n), [v_h])_{\Gamma\setminus\partial\Omega}+\sum_\Gamma(H(u_h,\hat u_h,n), v_h)_{\Gamma\cap\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
+ -\sum_K(F(u_h),\nabla v_h)_K+\sum_\Gamma(H(u_h,\hat u_h,n), [v_h])_{\Gamma\setminus\partial\Omega}+(H(u_h,g,n), v_h)_{\partial\Omega}=(f,v)_\Omega \quad\forall v_h\in V_h.
\end{equation}
-Whereas the outer unit normal
-$n|_{\partial K}$ is uniquely defined this is not so for $n_\Gamma$. $n_\Gamma$ might be the normal from
-either side of the face. Hence, we need to fix the normal $n$ on the face to
-be one of the two normals and denote the other normal by
-$-n$. This way we get $-n$ in the second face term in
+Whereas the outer unit normal $n|_{\partial K}$ is uniquely defined
+this is not so for $n_\Gamma$ as the latter might be the normal from
+either side of the face. Hence, we need to fix the normal $n$ on the
+face to be one of the two normals and denote the other normal by $-n$.
+This way we get $-n$ in the second face term in
\eqref{dg-scheme-faces-long} that finally produces the minus sign in
the jump $[v_h]$ in equation \eqref{dg-scheme-faces}.
\begin{equation}\label{dg-transport-gamma}
-\sum_K(u_h,\beta\cdot\nabla v_h)_K+\sum_\Gamma\left\{(\beta\cdot n\, u_h, [v_h])_{\Gamma_+\setminus\partial\Omega}+(\beta\cdot n\, \hat u_h, [v_h])_{\Gamma_-\setminus\partial\Omega}\right\}=(f,v)_\Omega-(\beta\cdot n\, g, v_h)_{\partial\Omega_-},
\end{equation}
-what will be refered to as \emph{second version} of the DG method.
+which will be refered to as \emph{second version} of the DG method.
\subsection{Implementation}
As already mentioned at the beginning of this example we will
-implement the assembling of the system matrix in two different ways.
+implement assembling the system matrix in two different ways.
The first one will be based on the first version \eqref{dg-transport}
of the DG method that includes a sum of integrals over all cell
boundaries $\partial K$. This is realized by a loop over all cells and
treats cell $\hat K$. This way some values like the shape function
values at quadrature points on faces need to be computed twice.
-To overcome this overhead and for comparison, we implement the
-assembling of a matrix also in a second and different way. This will
+To overcome this overhead and for comparison, we implement
+assembling of matrix also in a second and different way. This will
be based on the second version \eqref{dg-transport-gamma} that
includes a sum of integrals over all faces $\Gamma$. Here, several
-difficulties occur.
+difficulties occurs.
\begin{itemize}
-\item As degrees of freedom are associated to cells (and not to faces)
- and as a normal is only defined w.r.t a cell adjacent to the face we
+\item As degrees of freedom are associated with cells (and not to faces)
+ and as a normal is only defined w.r.t. a cell adjacent to the face we
cannot simply run over all faces of the triangulation but need to
perform the nested loop over all cells and all faces of each cell
like in the first implementation. This, because in \emph{deal.II}
faces are accessible from cells but not visa versa.
\item Due to the nested loop we arrive twice at each face. In order to
- assemble the face terms only once we either need to track which
+ assemble face terms only once we either need to track which
faces we have treated before, or we introduce a simple rule that decides
which of the two adjacent cells the face should be accessed and
treated from. Here, we employ the second approach and define the
following rule:
\begin{enumerate}
\renewcommand{\labelenumi}{\alph{enumi})}
- \item If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index.
+ \item If the two cells adjacent to a face are of the same refinement level we access and treat the face from the cell with lower index on this level.
\item If the two cells are of different refinement levels we access
and treat the face from the coarser cell.
\end{enumerate}
\item neighboring cell is of the same refinement level, and
\item neighboring cell is coarser
\end{enumerate}
-where `neighboring cell' denotes one of the two cells adjacent to the
-face which is not the current cell. In all these cases the assembling
-of the face terms are either the same or similar to each other. Hence,
-we can implement the assembling of the face terms either by `copy and
-paste' (the lazy way, that revenges not later than the scheme or the
-equation might want to be changed afterwards) or by calling a separate
-function that covers all four cases. To be kind of educational within
-this tutorial we perform the latter approach, of course. We go even
-further and encapsulate this function and everything that is needed
-for assembling the specific equation under consideration within a
-class called \emph{DGTransportEquation}. This class includes objects
-of all equation--specific functions, as there are the \emph{RHS} and
-the \emph{BoundaryValues} class, both derived from the \emph{Function}
+where the `neighboring cell' and the current cell have the mentioned
+faces in common. In last three cases the assembling of the face terms
+are almost the same. Hence, we can implement the assembling of the
+face terms either by `copy and paste' (the lazy way, whose
+disadvantages come up when the scheme or the equation might want to be
+changed afterwards) or by calling a separate function that covers all
+three cases. To be kind of educational within this tutorial we perform
+the latter approach, of course. We go even further and encapsulate
+this function and everything that is needed for assembling the
+specific equation under consideration within a class called
+\emph{DGTransportEquation}. This class includes objects of all
+equation--specific functions, the \emph{RHS} and the
+\emph{BoundaryValues} class, both derived from the \emph{Function}
class, and the \emph{Beta} class representing the vector field.
Furthermore, the \emph{DGTransportEquation} class comprises member
functions \emph{assemble\_face\_terms1} and
\emph{assemble\_system2} functions of the \emph{DGMethod} class will
be general enough to serve for different DG methods, different
equations, even for systems of equations (!) and, under small
-modifications, for nonlinear problems.
+modifications, for nonlinear problems. Finally we note that the
+program is dimension independent, i.e. after replacing
+\emph{DGMethod<2>} by \emph{DGMethod<3>} the code runs in 3d.
\end{document}
projects / dealii-svn.git / commitdiff