FENaPack - FEniCS Navier-Stokes preconditioning package

https://circleci.com/gh/blechta/fenapack.svg?style=svg

FENaPack is a package implementing preconditioners for Navier-Stokes problem using FEniCS and PETSc packages. In particular, variants of PCD (pressure-convection-diffussion) preconditioner from [1], [2] are implemented.

[1](1, 2, 3) Elman H. C., Silvester D. J., Wathen A. J., Finite Elements and Fast Iterative Solvers: With Application in Incompressible Fluid Dynamics. Oxford University Press 2005. 2nd edition 2014.
[2](1, 2) Olshanskii M. A., Vassilevski Y. V., Pressure Schur complement preconditioners for the discrete Oseen problem. SIAM J. Sci. Comput., 29(6), 2686-2704. 2007.

Usage

To use FENaPack matching version of FEniCS (version 2017.2.0) compiled with PETSc and petsc4py is needed. To install FENaPack from source do:

pip install [--user|--prefix=...] [-e] .

in the source/repository root dir. Editable install using -e allows to use FENaPack directly from source directory while editing it which is suitable for development.

Meshes for running demos can be downloaded from the FEniCS project website by executing download-meshes script. Demos can be run by navigating to a particular demo directory and typing:

NP=4
mpirun -n $NP python demo_foo-bar.py [-h]

Full documentation is available at https://fenapack.readthedocs.io/.

Authors

License

FENaPack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

FENaPack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License along with FENaPack. If not, see <http://www.gnu.org/licenses/>.

Acknowledgement

This work was supported by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project „IT4Innovations National Supercomputing Center – LM2015070“.

Mathematical background

Navier-Stokes equations

\[\begin{split}\left[\begin{array}{cc} -\nu\Delta + \mathbf{v}\cdot\nabla & \nabla \\ -\operatorname{div} & 0 \end{array}\right] \left[\begin{array}{c} \mathbf{u} \\ p \end{array}\right] = \left[\begin{array}{c} \mathbf{f} \\ 0 \end{array}\right]\end{split}\]

solved by GMRES preconditioned from right by

\[\begin{split}\mathbb{P} := \left[\begin{array}{cc} -\nu\Delta + \mathbf{v}\cdot\nabla & \nabla \\ & -\mathbb{S} \end{array}\right]\end{split}\]

with Schur complement \(\mathbb{S} = -\operatorname{div}\left(-\nu\Delta+\mathbf{v}\cdot\nabla\right)^{-1}\nabla\) would converge in two iterations. Unfortunately \(\mathbb{S}\) is dense. Possible trick is to approximate \(\mathbb{S}\) by swapping the order of the operators

\[\mathbb{S} \approx \mathbb{X}_\mathrm{BRM1} := (-\Delta) \left(-\nu\Delta+\mathbf{v}\cdot\nabla\right)^{-1}\]

or

\[\mathbb{S} \approx \mathbb{X}_\mathrm{BRM2} := \left(-\nu\Delta+\mathbf{v}\cdot\nabla\right)^{-1}(-\Delta).\]

This gives rise to the action of 11-block of preconditioner \(\mathbb{P}^{-1}\) given by

\[\mathbb{X}_\mathrm{BRM1}^{-1} := \left(-\nu\Delta+\mathbf{v}\cdot\nabla\right)(-\Delta)^{-1}.\]

or

\[\mathbb{X}_\mathrm{BRM2}^{-1} := (-\Delta)^{-1}\left(-\nu\Delta+\mathbf{v}\cdot\nabla\right).\]

Obviously additional artificial boundary condition for Laplacian solve \(-\Delta^{-1}\) is needed in the action of preconditioner. Modifying the approach from [2] we implement \(\mathbb{X}_\mathrm{BRM1}^{-1}\) as

\[\mathbb{X}_\mathrm{BRM1}^{-1} := \mathbb{M}_p^{-1} (\mathbb{I} + \mathbb{K}_p\mathbb{A}_p^{-1})\]

where \(\mathbb{M}_p\) is \(\nu^{-1}\)-multiple of mass matrix on pressure, \(\mathbb{K}_p \approx \nu^{-1}\mathbf{v}\cdot\nabla\) is a pressure convection matrix, and \(\mathbb{A}_p^{-1} \approx (-\Delta)^{-1}\) is a pressure Laplacian solve with zero boundary condition on inlet. This is implemented by fenapack.preconditioners.PCDPC_BRM1 and PCD preconditioner for Navier-Stokes equations.

Analogically we prefer to express BRM2 approach as

\[\mathbb{X}_\mathrm{BRM2}^{-1} := (\mathbb{I} + \mathbb{A}_p^{-1}\mathbb{K}_p) \mathbb{M}_p^{-1}\]

now with zero boundary condition on outlet for Laplacian solve and additional Robin term in convection matrix \(\mathbb{K}_p\) roughly as stated in [1], section 9.2.2. See also PCD preconditioner for Navier-Stokes equations and fenapack.preconditioners.PCDPC_BRM2.

Extension to time-dependent problems (PCDR preconditioners)

Time disretization applied in unsteady problems typically leads to the need to incorporate a reaction term into the preconditioner. Typically, we end up with

\[\mathbb{X}_\mathrm{BRM1}^{-1} := \left(\frac{1}{\tau}-\nu\Delta+\mathbf{v}\cdot\nabla\right)(-\Delta)^{-1},\]

or

\[\mathbb{X}_\mathrm{BRM2}^{-1} := (-\Delta)^{-1}\left(\frac{1}{\tau}-\nu\Delta+\mathbf{v}\cdot\nabla\right),\]

where \(\tau\) denotes a fixed time step and the original PCD preconditioner thus becomes PCDR (pressure-convection-diffusion-reaction) preconditioner. A straightforward way of how to implement the above actions is to update the pressure convection matrix \(\mathbb{K}_p\) by a contribution corresponding to the scaled pressure mass matrix, namely

\[\mathbb{X}_\mathrm{BRM1}^{-1} := \mathbb{M}_p^{-1} \left(\mathbb{I} + \left(\mathbb{K}_p + \tau^{-1} \mathbb{M}_p\right)\mathbb{A}_p^{-1}\right),\]

or

\[\mathbb{X}_\mathrm{BRM2}^{-1} := \left(\mathbb{I} + \mathbb{A}_p^{-1}\left(\mathbb{K}_p + \tau^{-1} \mathbb{M}_p\right)\right)\mathbb{M}_p^{-1}.\]

However, for unsteady problems we prefer to use the following elaborated implementation of PCDR preconditioners, namely

\[\mathbb{X}_\mathrm{BRM1}^{-1} := \mathbb{R}_p^{-1} + \mathbb{M}_p^{-1} (\mathbb{I} + \mathbb{K}_p\mathbb{A}_p^{-1}),\]

or

\[\mathbb{X}_\mathrm{BRM2}^{-1} := \mathbb{R}_p^{-1} + (\mathbb{I} + \mathbb{A}_p^{-1}\mathbb{K}_p) \mathbb{M}_p^{-1},\]

where \(\mathbb{R}_p^{-1} \approx \frac{1}{\tau} (-\Delta)^{-1}\), while \(\mathbb{R}_p\) itself is approximated and implemented as

\[\mathbb{R}_p := \mathbb{B} \left(\tau^{-1} \mathbb{D}_\mathrm{M}\right)^{-1} \mathbb{B}^T,\]

Here, \(\mathbb{D}_\mathrm{M}\) is the diagonal of the velocity mass matrix, \(\mathbb{D}_\mathrm{M} = \operatorname{diag}(\mathbb{M}_{\mathbf{u}})\), and \(\mathbb{B}^T\) corresponds to the discrete pressure gradient which is obtained as the 01-block of the original system matrix. Let us emphasize that this submatrix is extracted from the system matrix with velocity Dirichlet boundary conditions being applied on it.

The choice of \(\mathbb{R}_p\) as above can be justified especially in the case of \(\tau \rightarrow 0_+\), for which

\[\mathbb{S}^{-1} := \left(-\operatorname{div}\left(\frac{1}{\tau} - \nu\Delta+\mathbf{v}\cdot\nabla\right)^{-1}\nabla\right)^{-1} \approx \frac{1}{\tau}\left(\mathbb{B} \mathbb{M}_{\mathbf{u}}^{-1} \mathbb{B}^T\right)^{-1},\]

and simultaneously \(\mathbb{X}^{-1} \approx \mathbb{R}_p^{-1} = \frac{1}{\tau} \left(\mathbb{B} \mathbb{D}_\mathrm{M}^{-1} \mathbb{B}^T\right)^{-1}\). The same approximation of the minus Laplacian operator was previously used also in [1], see Remark 9.6 therein.

Developer's resources

Manual and API Reference