GridapSolvers.LinearSolvers

Krylov solvers

GridapSolvers.LinearSolvers.CGSolverType
struct CGSolver <: LinearSolver
  ...
end

CGSolver(Pl;maxiter=1000,atol=1e-12,rtol=1.e-6,flexible=false,verbose=0,name="CG")

Left-Preconditioned Conjugate Gradient solver.

source
GridapSolvers.LinearSolvers.MINRESSolverType
struct MINRESSolver <: LinearSolver 
  ...
end

MINRESSolver(m;Pl=nothing,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="MINRES")

MINRES solver, with optional left preconditioners Pl. The preconditioner must be symmetric and positive definite.

source
GridapSolvers.LinearSolvers.GMRESSolverType
struct GMRESSolver <: LinearSolver 
  ...
end

GMRESSolver(m;Pr=nothing,Pl=nothing,restart=false,m_add=1,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="GMRES")

GMRES solver, with optional right and left preconditioners Pr and Pl.

The solver starts by allocating a basis of size m. Then:

  • If restart=true, the basis size is fixed and restarted every m iterations.
  • If restart=false, the basis size is allowed to increase. When full, the solver allocates m_add new basis vectors.
source
GridapSolvers.LinearSolvers.FGMRESSolverType
struct FGMRESSolver <: LinearSolver 
  ...
end

FGMRESSolver(m,Pr;Pl=nothing,restart=false,m_add=1,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="FGMRES")

Flexible GMRES solver, with right-preconditioner Pr and optional left-preconditioner Pl.

The solver starts by allocating a basis of size m. Then:

  • If restart=true, the basis size is fixed and restarted every m iterations.
  • If restart=false, the basis size is allowed to increase. When full, the solver allocates m_add new basis vectors at a time.
source

Richardson iterative solver

Given a linear system $Ax = b$ and a left preconditioner $Pl$, this iterative solver performs the following iteration until the solution converges.

\[ x_{k+1} = x_k + \omega Pl^{-1} (b - A x_k)\]

GridapSolvers.LinearSolvers.RichardsonLinearSolverType
struct RichardsonLinearSolver <: LinearSolver 
  ...
end

RichardsonLinearSolver(ω,maxiter;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")

Richardson Iteration, with an optional left preconditioners Pl.

The relaxation parameter (ω) can either be of type Float64 or Vector{Float64}. This gives flexiblity in relaxation.

source

Smoothers

Given a linear system $Ax = b$, a smoother is an operator S that takes an iterative solution $x_k$ and its residual $r_k = b - A x_k$, and modifies them in place

\[ S : (x_k,r_k) \rightarrow (x_{k+1},r_{k+1})\]

such that $|r_{k+1}| < |r_k|$.

GridapSolvers.LinearSolvers.RichardsonSmootherType
struct RichardsonSmoother{A} <: LinearSolver
  M     :: A
  niter :: Int64
  ω     :: Float64
end

Iterative Richardson smoother. Given a solution x and a residual r, performs niter Richardson iterations with damping parameter ω using the linear solver M. A Richardson iteration is given by:

dx = ω * inv(M) * r
x  = x + dx
r  = r - A * dx

Updates both the solution x and the residual r in place.

source

Preconditioners

Given a linear system $Ax = b$, a preconditioner is an operator that takes an iterative residual $r_k$ and returns a correction $dx_k$.

GridapSolvers.LinearSolvers.GMGLinearSolverFunction
GMGLinearSolver(
  mh::ModelHierarchy,
  matrices::AbstractArray{<:AbstractMatrix},
  prolongations,
  restrictions;
  pre_smoothers   = Fill(RichardsonSmoother(JacobiLinearSolver(),10),num_levels(mh)-1),
  post_smoothers  = pre_smoothers,
  coarsest_solver = LUSolver(),
  mode::Symbol    = :preconditioner,
  maxiter = 100, atol = 1.0e-14, rtol = 1.0e-08, verbose = false,
)

Creates an instance of GMGLinearSolverFromMatrices from the underlying model hierarchy, the system matrices at each level and the transfer operators and smoothers at each level except the coarsest.

The solver has two modes of operation, defined by the kwarg mode:

  • :solver: The GMG solver takes a rhs b and returns a solution x.
  • :preconditioner: The GMG solver takes a residual r and returns a correction dx.
source
GMGLinearSolver(
  mh::ModelHierarchy,
  trials::FESpaceHierarchy,
  tests::FESpaceHierarchy,
  biforms::AbstractArray{<:Function},
  interp,
  restrict;
  pre_smoothers   = Fill(RichardsonSmoother(JacobiLinearSolver(),10),num_levels(mh)-1),
  post_smoothers  = pre_smoothers,
  coarsest_solver = Gridap.Algebra.LUSolver(),
  mode::Symbol    = :preconditioner,
  is_nonlinear    = false,
  maxiter = 100, atol = 1.0e-14, rtol = 1.0e-08, verbose = false,
)

Creates an instance of GMGLinearSolverFromMatrices from the underlying model hierarchy, the trial and test FEspace hierarchies, the weakform lhs at each level and the transfer operators and smoothers at each level except the coarsest.

The solver has two modes of operation, defined by the kwarg mode:

  • :solver: The GMG solver takes a rhs b and returns a solution x.
  • :preconditioner: The GMG solver takes a residual r and returns a correction dx.
source

Other

GridapSolvers.LinearSolvers.NullspaceSolverType
struct NullspaceSolver{A,B}
  solver :: A <: LinearSolver
  kernel :: B <: Union{AbstractMatrix, Vector{<:AbstractVector}}
  constrain_matrix :: Bool = true
end

Solver that computes the solution of a linear system A⋅x = b with a kernel constraint, i.e the returned solution is orthogonal to the provided kernel.

We assume the kernel is provided as a matrix K of dimensions (k,n) with n the dimension of the original system and k the number of kernel vectors.

Two modes of operation are supported:

  • If constrain_matrix is true, the solver will explicitly constrain the system matrix. I.e we consider the augmented system Â⋅x̂ = b̂ where

    • Ak = [A, K'; K, 0]
    • x̂ = [x; λ]
    • b̂ = [b; 0]

    This is often the only option for direct solvers, which require the system matrix to be invertible. This is only supported in serial (its performance bottleneck in parallel).

  • If constrain_matrix is false, the solver preserve the original system and simply project the initial guess and the solution onto the orthogonal complement of the kernel. This option is more suitable for iterative solvers, which usually do not require the system matrix to be invertible (e.g. GMRES, BiCGStab, etc).

source
GridapSolvers.LinearSolvers.CallbackSolverType
CallbackSolver(solver::LinearSolver,callback::Function)

A linear solver that runs a callback function after solving the linear system. The callback function should take the solution vector as its only argument and return nothing, i.e callback(x::AbstractVector) -> nothing.

This structure is useful to add functionality to any linear solver, such as:

  • Logging the solution, residuals, etc.
  • Monitoring properties of the solution, as it's divergence or mean.
  • Modifying the solution in-place after solving the linear system, to apply a correction, for example.
source