GridapSolvers.BlockSolvers
Many scalable preconditioners for multiphysics problems are based on (possibly partial) block factorizations. This module provides a simple interface to define and use block solvers for block-assembled systems.
Block types
In a same preconditioner, blocks can come from different sources. For example, in a Schur-complement-based preconditioner you might want to solve the eliminated block (which comes from the original matrix), while having an approximation for your Schur complement (which can come from a matrix assembled in your driver, or from a weakform).
For this reason, we define the following abstract interface:
GridapSolvers.BlockSolvers.SolverBlock
— Typeabstract type SolverBlock end
Abstract type representing a block in a block solver. More specifically, it indicates how a block is obtained from the original system matrix.
GridapSolvers.BlockSolvers.LinearSolverBlock
— Typeabstract type LinearSolverBlock <: SolverBlock end
SolverBlock that will not be updated between nonlinear iterations.
GridapSolvers.BlockSolvers.NonlinearSolverBlock
— Typeabstract type NonlinearSolverBlock <: SolverBlock end
SolverBlock that will be updated between nonlinear iterations.
On top of this interface, we provide some useful block implementations:
GridapSolvers.BlockSolvers.LinearSystemBlock
— Typestruct LinearSystemBlock <: LinearSolverBlock
SolverBlock representing a linear (i.e non-updateable) block that is directly taken from the system matrix. This block will not be updated between nonlinear iterations.
GridapSolvers.BlockSolvers.NonlinearSystemBlock
— Typestruct NonlinearSystemBlock <: LinearSolverBlock
I :: Vector{Int}
end
SolverBlock representing a nonlinear (i.e updateable) block that is directly taken from the system matrix. This block will be updated between nonlinear iterations.
Parameters:
- `ids::Vector{Int8}`: Block indices considered when updating the nonlinear block.
By default, all indices are considered.
GridapSolvers.BlockSolvers.MatrixBlock
— Typestruct MatrixBlock{A} <: LinearSolverBlock
SolverBlock representing an external, independent matrix.
Parameters:
mat::A
: The matrix.
GridapSolvers.BlockSolvers.BiformBlock
— Typestruct BiformBlock <: LinearSolverBlock
SolverBlock representing a linear block assembled from a bilinear form. This block will be not updated between nonlinear iterations.
Parameters:
f::Function
: The bilinear form, i.e f(du,dv) = ∫(...)dΩtrial::FESpace
: The trial space.test::FESpace
: The test space.assem::Assembler
: The assembler to use.
GridapSolvers.BlockSolvers.TriformBlock
— Typestruct TriformBlock <: NonlinearSolverBlock
SolverBlock representing a nonlinear block assembled from a trilinear form. This block will be updated between nonlinear iterations.
Parameters:
f::Function
: The trilinear form, i.e f(u,du,dv) = ∫(...)dΩparam::FESpace
: The parameter space, whereu
lives.trial::FESpace
: The trial space, wheredu
lives.test::FESpace
: The test space, wheredv
lives.assem::Assembler
: The assembler to use.ids::Vector{Int8}
: Block indices considered when updating the nonlinear block. By default, all indices are considered.
To create a new type of block, one needs to implement the following implementation (similar to the one for LinearSolver
):
GridapSolvers.BlockSolvers.block_symbolic_setup
— Functionblock_symbolic_setup(block::SolverBlock,solver::LinearSolver,mat::AbstractMatrix)
block_symbolic_setup(block::SolverBlock,solver::LinearSolver,mat::AbstractMatrix,x::AbstractVector)
Given a SolverBlock, returns the symbolic setup associated to the LinearSolver.
GridapSolvers.BlockSolvers.block_numerical_setup
— Functionblock_numerical_setup(block::SolverBlock,ss::BlockSS,mat::AbstractMatrix)
block_numerical_setup(block::SolverBlock,ss::BlockSS,mat::AbstractMatrix,x::AbstractVector)
Given a SolverBlock, returns the numerical setup associated to it.
GridapSolvers.BlockSolvers.block_numerical_setup!
— Functionblock_numerical_setup!(block::SolverBlock,ns::BlockNS,mat::AbstractMatrix)
block_numerical_setup!(block::SolverBlock,ns::BlockNS,mat::AbstractMatrix,x::AbstractVector)
Given a SolverBlock, updates the numerical setup associated to it.
GridapSolvers.BlockSolvers.block_offdiagonal_setup
— Functionblock_offdiagonal_setup(block::SolverBlock,mat::AbstractMatrix)
block_offdiagonal_setup(block::SolverBlock,mat::AbstractMatrix,x::AbstractVector)
Given a SolverBlock, returns the off-diagonal block of associated to it.
GridapSolvers.BlockSolvers.block_offdiagonal_setup!
— Functionblock_offdiagonal_setup!(cache,block::SolverBlock,mat::AbstractMatrix)
block_offdiagonal_setup!(cache,block::SolverBlock,mat::AbstractMatrix,x::AbstractVector)
Given a SolverBlock, updates the off-diagonal block of associated to it.
Block solvers
We can combine blocks to define a block solver. All block solvers take an array of blocks and a vector of solvers for the diagonal blocks (which need to be solved for). We provide two common types of block solvers:
BlockDiagonalSolvers
GridapSolvers.BlockSolvers.BlockDiagonalSolver
— Typestruct BlockDiagonalSolver{N} <: LinearSolver
Solver representing a block-diagonal solver, i.e
│ A11 0 0 │ │ x1 │ │ r1 │
│ 0 A22 0 │ ⋅ │ x2 │ = │ r2 │
│ 0 0 A33 │ │ x3 │ │ r3 │
where N
is the number of diagonal blocks.
Properties:
blocks::AbstractVector{<:SolverBlock}
: Matrix of solver blocks, indicating how each diagonal block of the preconditioner is obtained.solvers::AbstractVector{<:LinearSolver}
: Vector of solvers, one for each diagonal block.
GridapSolvers.BlockSolvers.BlockDiagonalSolver
— Methodfunction BlockDiagonalSolver(
blocks :: AbstractVector{<:SolverBlock},
solvers :: AbstractVector{<:LinearSolver}
)
Create and instance of BlockDiagonalSolver
from its underlying properties.
GridapSolvers.BlockSolvers.BlockDiagonalSolver
— Methodfunction BlockDiagonalSolver(
solvers::AbstractVector{<:LinearSolver};
is_nonlinear::Vector{Bool}=fill(false,length(solvers))
)
Create and instance of BlockDiagonalSolver
where all blocks are extracted from the linear system.
GridapSolvers.BlockSolvers.BlockDiagonalSolver
— Methodfunction BlockDiagonalSolver(
funcs :: AbstractArray{<:Function},
trials :: AbstractArray{<:FESpace},
tests :: AbstractArray{<:FESpace},
solvers :: AbstractArray{<:LinearSolver};
is_nonlinear::Vector{Bool}=fill(false,length(solvers))
)
Create and instance of BlockDiagonalSolver
where all blocks are given by integral forms.
BlockTriangularSolvers
GridapSolvers.BlockSolvers.BlockTriangularSolver
— Typestruct BlockTriangularSolver{T,N} <: LinearSolver
Solver representing a block-triangular (upper/lower) solver, i.e
│ A11 c12⋅A12 c13⋅A13 │ │ x1 │ │ r1 │
│ 0 A22 c23⋅A23 │ ⋅ │ x2 │ = │ r2 │
│ 0 0 A33 │ │ x3 │ │ r3 │
where N
is the number of diagonal blocks and T
is either Val{:upper}
or Val{:lower}
.
Properties:
blocks::AbstractMatrix{<:SolverBlock}
: Matrix of solver blocks, indicating how each block of the preconditioner is obtained.solvers::AbstractVector{<:LinearSolver}
: Vector of solvers, one for each diagonal block.coeffs::AbstractMatrix{<:Real}
: Matrix of coefficients, indicating the contribution of the off-diagonal blocks to the right-hand side of each diagonal. In particular, blocks can be turned off by setting the corresponding coefficient to zero.
GridapSolvers.BlockSolvers.BlockTriangularSolver
— Methodfunction BlockTriangularSolver(
blocks :: AbstractMatrix{<:SolverBlock},
solvers :: AbstractVector{<:LinearSolver},
coeffs = fill(1.0,size(blocks)),
half = :upper
)
Create and instance of BlockTriangularSolver
from its underlying properties. The kwarg half
can have values :upper
or :lower
.
GridapSolvers.BlockSolvers.BlockTriangularSolver
— Methodfunction BlockTriangularSolver(
solvers::AbstractVector{<:LinearSolver};
is_nonlinear::Matrix{Bool}=fill(false,(length(solvers),length(solvers))),
coeffs=fill(1.0,size(is_nonlinear)),
half=:upper
)
Create and instance of BlockTriangularSolver
where all blocks are extracted from the linear system. The kwarg half
can have values :upper
or :lower
.
Staggered FEOperators
GridapSolvers.BlockSolvers.StaggeredFESolver
— Typestruct StaggeredFESolver{NB} <: FESpaces.FESolver
solvers :: Vector{<:Union{LinearSolver,NonlinearSolver}}
end
Solver for staggered problems. See StaggeredFEOperator
for more details.
GridapSolvers.BlockSolvers.StaggeredFEOperator
— Typeabstract type StaggeredFEOperator{NB,SB} <: FESpaces.FEOperator end
Staggered operator, used to solve staggered problems.
We define a staggered problem as a multi-variable non-linear problem where the equation for the k-th variable u_k
only depends on the previous variables u_1,...,u_{k-1}
(and itself).
Such a problem can then be solved by solving each variable sequentially, using the previous variables as input. The most common examples of staggered problems are one-directional coupling problems, where the variables are coupled in a chain-like manner.
Two types of staggered operators are currently supported:
StaggeredAffineFEOperator
: when the k-th equation is linear inu_k
.StaggeredNonlinearFEOperator
: when the k-th equation is non-linear inu_k
.
GridapSolvers.BlockSolvers.StaggeredAffineFEOperator
— Typestruct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB}
...
end
Affine staggered operator, used to solve staggered problems where the k-th equation is linear in u_k
.
Such a problem is formulated by a set of bilinear/linear form pairs:
a_k((u_1,...,u_{k-1}),u_k,v_k) = ∫(...)
l_k((u_1,...,u_{k-1}),v_k) = ∫(...)
than cam be assembled into a set of linear systems:
A_k u_k = b_k
where A_k
and b_k
only depend on the previous variables u_1,...,u_{k-1}
.
GridapSolvers.BlockSolvers.StaggeredNonlinearFEOperator
— Typestruct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB}
...
end
Nonlinear staggered operator, used to solve staggered problems where the k-th equation is nonlinear in u_k
.
Such a problem is formulated by a set of residual/jacobian pairs:
jac_k((u_1,...,u_{k-1}),u_k,du_k,dv_k) = ∫(...)
res_k((u_1,...,u_{k-1}),u_k,v_k) = ∫(...)
Gridap.Algebra.solve!
— Methodsolve(solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB})
solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, cache::Nothing) where NB
solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, cache) where NB