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 — Type

abstract 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 — Type

abstract type LinearSolverBlock <: SolverBlock end

SolverBlock that will not be updated between nonlinear iterations.

GridapSolvers.BlockSolvers.NonlinearSolverBlock — Type

abstract 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 — Type

struct 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 — Type

struct 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 — Type

struct MatrixBlock{A} <: LinearSolverBlock

SolverBlock representing an external, independent matrix.

Parameters:

mat::A: The matrix.

GridapSolvers.BlockSolvers.BiformBlock — Type

struct 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 — Type

struct 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, where u lives.
trial::FESpace: The trial space, where du lives.
test::FESpace: The test space, where dv 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 — Function

block_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 — Function

block_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! — Function

block_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 — Function

block_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! — Function

block_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 — Type

struct 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 — Method

function BlockDiagonalSolver(
  blocks  :: AbstractVector{<:SolverBlock},
  solvers :: AbstractVector{<:LinearSolver}
)

Create and instance of BlockDiagonalSolver from its underlying properties.

GridapSolvers.BlockSolvers.BlockDiagonalSolver — Method

function 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 — Method

function 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 — Type

struct 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 — Method

function 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 — Method

function 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.