Filter techniques

The linear solver library is designed as a black box library. Once linsol_solveAdaptively or linsol_precondDefect is called, the solver does not interact with the main program until a solution is computed or the solver broke down for any reason. Hereby, the basic solution process is based on a system matrix and a right-hand side vector.

Filters provide a flexible way to impose additional restrictions on the matrix which would otherwise be hard to be implemented. They are mainly designed for imposing some integral-mean-value condition into a vector, but their use is also rather advantageous for the implementation of boundary conditions in multigrid-based solvers. In particular, the implementation of Dirichlet boundary conditions is completely handled via vector filters.

What is a vector filter?
A vector filter is a numerical algorithm that forces a solution, right-hand side, defect or correction vector into a proper subspace of a finite element space. As a special property, vector filters can be specified in iterative solvers to modify the iteration.
What is a matrix filter?
A matrix filter is actually not a filter. It is a matrix modification routine that prepares a matrix to work with a specific vector filter. Matrix filters always belong to vector filters.

Involved modules

The following modules are involved in the use of filters:

Module Description
vectorfilters.f90 Implementation of vector filters
matrixfilters.f90 Implementation of matrix filters
filtersupport.f90 Support structures for filters being used in linear solvers
linearsolver.f90 The linear solver library uses filters during the solution process

Examples

In the following, we give a fundamental example that demonstrate the use of vector filters.

Example: The integral-mean-value filter

We consider the pure Neumann Poisson problem $$ \begin{align} -\Delta u &= f \qquad\text{in $\Omega$}\\ \partial_n u &= 0\qquad\text{on $\partial\Omega$} \end{align}$$ on a domain $\Omega\subset\mathbb{R}^d$ for a proper right­-hand side $f$ which fulfils the (necessary) compatibility condition $$\int_{\Omega}f=0.$$ The solution $u$ is to be found in $$ u \in L^2_0(\Omega) := \left\{ v\in L_2(\Omega)\,\bigg|\, \int_\Omega v = 0 \right\} \subset L_2(\Omega). $$

Let $V_h\subset H^1(\Omega)$ be a finite element subspace. The discrete counterpart to this problem is: Find $u_h\in V_h^0 := V_h\cap L^2_0$, such that $$ a(u_h,v_h) = (\nabla u_h, \nabla v_h) = (f, v_h) \qquad\forall\,v_h\in H^1(\Omega) $$ The corresponding discrete matrix-vector system takes the usual form $$ A \vec u = \vec f $$ for the system matrix $A$, the vector of degrees of freedom $\vec u$ corresponding to $u_h$ and the vector of right-hand side entries $\vec f$.

The additional restriction $u_h\in L^2_0$ drives the problem definite. However, the implementation of this condition into the corresponding system matrix is not straightforward as the condition $\int_{\Omega} u_h=0$ is an algebraic side condition through the complete domain. One method to circumvent this difficulty is the use of a so-called "integral mean value = zero" filter which is applied during an iteration.

A "defect vector filter" $F: \mathbb{R}^n \to \mathbb{R}^n$ in an iterative linear solver can be understood as an additional preconditioner. A very simple iterative solver is the preconditioned defect correction loop with a preconditioner $C:\mathbb{R}^n \to \mathbb{R}^n$:

$$\vec u_{n+1} = \vec u_n + C^{-1} F \underbrace{(\vec f - A \vec u_n)}_{:=\vec d_n} $$ The idea is as follows:

  • Let $\int_{\Omega} u_0=0$. For example, set $\vec u_0=0$.

  • Let $F d_n$ denote the continuous counterpart to $F\vec d_n$. The "defect filter" $F$ modifies the defect $\vec d_n$ to ensure $$\int_\Omega F d_n = 0$$ (which does not necessarily hold for $d_n$ in all cases). This assures the above "compatibility condition" for $Fd_n$ which is the right-hand side for the preconditioner: $$ C \vec y_n = F \vec d_n $$

  • The preconditioner $C^{-1}$ is assumed to be an approximation to $A^{-1}$ producing a correction $y_n\in V_h^0$.

  • As a consequence, there is $\int_\Omega u_{n+1}=0$, and thus, there is $u_n\in V_h^0$ for all $n\in\mathbb{N}$.

Example: The boundary-condition filter

We consider the pure homogeneous Poisson problem $$ \begin{align} -\Delta u &= f \\ u_{|\partial\Omega} &= 0 \end{align}$$ on a domain $\Omega\subset\mathbb{R}^d$ for a proper right­hand side $f$ and a function $g:\partial\Omega\to\mathbb{R}$ an appropriate boundary function. The corresponding discrete problem reads: Find $u_h\in V_h$, with $V_h=H^1_0(\Omega)$, such that $$ a(u_h,v_h) = (\nabla u_h, \nabla v_h) = (f,v_h) \qquad\forall\, v_h\in V_h.$$

We consider two meshes $\Omega_h$ and $\Omega_{2h}$ as approximations to the domain $\Omega$. On each of these meshes, there exist a corresponding linear system, solution and right-hand side vector: $$\begin{alignat}{3} A_{h} \vec u_h &= f_h &\qquad&\text{for mesh size $h$}\\ A_{2h} \vec u_{2h} &= f_{2h} &\qquad&\text{for mesh size $2h$} \end{alignat}$$ A two-grid algorithm applies a defect correction using a solution of the coarse grid problem. Here some notations:

Prolongation
Let $P:\vec u_{2h} \mapsto \vec u_h$ be a prolongation operator that maps a solution from the coarse to the fine grid b a proper interpolation. I.e., if there is $u_h=u_{2h}$ in $V_h$, $P$ maps the degrees of freedom $\vec u_{2h}$ to those of $\vec u_h$.
Restriction
Let $R:\vec f_h \mapsto \vec f_{2h}$ be a restriction operator which maps right-hand side vectors from the fine mesh to those on the coarse mesh. Such an operator can be realised as the adjoint of the prolongation $P$.

A coarse grid correction in a two-grid algorithm is a defect correction step with the preconditioner $C^{-1}$ being the solution of a coarse grid problem: $$\vec u_{n+1} = \vec u_{n} + \underbrace{PA_{2h}^{-1}R}_{C^{-1}} \underbrace{(\vec f - A_h \vec u_n)}_{\vec d_h}$$

The restriction operator maps the defect vector to the coarse mesh, $R: \vec d_h\mapsto \vec d_{2h}$. The problem is that $\vec d_{2h}$ is not homogeneous in general:

The degrees of freedom in $\vec d_{2h}$ corresponding to the boundary should be zero, as $\vec d_{2h}$ serves as the right-hand side for $A_{2h}^{-1}$ which solves a homogeneous boundary problem on the coarse mesh.

A simple restriction does not fulfil this; it just does a weighted mean of the DOFs and thus, $\vec d_{2h}$ does not have "homogeneous boundary conditions".

The trick at this point is to introduce a special "boundary condition defect filter" after the restriction, which we call $F$ here: $$\vec u_{n+1} = \vec u_{n} + PA_{2h}^{-1}FR (\vec f - A_h \vec u_n)$$ The filter has to ensure that $F \vec d_{2h}$ fulfils "homogeneous boundary conditions", i.e., $$F d_{2h} \in V_{2h}$$ (with $F d_{2h}$ denoting the continuous counterpart of $F \vec d_{2h}$). This trick gracefully improves the performance of a multilevel-based algorithm.

Filter targets and types

Consider a general linear equation $$Ax=b$$ and a preconditioned defect correction loop $$x_{n+1}=x_n + C^{-1} (b-Ax).$$ This setting covers most iterative algorithms. There are four types of vectors, filters can possibly be applied to:

  • Right-hand side vectors (e.g., $b$)
  • Solution vectors (e.g., $x$ with $Ax=b$)
  • Defect vectors (e.g., $d_n:=(b-Ax_n)$)
  • Correction vectors (e.g., $y_n:=C^{-1} d_n$, to be used in a correction $x_{n+1} = x_n + y_n$)

The linear solver library is implemented in a preconditioning style, i.e., all solvers are preconditioners and all vectors to be preconditioned are defect vectors. Hence, filters in the linear solver are always defect filters. Right-hand side and solution filters are, however, used to prepare a linear system for being solved.

The most important filters provided in the FEAT2 library are:

Boundary-condition filter

Used for the implementation of boundary conditions.

Target Filter Description
Solution vecfil_discreteBCsol Imposes boundary conditions into the solution
RHS vecfil_discreteBCrhs Imposes boundary conditions into the RHS
Defect vecfil_discreteBCdef Forces DOFs on the boundary to zero
Correction vecfil_discreteBCdef Forces DOFs on the boundary to zero
Integral-mean-value-zero filter

Used for imposing $\int_\Omega v_h=0$ for a finite element function $v_h$.

Target Filter Description
Solution vecfil_subvectorToL20 Normalise a FEM function to be in the space $L^2_0$.
RHS vecfil_subvectorToL20 Normalise a FEM function to be in the space $L^2_0$.
Defect vecfil_subvectorToL20 Normalise a FEM function to be in the space $L^2_0$.
Correction vecfil_subvectorToL20 Normalise a FEM function to be in the space $L^2_0$.
$l_1$-zero filter

Used as a fast, approximative alternative to the "Integral-mean-value-zero filter". Equivalent for some finite elements.

Target Filter Description
Solution vecfil_subvectorSmallL1To0 Normalise a vector to be $||.||_{l_1}=0$
RHS vecfil_subvectorSmallL1To0 Normalise a vector to be $||.||_{l_1}=0$
Defect vecfil_subvectorSmallL1To0 Normalise a vector to be $||.||_{l_1}=0$
Correction vecfil_subvectorSmallL1To0 Normalise a vector to be $||.||_{l_1}=0$
Overwrite-entry-by-zero filter

Used for imposing $\int_\Omega v_h=0$ for a finite element function $v_h$. Lumped mass-matrix approach.

Target Filter Description
Solution vecfil_solL1To0ByLMass Normalise the DOFs of a FEM solution $u_h$ such that $\int_{\Omega} u_h = 0$ by using a lumped mass matrix
RHS vecfil_rhsL1To0ByLmass Normalise the DOFs of a FEM RHS function $f_h$ such that $\int_{\Omega} f_h = 0$ by using a lumped mass matrix
Defect vecfil_oneEntryZero Replace one entry(row) of a subvector with zero.
Correction vecfil_oneEntryZero Replace one entry(row) of a subvector with zero.

Matrix filters

Matrix filters are actually not really filters. Routines of this type are matrix modification routines that change the properties of a matrix according to the need of a vector filter. A simple example for this is the boundary condition filter:

Example: The routines vecfil_discreteBCsol and vecfil_discreteBCrhs implement Dirichlet boundary conditions into a solution and defect vector. This is done by overwriting the DOFs on the boundary by the values described by the boundary conditions. The corresponding "matrix filter", matfil_discreteBC replaces the rows corresponding to the DOFs on the boundary by unit vectors. As a consequence, the matrix is the identity operator for boundary DOFs and the solution vector already shows the correct values.

Routines for "filtering" of a matrix are found in matrixfilters.f90.

Filter chains

A filter chain is a set of vector filters to be applied to defect and/or correction vectors during a linear solver. Many iterative solution algorithms support filter chains, while direct solvers like UMFPACK do not. Thus, if the filtering technique should be used, the user has to take care only to choose those solvers for the solution of linear systems that support filter chains.

The filter chain is set up as an array of structures of type t_filterChain. The array is passed to the initialisation routine of a linear solver, which then applies it to its intermediate vectors.

The user can theoretically add all types of filters (for solution, RHS, defect, correction) to a filter chain. If a filter chain is applied to a vector, FEAT2 automatically decides on the type of the vector which filters of a filter chain are applied to the vector: Defect filters are applied only to defect vectors, correction filters are only applied to correction vectors, etc. If a filter does not match the type of the vector it should be applied to, it is not applied.

Example: The following example demonstrates how to set up a filter chain that contains a Dirichlet boundary condition filter vor defect vectors, which is the host common type of usage:

type(t_discreteBC), target :: rdiscreteBC
type(t_linsolNode), pointer :: p_rsolver

integer :: nfilters
type(t_filterChain), dimension(1), target :: RmyFilters
...
call filter_clearFilterChain (RmyFilters,nfilters)
call filter_newFilterDiscBCDef (RmyFilters,nfilters,rdiscreteBC)

call linsol_initBiCGStab (p_rsolver,RfilterChain = RmyFilters)
...

At first, an array RmyFilters is declared which takes hold of all filters. With filter_clearFilterChain, the filter chain is initialised. filter_newFilterDiscBCDef adds a filter for implementing Dirichlet boundary conditions into defect vectors. The filter chain RmyFilters is then passed as parameter to linsol_initBiCGStab which will apply it to its defect vectors.

Example: The following example demonstrates how to set up a filter chain that filters defect vectors to integral mean value zero, used for pure Neumann problems with the Poisson equation

type(t_discreteBC), target :: rdiscreteBC
type(t_linsolNode), pointer :: p_rsolver

integer :: nfilters
type(t_filterChain), dimension(1), target :: RmyFilters
...
call filter_clearFilterChain (RmyFilters,nfilters)
call filter_newFilterToL20 (RmyFilters,nfilters,1)

call linsol_initBiCGStab (p_rsolver,RfilterChain = RmyFilters)
...

Similar as above, an array RmyFilters is set up, but this time with an integral mean value filter. The filter is passed to BiCGStab and will ensure a zero integral mean value for defect vectors during the iteration. If the initial solution and the right-hand side passed to the solver have both integral mean value zero, the final solution will fulfil this as well.

Using filters

How and where filters are used depends on the type of the filter. The following list gives an overview about the usage of each of the filters.

Dirichlet boundary-condition filter

This filter has to be used whereever boundary conditions come into play:

  • vecfil_discreteBCsol: To be used to make sure, Dirichlet boundary conditions are imposed into a solution vector.
  • vecfil_discreteBCrhs: To be used to filter a right-hand side vector in advance of solving a linear system.
  • matfil_discreteBC: To be used to filter a matrix in advance of solving a linear system.
  • vecfil_discreteBCdef: To be used to filter a defect vector, e.g., in a nonlinear loop.
  • The linear solver should apply a filter chain to apply vecfil_discreteBCdef internally on defect vectors.

The following example demonstrates how to prepare a linear system for Dirichlet boundary conditions:

! Solution/RHS/Matrix/boundary conditions
type(t_vectorBlock) :: rx, rrhs
type(t_matrixBlock) :: rmatrix
type(t_discreteBC), target :: rdiscreteBC

! Linear solver
type(t_linsolNode), pointer :: p_rsolver

... (assemble matrix, RHS, prepare BC)

! Clear solution
call lsysbl_clearVector (rx)

! Implement BC into RHS, solution, matrix
call matfil_discreteBC (rmatSystem,rdiscreteBC)
call vecfil_discreteBCrhs (rvecRhs,rdiscreteBC)
call vecfil_discreteBCsol (rvecSol,rdiscreteBC)

! Create a solver
call linsol_initUmfpack4 (p_rsolver)

! Attach the system matrix to the solver.
call linsol_setMatrix(p_rsolver, rmatSystem)

... (solve the system)

If an iterative solution algorithm is to be used, it is advisable to additionally set up a filter chain containing the defect filter. This is demonstrated in the following code:

! Solution/RHS/Matrix/boundary conditions
type(t_vectorBlock) :: rvecSol, rvecRhs
type(t_matrixBlock) :: rmatrix
type(t_discreteBC), target :: rdiscreteBC

! Linear solver
type(t_linsolNode), pointer :: p_rsolver

! Filter chain
integer :: nfilters
type(t_filterChain), dimension(1), target :: RmyFilters

... (assemble matrix, RHS, prepare BC)

! Clear solution
call lsysbl_clearVector (rvecSol)

! Implement BC into RHS, solution, matrix
call matfil_discreteBC (rmatSystem,rdiscreteBC)
call vecfil_discreteBCrhs (rvecRhs,rdiscreteBC)
call vecfil_discreteBCsol (rvecSol,rdiscreteBC)

! Set up a filter chain for filtering the defect
call filter_clearFilterChain (RmyFilters,nfilters)
call filter_newFilterDiscBCDef (RmyFilters,nfilters,rdiscreteBC)

! Create a solver
call linsol_initBiCGStab (p_rsolver,RfilterChain = RmyFilters)

! Attach the system matrix to the solver.
call linsol_setMatrix(p_rsolver, rmatSystem)

... (solve the system)
Zero-integral-mean-value filter

This filter has to be used to filter solution, defect and right-hand side to integral mean value zero:

  • vecfil_subvectorToL20: To be used to make sure, a vector has integral mean value zero
  • The linear solver shall apply a filter chain to apply vecfil_subvectorToL20 internally on defect vectors.

The following example demonstrates how to prepare a linear solver for a pure Neumann Poisson problem:

! Solution/RHS/Matrix/boundary conditions
type(t_vectorBlock) :: rvecSol, rvecRhs
type(t_matrixBlock) :: rmatrix

! Linear solver
type(t_linsolNode), pointer :: p_rsolver

! Filter chain
integer :: nfilters
type(t_filterChain), dimension(1), target :: RmyFilters

... (assemble matrix, RHS)

! Clear solution
call lsysbl_clearVector (rvecSol)

! Implement zero integral mean value into RHS, solution, matrix
call vecfil_subvectorToL20 (rvecRhs,1)
call vecfil_subvectorToL20 (rvecSol,1)

! Set up a filter chain for filtering the defect
call filter_clearFilterChain (RmyFilters,nfilters)
call filter_newFilterToL20 (RmyFilters,nfilters,1)

! Create a solver
call linsol_initBiCGStab (p_rsolver,RfilterChain = RmyFilters)

! Attach the system matrix to the solver.
call linsol_setMatrix(p_rsolver, rmatSystem)

... (solve the system)