]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
First version of intro.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 11 Feb 2006 00:26:02 +0000 (00:26 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Sat, 11 Feb 2006 00:26:02 +0000 (00:26 +0000)
git-svn-id: https://svn.dealii.org/trunk@12310 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/tutorial/chapter-2.step-by-step/step-20.data/intro.tex [new file with mode: 0644]

diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-20.data/intro.tex b/deal.II/doc/tutorial/chapter-2.step-by-step/step-20.data/intro.tex
new file mode 100644 (file)
index 0000000..2d5851b
--- /dev/null
@@ -0,0 +1,366 @@
+\documentclass{article}
+\usepackage{amsmath}
+\renewcommand{\vec}[1]{\mathbf{#1}}
+\begin{document}
+
+This program is devoted to two aspects: the use of mixed finite elements -- in
+particular Raviart-Thomas elements -- and using block matrices to define
+solvers and preconditioners that use the substructure of the system
+matrix. The equation we are going to solve is again the Laplace equation,
+though with a matrix-valued coefficient:
+\begin{align*}
+  -\nabla \cdot K(\vec x) \nabla p &= f \qquad && \text{in $\Omega$}, \\
+  p &= g && \text{on $\partial\Omega$}.
+\end{align*}
+$K(\vec x)$ is assumed to be uniformly positive definite, i.e. there is
+$\alpha>0$ such that the eigenvalues $\lambda_i(\vec x)$ of $K(x)$ satisfy
+$\lambda_i(\vec x)\ge \alpha$. The use of the symbol $p$ instead of the usual
+$u$ for the solution variable will become clear in the next section.
+
+After discussing the equation and the formulation we are going to use to solve
+it, this introduction will cover the use of block matrices and vectors, the
+definition of solvers and preconditioners, and finally the actual testcase we
+are going to solve.
+
+\subsection*{Formulation, weak form, and discrete problem}
+
+In the form above, the Laplace equation is considered a good model equation
+for fluid flow in porous media. In particular, if flow is so slow that all
+dynamic effects such as the acceleration terms in the Navier-Stokes equation
+become irrelevant, and if the flow pattern is stationary, then the Laplace
+equation models the pressure that drives the flow reasonable well. Because the
+solution variable is a pressure, we here use the name $p$ instead of the
+name $u$ more commonly used for the solution of partial differential equations.
+
+Typical applications of this view of the Laplace equation are then modeling
+groundwater flow, or the flow of hydrocarbons in oil reservoirs. In these
+applications, $K$ is then the permeability tensor, i.e. a measure for how much
+resistence the soil or rock matrix asserts on the fluid flow. In the
+applications just named, a desirable feature is that the numerical scheme is
+locally conservative, i.e. that whatever flows into a cell also flows out of
+it (or the difference is equal to the integral over the source terms over each
+cell, if the sources are nonzero). However, as it turns out, the usual
+discretizations of the Laplace equation do not satisfy this property. On the
+other hand, one can achieve this by choosing a different formulation.
+
+To this end, one first introduces a second variable, called the flux, $\vec
+u=-K\nabla p$. By its definition, the flux is a vector in the negative
+direction of the pressure gradient, multiplied by the permeability tensor. If
+the permeability tensor is proportional to the unit matrix, this equation is
+easy to understand and intuitive: the higher the permeability, the higher the
+flux; and the flux is proportional to the gradient of the pressure, going from
+areas of high pressure to areas of low pressure.
+
+With this second variable, one then finds an alternative version of the
+Laplace equation, called the mixed formulation:
+\begin{align*}
+  K^{-1} \vec u - \nabla p &= 0 \qquad && \text{in $\Omega$}, \\
+  -\text{div}\ \vec u &= 0 \qquad && \text{in $\Omega$}, \\
+  p &= g \qquad && \text{on $\partial\Omega$}.
+\end{align*}
+
+The weak formulation of this problem is found by multiplying the two
+equations with test functions and integrating some terms by parts:
+\begin{align*}
+  A(\{\vec u,p\},\{\vec v,q\}) = F(\{\vec v,q\}),
+\end{align*}
+where
+\begin{align*}
+  A(\{\vec u,p\},\{\vec v,q\})
+  &=
+  (\vec v, K^{-1}\vec u)_\Omega - (\text{div}\ \vec v, p)_\Omega
+  - (q,\text{div}\ \vec u)_\Omega
+  \\
+  F(\{\vec v,q\}) &= -(g,\vec v\cdot \vec n)_{\partial\Omega} - (f,q)_\Omega.
+\end{align*}
+Here, $\vec n$ is the outward normal vector at the boundary. Note how in this
+formulation, Dirichlet boundary values of the original problem are
+incorporated in the weak form.
+
+To be well-posed, we have to look for solutions and test functions in the
+space $H(\text{div})=\{\vec w\in L^2(\Omega)^d:\ \text{div}\ \vec w\in L^2\}$
+for $\vec u,\vec v$, and $L^2$ for $p,q$. It is a well-known fact stated in
+almost every book on finite element theory that if one chooses discrete finite
+element spaces for the approximation of $\vec u,p$ inappropriately, then the
+resulting discrete saddle-point problem is instable and the discrete solution
+will not converge to the exact solution.
+
+To overcome this, a number of different finite element pairs for $\vec u,p$
+have been developed that lead to a stable discrete problem. One such pair is
+to use the Raviart-Thomas spaces $RT(k)$ for the velocity $\vec u$ and
+discontinuous elements of class $DQ(k)$ for the pressure $p$. For details
+about these spaces, we refer in particular to the book on mixed finite element
+methods by Brezzi and Fortin, but many other books on the theory of finite
+elements, for example the classic book by Brenner and Scott, also state the
+relevant results.
+
+
+\subsection*{Assembling the linear system}
+
+The deal.II library (of course) implements Raviart-Thomas elements $RT(k)$ of
+arbitrary order $k$, as well as discontinuous elements $DG(k)$. If we forget
+about their particular properties for a second, we then have to solve a
+discrete problem
+\begin{align*}
+  A(x_h,w_h) = F(w_h),
+\end{align*}
+with the bilinear form and right hand side as stated above, and $x_h=\{\vec
+u_h,p_h\}$, $w_h=\{\vec v_h,q_h\}$. Both $x_h$ and $w_h$ are from the space
+$X_h=RT(k)\times DQ(k)$, where $RT(k)$ is itself a space of $dim$-dimensional
+functions to accomodate for the fact that the flow velocity is vector-valued.
+The necessary question then is: how do we do this in a program?
+
+Vector-valued elements have already been discussed in previous tutorial
+programs, the first time and in detail in step-8. The main difference there
+was that the vector-valued space $V_h$ is uniform in all its components: the
+$dim$ components of the displacement vector are all equal and from the same
+function space. What we could therefore do was to build $V_h$ as the outer
+product of the $dim$ times the usual $Q(1)$ finite element space, and by this
+make sure that all our shape functions have only a single non-zero vector
+component. Instead of dealing with vector-valued shape functions, all we did
+in step-8 was therefore to look at the (scalar) only non-zero component and
+use the \texttt{fe.system\_to\_component\_index(i).first} call to figure out
+which component this actually is.
+
+This doesn't work with Raviart-Thomas elements: following from their
+construction to satisfy certain regularity properties of the space
+$H(\text{div})$, the shape functions of $RT(k)$ are usually nonzero in all
+their vector components at once. For this reason, were
+\texttt{fe.system\_to\_component\_index(i).first} applied to determine the only
+nonzero component of shape function $i$, an exception would be generated. What
+we really need to do is to get at \textit{all} vector components of a shape
+function. In deal.II diction, we call such finite elements
+\textit{non-primitive}, whereas finite elements that are either scalar or for
+which every vector-valued shape function is nonzero only in a single vector
+component are called \textit{primitive}.
+
+So what do we have to do for non-primitive elements? To figure this out, let
+us go back in the tutorial programs, almost to the very beginnings. There, we
+learned that we use the \texttt{FEValues} class to determine the values and
+gradients of shape functions at quadrature points. For example, we would call
+\texttt{fe\_values.shape\_value(i,q\_point)} to obtain the value of the
+\texttt{i}th shape function on the quadrature point with number
+\texttt{q\_point}. Later, in step-8 and other tutorial programs, we learned
+that this function call also works for vector-valued shape functions (of
+primitive finite elements), and that it returned the value of the only
+non-zero component of shape function \texttt{i} at quadrature point
+\texttt{q\_point}.
+
+For non-primitive shape functions, this is clearly not going to work: there is
+no single non-zero vector component of shape function \texttt{i}, and the call 
+to \texttt{fe\_values.shape\_value(i,q\_point)} would consequently not make
+much sense. However, deal.II offers a second function call,
+\texttt{fe\_values.shape\_value\_component(i,q\_point,comp)} that returns the
+value of the \texttt{comp}th vector component of shape function  \texttt{i} at
+quadrature point \texttt{q\_point}, where \texttt{comp} is an index between
+zero and the number of vector components of the present finite element; for
+example, the element we will use to describe velocities and pressures is going
+to have $dim+1$ components. It is worth noting that this function call can
+also be used for primitive shape functions: it will simply return zero for all
+components except one; for non-primitive shape functions, it will in general
+return a non-zero value for more than just one component.
+
+We could now attempt to rewrite the bilinear form above in terms of vector
+components. For example, in 2d, the first term could be rewritten like this
+(note that $u_0=x_0, u_1=x_1, p=x_2$):
+\begin{align*}
+  (\vec u_h^i, K^{-1}\vec u_h^j)
+  =
+  &\left((x_h^i)_0, K^{-1}_{00} (x_h^j)_0\right) +
+   \left((x_h^i)_0, K^{-1}_{01} (x_h^j)_1\right) + \\
+  &\left((x_h^i)_1, K^{-1}_{10} (x_h^j)_0\right) +
+   \left((x_h^i)_1, K^{-1}_{11} (x_h^j)_1\right).
+\end{align*}
+If we implemented this, we would get code like this:
+\begin{verbatim}
+  for (unsigned int q=0; q<n_q_points; ++q) 
+    for (unsigned int i=0; i<dofs_per_cell; ++i)
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+        local_matrix(i,j) += (Kinverse[q]][0][0] *
+                              fe_values.shape_value_component(i,q,0) *
+                              fe_values.shape_value_component(j,q,0) 
+                              +
+                              Kinverse[q]][0][1] *
+                              fe_values.shape_value_component(i,q,0) *
+                              fe_values.shape_value_component(j,q,1) 
+                              +
+                              Kinverse[q]][1][0] *
+                              fe_values.shape_value_component(i,q,1) *
+                              fe_values.shape_value_component(j,q,0) 
+                              +
+                              Kinverse[q]][1][1] *
+                              fe_values.shape_value_component(i,q,1) *
+                              fe_values.shape_value_component(j,q,1) 
+                             )
+                             *
+                             fe_values.JxW(q);
+\end{verbatim}
+This is, at best, tedious, error prone, and not dimension independent. There
+are obvious ways to make things dimension independent, but in the end, the
+code is simply not pretty. What would be much nicer is if we could simply
+extract the $\vec u$ and $p$ components of a shape function $x_h^i$. In the
+program we do that, by writing functions like this one:
+\begin{verbatim}
+template <int dim>
+Tensor<1,dim>
+extract_u (const FEValuesBase<dim> &fe_values,
+           const unsigned int i,
+           const unsigned int q)
+{
+  Tensor<1,dim> tmp;
+
+  for (unsigned int d=0; d<dim; ++d)
+    tmp[d] = fe_values.shape_value_component (i,q,d);
+
+  return tmp;
+}
+\end{verbatim}
+
+What this function does is, given an \texttt{fe\_values} object, to extract
+the values of the first $dim$ components of shape function \texttt{i} at
+quadrature points \texttt{q}, that is the velocity components of that shape
+function. Put differently, if we write shape functions $x_h^i$ as the tuple
+$\{\vec u_h^i,p_h^i\}$, then the function returns the velocity part of this
+tuple. Note that the velocity is of course a $dim$-dimensional tensor, and
+that the function returns a corresponding object.
+
+Likewise, we have a function that extracts the pressure component of a shape
+function:
+\begin{verbatim}
+template <int dim>
+double extract_p (const FEValuesBase<dim> &fe_values,
+                  const unsigned int i,
+                  const unsigned int q)
+{
+  return fe_values.shape_value_component (i,q,dim);
+}
+\end{verbatim}
+Finally, the bilinear form contains terms involving the gradients of the
+velocity component of shape functions. To be more precise, we are not really
+interested in the full gradient, but only the divergence of the velocity
+components, i.e. $\text{div}\ \vec u_h^i=\sum_{d=0}^{dim-1} \frac \partial
+{\partial x_d} (\vec
+u_h^i)_d$. Here's a function that returns this quantity:
+\begin{verbatim}
+template <int dim>
+double
+extract_div_u (const FEValuesBase<dim> &fe_values,
+               const unsigned int i,
+               const unsigned int q)
+{
+  double divergence = 0;
+  for (unsigned int d=0; d<dim; ++d)
+    divergence += fe_values.shape_grad_component (i,q,d)[d];
+
+  return divergence;
+}
+\end{verbatim}
+
+With these three functions, all of which are completely dimension independent
+and will therefore also work in 3d, assembling the local matrix and right hand
+side contributions becomes a charm:
+\begin{verbatim}
+for (unsigned int q=0; q<n_q_points; ++q) 
+  for (unsigned int i=0; i<dofs_per_cell; ++i)
+    {
+      const Tensor<1,dim> phi_i_u = extract_u (fe_values, i, q);
+      const double div_phi_i_u    = extract_div_u (fe_values, i, q);
+      const double phi_i_p        = extract_p (fe_values, i, q);
+           
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+        {
+          const Tensor<1,dim> phi_j_u = extract_u (fe_values, j, q);
+          const double div_phi_j_u    = extract_div_u (fe_values, j, q);
+          const double phi_j_p        = extract_p (fe_values, j, q);
+               
+          local_matrix(i,j) += (phi_i_u * Kinverse[q] * phi_j_u
+                                - div_phi_i_u * phi_j_p
+                                - phi_i_p * div_phi_j_u)
+                               * fe_values.JxW(q);
+        }
+
+      local_rhs(i) += -(phi_i_p *
+                        rhs_values[q] *
+                        fe_values.JxW(q));
+    }
+\end{verbatim}
+This very closely resembles the form we have originally written down the
+bilinear form and right hand side.
+
+There is one final term that we have to take care of: the right hand side
+contained the term $(g,\vec v\cdot \vec n)_{\partial\Omega}$, constituting the
+weak enforcement of pressure boundary conditions. We have already seen in
+step-7 how to deal with face integrals: essentially exactly the same as with
+domain integrals, except that we have to use the \texttt{FEFaceValues} class
+instead of \texttt{FEValues}. To compute the boundary term we then simply have
+to loop over all boundary faces and integrate there. If you look closely at
+the definitions of the \texttt{extract\_*} functions above, you will realize
+that it isn't even necessary to write new functions that extract the velocity
+and pressure components of shape functions from \texttt{FEFaceValues} objects:
+both \texttt{FEValues} and \texttt{FEFaceValues} are derived from a common
+base class, \texttt{FEValuesBase}, and the extraction functions above can
+therefore deal with both in exactly the same way. Assembling the missing
+boundary term then takes on the following form:
+\begin{verbatim}
+for (unsigned int face_no=0;
+     face_no<GeometryInfo<dim>::faces_per_cell;
+     ++face_no)
+  if (cell->at_boundary(face_no))
+    {
+      fe_face_values.reinit (cell, face_no);
+    
+      pressure_boundary_values
+        .value_list (fe_face_values.get_quadrature_points(),
+                     boundary_values);
+
+      for (unsigned int q=0; q<n_face_q_points; ++q) 
+        for (unsigned int i=0; i<dofs_per_cell; ++i)
+          {
+            const Tensor<1,dim>
+              phi_i_u = extract_u (fe_face_values, i, q);
+                
+            local_rhs(i) += -(phi_i_u *
+                              fe_face_values.normal_vector(q) *
+                              boundary_values[q] *
+                              fe_face_values.JxW(q));
+        }
+  }
+\end{verbatim}
+
+You will find the exact same code as above in the sources for the present
+program. We will therefore not comment much on it below.
+
+
+\subsection*{Linear solvers and preconditioners}
+
+
+
+\subsection*{Definition of the testcase}
+
+In this tutorial program, we will solve the Laplace equation in mixed
+formulation as stated above. Since we want to monitor convergence of the
+solution inside the program, we choose right hand side, boundary conditions,
+and the coefficient so that we recover a solution function known to us. In
+particular, we choose the pressure solution
+\begin{align*}
+  p = -\left(\frac \alpha 2 xy^2 + \beta x - \frac \alpha 6 x^2\right),
+\end{align*}
+and for the coefficient we choose the unit matrix $K_{ij}=\delta_{ij}$ for
+simplicity. Consequently, the exact velocity satisfies
+\begin{align*}
+  \vec u = 
+  \begin{pmatrix}
+    \frac \alpha 2 y^2 + \beta - \frac \alpha 2 x^2 \\
+    \alpha xy
+  \end{pmatrix}.
+\end{align*}
+This solution was chosen since it is exactly divergence free, making it a
+realistic testcase for incompressible fluid flow. By consequence, the right
+hand side equals $f=0$, and as boundary values we have to choose
+$g=p|_{\partial\Omega}$.
+
+For the computations in this program, we choose $\alpha=0.3,\beta=1$. You can
+find the resulting solution in the ``Results'' section below, after the
+commented program.
+
+\end{document}

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.