From: Wolfgang Bangerth Date: Thu, 7 Nov 2019 16:34:03 +0000 (-0700) Subject: Add a discussion about Schur complements to step-22. X-Git-Tag: v9.2.0-rc1~881^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F9018%2Fhead;p=dealii.git Add a discussion about Schur complements to step-22. --- diff --git a/examples/step-22/doc/intro.dox b/examples/step-22/doc/intro.dox index 74a11a1dde..06c0c2beec 100644 --- a/examples/step-22/doc/intro.dox +++ b/examples/step-22/doc/intro.dox @@ -571,6 +571,120 @@ classes, we will be able to use them interchangeably using the same syntax in all places. +

A note on the structure of the linear system

+ +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. + +

The testcase

The domain, right hand side and boundary conditions we implement below relate