all places.
+<h4> A note on the structure of the linear system </h4>
+
+Above, we have claimed that the linear system has the form
+@f{eqnarray*}
+ \left(\begin{array}{cc}
+ A & B^T \\ B & 0
+ \end{array}\right)
+ \left(\begin{array}{cc}
+ U \\ P
+ \end{array}\right)
+ =
+ \left(\begin{array}{cc}
+ F \\ G
+ \end{array}\right),
+@f}
+i.e., in particular that there is a zero block at the bottom right of the
+matrix. This then allowed us to write the Schur complement as
+$S=B A^{-1} B^T$. But this is not quite correct.
+
+Think of what would happen if there are constraints on some
+pressure variables (see the
+@ref constraints "Constraints on degrees of freedom" documentation
+module), for example because we use adaptively
+refined meshes and continuous pressure finite elements so that there
+are hanging nodes. Another cause for such constraints are Dirichlet
+boundary conditions on the pressure. Then the AffineConstraints
+class, upon copying the local contributions to the matrix into the
+global linear system will zero out rows and columns corresponding
+to constrained degrees of freedom and put a positive entry on
+the diagonal. (You can think of this entry as being one for
+simplicity, though in reality it is a value of the same order
+of magnitude as the other matrix entries.) In other words,
+the bottom right block is really not empty at all: It has
+a few entries on the diagonal, one for each constrained
+pressure degree of freedom, and a correct description
+of the linear system we have to solve is that it has the
+form
+@f{eqnarray*}
+ \left(\begin{array}{cc}
+ A & B^T \\ B & D_c
+ \end{array}\right)
+ \left(\begin{array}{cc}
+ U \\ P
+ \end{array}\right)
+ =
+ \left(\begin{array}{cc}
+ F \\ G
+ \end{array}\right),
+@f}
+where $D_c$ is the zero matrix with the exception of the
+positive diagonal entries for the constrained degrees of
+freedom. The correct Schur complement would then in fact
+be the matrix $S = B A^{-1} B^T - D_c $ instead of the one
+stated above.
+
+Thinking about this makes us, first, realize that the
+resulting Schur complement is now indefinite because
+$B A^{-1} B^T$ is symmetric and positive definite whereas
+$D_c$ is a positive semidefinite, and subtracting the latter
+from the former may no longer be positive definite. This
+is annoying because we could no longer employ the Conjugate
+Gradient method on this true Schur complement. That said, we could
+fix the issue in AffineConstraints::distribute_local_to_global() by
+simply putting *negative* values onto the diagonal for the constrained
+pressure variables -- because we really only put something nonzero
+to ensure that the resulting matrix is not singular; we really didn't
+care whether that entry is positive or negative. So if the entries
+on the diagonal of $D_c$ were negative, then $S$ would again be a
+symmetric and positive definite matrix.
+
+But, secondly, the code below doesn't actually do any of that: It
+happily solves the linear system with the wrong Schur complement
+$S = B A^{-1} B^T$ that just ignores the issue altogether. Why
+does this even work? To understand why this is so, recall that
+when writing local contributions into the global matrix,
+AffineConstraints::distribute_local_to_global() zeros out the
+rows and columns that correspond to constrained degrees of freedom.
+This means that $B$ has some zero rows, and $B^T$ zero columns.
+As a consequence, if one were to multiply out what the entries
+of $S$ are, one would realize that it has zero rows and columns
+for all constrained pressure degrees of freedom, including a
+zero on the diagonal. The nonzero entries of $D_c$ would fit
+into exactly those zero diagonal locations, and ensure that $S$
+is invertible. Not doing so, strictly speaking, means that $S$
+remains singular: It is symmetric and positive definite on the
+subset of non-constrained pressure degrees of freedom, and
+simply the zero matrix on the constrained pressures. Why
+does the Conjugate Gradient method work for this matrix?
+Because AffineConstraints::distribute_local_to_global()
+also makes sure that the right hand side entries that
+correspond to these zero rows of the matrix are *also*
+zero, i.e., the right hand side is compatible.
+
+What this means is that whatever the values of the solution
+vector for these constrained pressure degrees of freedom,
+these rows will always have a zero residual and, if one
+were to consider what the CG algorithm does internally, just
+never produce any updates to the solution vector. In other
+words, the CG algorithm just *ignores* these rows, despite the
+fact that the matrix is singular. This only works because these
+degrees of freedom are entirely decoupled from the rest of the
+linear system (because the entire row and corresponding column
+are zero). At the end of the solution process, the constrained
+pressure values in the solution vector therefore remain exactly
+as they were when we started the call to the solver; they are
+finally overwritten with their correct values when we call
+AffineConstraints::distribute() after the CG solver is done.
+
+The upshot of this discussion is that the assumption that the
+bottom right block of the big matrix is zero is a bit
+simplified, but that just going with it does not actually lead
+to any practical problems worth addressing.
+
+
<h3>The testcase</h3>
The domain, right hand side and boundary conditions we implement below relate
projects / dealii.git / commitdiff