From f1790ee2013ee85a4a052f1daf08da4fef775040 Mon Sep 17 00:00:00 2001 From: wolf Date: Sat, 11 Feb 2006 00:26:02 +0000 Subject: [PATCH] First version of intro. git-svn-id: https://svn.dealii.org/trunk@12310 0785d39b-7218-0410-832d-ea1e28bc413d --- .../step-20.data/intro.tex | 366 ++++++++++++++++++ 1 file changed, 366 insertions(+) create mode 100644 deal.II/doc/tutorial/chapter-2.step-by-step/step-20.data/intro.tex 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 index 0000000000..2d5851b334 --- /dev/null +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-20.data/intro.tex @@ -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 +Tensor<1,dim> +extract_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + Tensor<1,dim> tmp; + + for (unsigned int d=0; d +double extract_p (const FEValuesBase &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 +double +extract_div_u (const FEValuesBase &fe_values, + const unsigned int i, + const unsigned int q) +{ + double divergence = 0; + for (unsigned int d=0; d 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 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::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 + 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} -- 2.39.5