-
Notifications
You must be signed in to change notification settings - Fork 5
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
1 changed file
with
25 additions
and
26 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,8 +1,8 @@ | ||
Boffi-Lovadina Augmentation for Bernardi-Raugel discretization of the Stokes problem | ||
Boffi-Lovadina Augmentation for Bernardi-Raugel discretisation of the Stokes problem | ||
====================================================================================== | ||
|
||
In this tutorial, we explore constructing preconditioners for Bernardi-Raugel discretizations of the Stokes problem with Boffi-Lovadina augmentation. | ||
In particular, we will consider a Bernardi-Raugel inf-sup stable discretization of the Stokes problem, i.e. | ||
In this tutorial, we explore constructing preconditioners for Bernardi-Raugel discretisations of the Stokes problem with Boffi-Lovadina augmentation. | ||
In particular, we will consider a Bernardi-Raugel inf-sup stable discretisation of the Stokes problem, i.e. | ||
|
||
.. math:: | ||
|
@@ -13,7 +13,7 @@ In particular, we will consider a Bernardi-Raugel inf-sup stable discretization | |
(\nabla\cdot \vec{u},q)_{L^2(\Omega)} = 0 \qquad \forall q\in L^2(\Omega) | ||
\end{cases} | ||
Such a discretization can easily be constructed using NGSolve as follows: :: | ||
Such a discretisation can easily be constructed using NGSolve as follows: :: | ||
|
||
from ngsolve import * | ||
from netgen.occ import * | ||
|
@@ -46,8 +46,8 @@ Such a discretization can easily be constructed using NGSolve as follows: :: | |
f = LinearForm(V).Assemble() | ||
g = LinearForm(Q).Assemble(); | ||
|
||
We can now explore what happens if we use a Schour complement preconditioner for the saddle point problem. | ||
We can construct the Schur complement preconditioner using the following code: :: | ||
We can now explore what happens if we use a Schur complement preconditioner for the saddle point problem. | ||
We can construct the (dense!) Schur complement preconditioner using the following code: :: | ||
|
||
K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, None] ] ) | ||
from ngsPETSc.pc import * | ||
|
@@ -72,9 +72,9 @@ We can construct the Schur complement preconditioner using the following code: : | |
printrates=True, initialize=False) | ||
Draw(gfu) | ||
|
||
The Schur complement preconditioner converges in a few iterations, but it is not very efficient since we need to invert the Schur complement. | ||
The Schur complement preconditioner converges in a few iterations, but as written it is not very efficient since we need to assemble and invert the dense Schur complement. | ||
|
||
.. list-table:: Preconditioners performance | ||
.. list-table:: Preconditioner performance | ||
:widths: auto | ||
:header-rows: 1 | ||
|
||
|
@@ -83,10 +83,9 @@ The Schur complement preconditioner converges in a few iterations, but it is not | |
* - Schur complement | ||
- 4 (9.15e-14) | ||
|
||
In fact, the Schur complement is dense hence inverting it is not a good idea. Furthermore, to perform the inversion of the Schur complement we had to write a lot of "boilerplate" code. | ||
Since our discretization is inf-sup stable it is possible to prove that the mass matrix of the pressure space is spectrally equivalent to the Schur complement. | ||
Since our discretisation is inf-sup stable it is possible to prove that the mass matrix of the pressure space is spectrally equivalent to the Schur complement. | ||
This means that we can use the mass matrix of the pressure space as a preconditioner for the Schur complement. | ||
Notice that we still need to invert the mass matrix and we will do so using a `PETSc PC` of type Jacobi, which is the exact inverse since we are using `P0` elements. | ||
Notice that we still need to invert the mass matrix and we will do so using a `PETSc PC` of type Jacobi, which is the exact inverse since we are using `P0` pressure elements. | ||
We will also invert the Laplacian block using a `PETSc PC` of type `LU`. :: | ||
|
||
m = BilinearForm((1/nu)*p*q*dx).Assemble() | ||
|
@@ -103,7 +102,7 @@ We will also invert the Laplacian block using a `PETSc PC` of type `LU`. :: | |
maxsteps=100, printrates=True, initialize=False) | ||
Draw(gfu) | ||
|
||
.. list-table:: Preconditioners performance | ||
.. list-table:: Preconditioner performance | ||
:widths: auto | ||
:header-rows: 1 | ||
|
||
|
@@ -145,7 +144,7 @@ We can also construct a multi-grid preconditioner for the top left block of the | |
solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-8, | ||
maxsteps=100, printrates=True, initialize=False) | ||
.. list-table:: Preconditioners performance | ||
.. list-table:: Preconditioner performance | ||
:widths: auto | ||
:header-rows: 1 | ||
|
||
|
@@ -158,7 +157,7 @@ We can also construct a multi-grid preconditioner for the top left block of the | |
* - Mass & Two Level Additive Schwarz | ||
- 100 (4.68e-06) | ||
|
||
The preconditioner constructed so far doesn't seem to be ideal, in fact, our Krylov solver took many iterations to converge with a direct LU factorization of the velocity block and did not converge at all with `HYPRE`. | ||
The preconditioner constructed so far doesn't seem to be ideal, in fact, our Krylov solver took many iterations to converge with a direct LU factorisation of the velocity block and did not converge at all with `HYPRE`. | ||
To resolve this issue we resort to an augmented Lagrangian formulation, i.e. | ||
|
||
.. math:: | ||
|
@@ -191,7 +190,7 @@ This formulation can easily be constructed in NGSolve, as follows: :: | |
Using an augmented Lagrangian formulation and a direct inversion of the velocity block, we were able to converge in only two iterations. | ||
This is because the augmented Lagrangian formulation improves the spectral equivalence between the mass matrix of the pressure space and the Schur complement. | ||
|
||
.. list-table:: Preconditioners performance | ||
.. list-table:: Preconditioner performance | ||
:widths: auto | ||
:header-rows: 1 | ||
|
||
|
@@ -206,8 +205,8 @@ This is because the augmented Lagrangian formulation improves the spectral equiv | |
* - Augmented Lagrangian LU | ||
- 2 (9.24e-8) | ||
|
||
Notice that so far we have been inverting the matrix corresponding to the Laplacian block using a direct LU factorization. | ||
This is not ideal for large problems, and we can try to use a `Hypre` preconditioner for the Laplacian block. :: | ||
Notice that so far we have been inverting the matrix corresponding to the Laplacian block using a direct LU factorisation. | ||
This is not ideal for large problems, and we can try to use a `HYPRE` preconditioner for the Laplacian block. :: | ||
|
||
smoother = aG.mat.CreateBlockSmoother(blocks) | ||
preHG = PETScPreconditioner(aG.mat, vertexdofs, solverParameters={"pc_type":"hypre"}) | ||
|
@@ -221,7 +220,7 @@ This is not ideal for large problems, and we can try to use a `Hypre` preconditi | |
Draw(gfu) | ||
|
||
|
||
.. list-table:: Preconditioners performance | ||
.. list-table:: Preconditioner performance | ||
:widths: auto | ||
:header-rows: 1 | ||
|
||
|
@@ -241,7 +240,7 @@ This is not ideal for large problems, and we can try to use a `Hypre` preconditi | |
|
||
Our first attempt at using a `HYPRE` preconditioner for the Laplacian block did not converge. | ||
This is because the top left block of the saddle point problem now contains the augmentation term, which has a very large kernel. | ||
It is well known that algebraic multi-grid methods do not work well with indefinite problems, and this is what we are observing here. | ||
It is well known that algebraic multi-grid methods do not work well for such problems, and this is what we are observing here. | ||
We begin by constructing the augmented Lagrangian formulation in more numerical linear algebra terms, i.e. | ||
|
||
.. math:: | ||
|
@@ -284,12 +283,12 @@ We can construct this linear algebra problem inside NGSolve as follows: :: | |
printrates=True, initialize=False) | ||
Draw(gfu) | ||
|
||
Since the augmentation block has a lower rank than the Laplacian block, we can use the Sherman-Morrisson-Woodbory formula to invert the augmentation block. | ||
Since the augmentation block has a lower rank than the Laplacian block, we can use the Sherman-Morrison-Woodbury formula to invert the augmentation block. | ||
|
||
.. math:: | ||
(A + B^T(\gamma M^{-1})B)^{-1} = A^{-1} - A^{-1}B^T(\frac{1}{\gamma}M^{-1} + BA^{-1}B^T)^{-1}BA^{-1} | ||
We will do this in two different ways first we will invert the :math:`(\frac{1}{\gamma}M^{-1} + BA^{-1}B^T)` block using a direct LU factorization. | ||
We will do this in two different ways first we will invert the :math:`(\frac{1}{\gamma}M^{-1} + BA^{-1}B^T)` block using a direct LU factorisation. | ||
Then we will notice that since the penalisation parameter is large we can ignore the :math:`\frac{1}{\gamma}M^{-1}` term and use the mass matrix since it is spectrally equivalent to the Schur complement. :: | ||
|
||
SM = (d.mat + b.mat@[email protected]).ToDense().NumPy() | ||
|
@@ -308,7 +307,7 @@ Then we will notice that since the penalisation parameter is large we can ignore | |
gfu.Set(uin, definedon=mesh.Boundaries("inlet")) | ||
sol = BlockVector( [gfu.vec, gfp.vec] ) | ||
|
||
print("-----------|Boffi--Lovadina Augmentation Sherman-Morrisson-Woodbory|-----------") | ||
print("-----------|Boffi--Lovadina Augmentation Sherman-Morrison-Woodbury|-----------") | ||
solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-10, | ||
printrates=True, initialize=False) | ||
Draw(gfu) | ||
|
@@ -319,14 +318,14 @@ Then we will notice that since the penalisation parameter is large we can ignore | |
gfu.Set(uin, definedon=mesh.Boundaries("inlet")) | ||
sol = BlockVector( [gfu.vec, gfp.vec] ) | ||
|
||
print("-----------|Boffi--Lovadina Augmentation Sherman-Morrisson-Woodbory|-----------") | ||
print("-----------|Boffi--Lovadina Augmentation Sherman-Morrison-Woodbury|-----------") | ||
solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, tol=1e-13, | ||
printrates=True, initialize=False) | ||
Draw(gfu) | ||
|
||
We see that a purely algebraic approach based on the Sherman-Morrisson-Woodbory formula is more efficient for the augmented Lagrangian formulation than a naive two-level additive Schwarz approach. | ||
We see that a purely algebraic approach based on the Sherman-Morrison-Woodbury formula is more efficient for the augmented Lagrangian formulation than a naive two-level additive Schwarz approach. | ||
|
||
.. list-table:: Preconditioners performance | ||
.. list-table:: Preconditioner performance | ||
:widths: auto | ||
:header-rows: 1 | ||
|
||
|
@@ -342,5 +341,5 @@ We see that a purely algebraic approach based on the Sherman-Morrisson-Woodbory | |
- 2 (9.24e-08) | ||
* - Augmented Two Level Additive Schwarz | ||
- 100 (1.06e-03) | ||
* - Augmentation Schermon-Morrisson-Woodbory | ||
* - Augmentation Schermon-Morrison-Woodbury | ||
- 84 (1.16e-07) |