<a name="Intro"></a>
<h1>Introduction</h1>
-Please note: This is work in progress and will be an example for block
-smoothers in geometric multigrid.
+This program solves an advection-diffusion problem using a geometric
+multigrid preconditioner. The basics of this preconditioner are
+discussed in step-16; here we discuss the necessary changes needed for
+a non-symmetric PDE, we introduce the idea of block smoothing (as
+compared to point smoothing), and examine the effects of DoF
+renumbering for additive and multiplicative methods.
-OUTLINE
+<h2>Equation</h2>
+The advection-diffusion equation is given by
-- based on step-16 (GMG), and step-9 (advection)
+@f{align*}{
+-\varepsilon \Delta u + \vec{\beta}\cdot \nabla u & =& f &
+\text{ on } \Omega\\
+u &=& g & \text{ on } \partial\Omega
+@f}
-differences to GMG in step-16:
+where $\varepsilon>0$, $\vec{\beta}$ is the \textit{advection
+direction}, and $f$ is a source. A few notes:
-- non-symmetric problem (interface matrices)
-- renumbering of dofs on each level
-- different smoothers (include block)
+1. If $\vec{\beta}=\vec{0}$, this is the Laplace equation solved in
+step-16 (and many other places).
+2. If $\varepsilon=0$ then this is the advection equation solved in
+step-9.
-<h2>Equation</h2>
+3. If $\varepsilon \ll \norm{\vec{\beta}},$ we say the problem is
+\textit{advection-dominated}, else we say the problem is
+\textit{diffusion-dominated}.
+
+For the discussion in this tutorial we will be concerned with
+advection-dominated flow.
+
+Using the standard Galerkin finite element method, for suitable test
+function $v_h$, the discretized weak form of the PDE reads
+
+@f{align*}{
+a(u_h,\v_h) = F(v_h)
+@f}
+
+where
+
+@f{align*}{
+a(u_h,\v_h) &= (\nabla u_h,\, \nabla v_h) + (\vec{\beta}\cdot\nabla u_h,\,v_h)\\
+F(v_h) &= (f,\,v)
+@f}
-* Statement of PDE and weak-form with physical description
+<h3>Streamline diffusion</h3>
-* SUPG (streamline diffusion): State modifications to weak-form.
+The following error estimate can be shown for this PDE:
+
+@f{align*}{
+\norm{\nabla (u-u_h)} \leq (1+\mathcal{P}) \inf_{v_h} \norm (\nabla (u-v_h)}
+@f}
+
+where $\mathcal{P} \sim \norm{\vec{\beta}}/\varepsilon$ (referred
+to as the \textit{Peclet} number). This implies that we may have poor
+numerical solutions when $\varepsilon \ll \norm{\vec{\beta}}$. To
+combat this, we will consider the new weak form
+
+@f{align*}{
+a(u_h,\,v_h) + (-\varepsilon \Delta u_h + \vec{\beta}\cdot
+\nabla u_h-f,\,\delta \vec{\beta}\cdot \nabla v_h) = F(v_h)
+@f}
+
+where $\delta$ is a stabilization parameter defined in <a
+href="https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27">
+On Discontinuity—Capturing Methods for Convection—Diffusion Equations
+by Volker and Petr</a>. Essentially, adding in the discrete strong
+form residual enhances the coercivity of the bilinear form
+$a(\cdot,\cdot)$ which increases the stability of the discrete
+solution. This method is commonly referred to as \textit{streamline
+diffusion} or \textit{SUPG} (streamline upwind/Petrov-Galerkin).
<h2>Smoothers</h2>
One of the goals of this tutorial is to expand from using a simple
-(point-wise) Gauss-Seidel (SOR) smoother that is used in step-16 (class
-PreconditionSOR) on each level of the multigrid hierarchy. Here, we consider
-point-wise smoothers (Jacobi and SOR) and cell-based smoothers (Block Jacobi
-and Block SOR). The cell-based smoothers can best be explained within the framework
-of additive and multiplicative Schwarz methods.
+(point-wise) Gauss-Seidel (SOR) smoother that is used in step-16
+(class PreconditionSOR) on each level of the multigrid hierarchy.
+Here, we consider point-wise smoothers (Jacobi and SOR) and cell-based
+smoothers (Block Jacobi and Block SOR). The cell-based smoothers can
+best be explained within the framework of additive and multiplicative
+Schwarz methods.
In contrast to step-16, our test problem contains an advective
term. Especially with small viscosity, information is transported along
Finally, let us note that the standard Gauss-Seidel (or SOR) method can be
seen as a multiplicative Schwarz method with a subproblem for each DoF.
+
<h2>Test problem</h2>
- Show
- image of solution without and with SUPG.
+We will be considering the following test problem: $\Omega =
+[-1,\,1]\times[-1,\,1]$ with circle of radius 0.3 centered at the
+origin removed, $\varepsilon=0.005$, $\vec{\beta} =
+[-\sin(\pi/6),\,\cos(\pi/6)]$, $f=0$, and the boundary function
+
+@f{align*}{
+g = \left\{\begin{array}{ll} 1 & x=-1 \text{ OR } y=-1,\,x\geq 0.5 \\
+0 & \text{else} \end{array}\right.
+@f}
+
+The following figures depict the solutions with (left) and without
+(right) streamline diffusion. Without streamline diffusion we see high
+oscillations around the boundary layer, demonstrating the instability
+of the standard Galerkin finite elements for this problem.
+
+<TABLE WIDTH="60%" ALIGN="center">
+ <tr>
+ <td ALIGN="center">
+ <img src="https://www.dealii.org/images/steps/developer/step-63-solution.png" alt="">
+ </td>
+
+ <td ALIGN="center">
+ <img src="https://www.dealii.org/images/steps/developer/step-63-solution-no-sd.png" alt="">
+ </td>
+ </tr>
+</table>
able to compute the inverse of the cell matrices much cheaper than
what is currently being done inside deal.II. This research is based on
the fast diagonalization method (dating back to the 1960s) and has
-been used in the spectral community for around 20 years (see, e.g.,
-https://doi.org/10.1007/s10915-004-4787-3). There are currently
-efforts to generalize these methods to DG and make them more
-robust. Also, it seems that one should be able to take advantage of
-matrix-free implementations and the fact that, in the interior of
-the domain, cell matrices tend to look very similar, allowing fewer
-matrix inverse computations.
+been used in the spectral community for around 20 years (see, e.g., <a
+href="https://doi.org/10.1007/s10915-004-4787-3"> Hybrid
+Multigrid/Schwarz Algorithms for the Spectral Element Method by Lottes
+and Fischer</a>). There are currently efforts to generalize these
+methods to DG and make them more robust. Also, it seems that one
+should be able to take advantage of matrix-free implementations and
+the fact that, in the interior of the domain, cell matrices tend to
+look very similar, allowing fewer matrix inverse computations.
Combining 1. and 2. gives a good reason for expecting that a method
like block Jacobi could become very powerful in the future, even
though currently for these examples it is quite slow.
-
-
<h3> Possible Extensions </h3>
<h4> Constant iterations for Q<sub>5</sub> </h4>
// @sect3{Right-hand Side and Boundary Values}
// The problem solved in this tutorial is an adaptation of Ex. 3.1.3
- // found on pg. 118 of Finite Elements and Fast Iterative Solvers:
- // with Applications in Incompressible Fluid Dynamics by Elman, Silvester,
- // and Wathen. The main difference being that we add a hole in the center
- // of our domain with zero Dirichlet boundary.
+ // found on pg. 118 of
+ // <a
+ // href="https://global.oup.com/academic/product/finite-elements-and-fast-iterative-solvers-9780199678808">
+ // Finite Elements and Fast Iterative Solvers: with Applications in
+ // Incompressible Fluid Dynamics by Elman, Silvester, and Wathen</a> with
+ // Applications in Incompressible Fluid Dynamics by Elman, Silvester, and
+ // Wathen. The main difference being that we add a hole in the center of our
+ // domain with zero Dirichlet boundary.
// We have a zero right-hand side.
template <int dim>
// @sect3{Streamline Diffusion}
// Streamline diffusion stabilization term. Value is defined in
- // "On Discontinuity—Capturing Methods for Convection—Diffusion
- // Equations" by Volker and Petr
- // (https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27).
+ // <a href="https://link.springer.com/chapter/10.1007/978-3-540-34288-5_27">
+ // On Discontinuity—Capturing Methods for Convection—Diffusion Equations
+ // by Volker and Petr</a>
template <int dim>
double compute_stabilization_delta(const double hk,
const double eps,
copy_data.cell_matrix(i, j) +=
// Galerkin contribution:
(settings.epsilon *
- scratch_data.fe_values.shape_grad(j, q_point) *
scratch_data.fe_values.shape_grad(i, q_point) *
+ scratch_data.fe_values.shape_grad(j, q_point) *
scratch_data.fe_values.JxW(q_point)) +
- ((advection_direction *
+ (scratch_data.fe_values.shape_value(i, q_point) *
+ (advection_direction *
scratch_data.fe_values.shape_grad(j, q_point)) *
- scratch_data.fe_values.shape_value(i, q_point)) *
- scratch_data.fe_values.JxW(q_point) +
+ scratch_data.fe_values.JxW(q_point)) +
// Streamline diffusion contribution:
delta *
(advection_direction *
projects / dealii.git / commitdiff