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 righthand 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 integralmeanvalue condition into a vector, but their use is also rather advantageous for the implementation of boundary conditions in multigridbased 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, righthand 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 integralmeanvalue 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 righthand 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 matrixvector 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 righthand 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 socalled "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 righthand 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 boundarycondition 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 righthand 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 righthand 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 twogrid 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 righthand 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 twogrid 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 righthand 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 multilevelbased 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} (bAx).$$ This setting covers most iterative algorithms. There are four types of vectors, filters can possibly be applied to:
 Righthand side vectors (e.g., $b$)
 Solution vectors (e.g., $x$ with $Ax=b$)
 Defect vectors (e.g., $d_n:=(bAx_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. Righthand 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:
 Boundarycondition 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  Integralmeanvaluezero 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 "Integralmeanvaluezero 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$  Overwriteentrybyzero filter

Used for imposing $\int_\Omega v_h=0$ for a finite element function $v_h$. Lumped massmatrix 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
andvecfil_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 righthand 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 boundarycondition 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 righthand 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)
 Zerointegralmeanvalue filter

This filter has to be used to filter solution, defect and righthand 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)