From aec7946d38091e650e44eb2f34d404d27378797e Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 17 Jan 2024 16:24:54 +1100 Subject: [PATCH] Added docs for LinearSolvers --- docs/src/LinearSolvers.md | 41 +++++++++++++- src/LinearSolvers/IterativeLinearSolvers.jl | 36 +++++++++++-- src/LinearSolvers/Krylov/CGSolvers.jl | 8 +++ src/LinearSolvers/Krylov/FGMRESSolvers.jl | 16 +++++- src/LinearSolvers/Krylov/GMRESSolvers.jl | 16 +++++- src/LinearSolvers/Krylov/KrylovUtils.jl | 20 +++++-- src/LinearSolvers/Krylov/MINRESSolvers.jl | 10 +++- src/LinearSolvers/PETSc/ElasticitySolvers.jl | 7 +++ src/LinearSolvers/PETSc/PETScCaches.jl | 14 ++--- src/LinearSolvers/PETSc/PETScUtils.jl | 2 + src/LinearSolvers/RichardsonSmoothers.jl | 57 ++++++++++++-------- src/LinearSolvers/SchurComplementSolvers.jl | 1 + 12 files changed, 186 insertions(+), 42 deletions(-) diff --git a/docs/src/LinearSolvers.md b/docs/src/LinearSolvers.md index b60c6131..3633d380 100644 --- a/docs/src/LinearSolvers.md +++ b/docs/src/LinearSolvers.md @@ -5,6 +5,43 @@ CurrentModule = GridapSolvers.LinearSolvers # GridapSolvers.LinearSolvers -```@autodocs -Modules = [LinearSolvers,] +## Krylov solvers + +```@docs + CGSolver + MINRESSolver + GMRESSolver + FGMRESSolver + krylov_mul! + krylov_residual! +``` + +## Smoothers + +```@docs + RichardsonSmoother +``` + +## Wrappers + +### PETSc + +Building on top of [GridapPETSc.jl](https://github.com/gridap/GridapPETSc.jl), GridapSolvers provides specific solvers for some particularly complex PDEs: + +```@docs + ElasticitySolver + CachedPETScNS + get_dof_coordinates +``` + +### IterativeSolvers.jl + +GridapSolvers provides wrappers for some iterative solvers from the package [IterativeSolvers.jl](https://iterativesolvers.julialinearalgebra.org/dev/): + +```@docs + IterativeLinearSolver + IS_ConjugateGradientSolver + IS_GMRESSolver + IS_MINRESSolver + IS_SSORSolver ``` diff --git a/src/LinearSolvers/IterativeLinearSolvers.jl b/src/LinearSolvers/IterativeLinearSolvers.jl index b3d9d4b6..ab3cd757 100644 --- a/src/LinearSolvers/IterativeLinearSolvers.jl +++ b/src/LinearSolvers/IterativeLinearSolvers.jl @@ -8,13 +8,21 @@ struct SSORIterativeSolverType <: IterativeLinearSolverType end # Constructors """ + struct IterativeLinearSolver <: LinearSolver + ... + end + Wrappers for [IterativeSolvers.jl](https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl) krylov-like iterative solvers. - Currently supported: - - ConjugateGradientSolver - - GMRESSolver - - MINRESSolver + All wrappers take the same kwargs as the corresponding solver in IterativeSolvers.jl. + + The following solvers are available: + + - [`IS_ConjugateGradientSolver`](@ref) + - [`IS_GMRESSolver`](@ref) + - [`IS_MINRESSolver`](@ref) + - [`IS_SSORSolver`](@ref) """ struct IterativeLinearSolver{A} <: Gridap.Algebra.LinearSolver args @@ -28,24 +36,44 @@ end SolverType(::IterativeLinearSolver{T}) where T = T() +""" + IS_ConjugateGradientSolver(;kwargs...) + + Wrapper for the [Conjugate Gradient solver](https://iterativesolvers.julialinearalgebra.org/dev/linear_systems/cg/). +""" function IS_ConjugateGradientSolver(;kwargs...) options = [:statevars,:initially_zero,:Pl,:abstol,:reltol,:maxiter,:verbose,:log] @check all(map(opt -> opt ∈ options,keys(kwargs))) return IterativeLinearSolver(CGIterativeSolverType(),nothing,kwargs) end +""" + IS_GMRESSolver(;kwargs...) + + Wrapper for the [GMRES solver](https://iterativesolvers.julialinearalgebra.org/dev/linear_systems/gmres/). +""" function IS_GMRESSolver(;kwargs...) options = [:initially_zero,:abstol,:reltol,:restart,:maxiter,:Pl,:Pr,:log,:verbose,:orth_meth] @check all(map(opt -> opt ∈ options,keys(kwargs))) return IterativeLinearSolver(GMRESIterativeSolverType(),nothing,kwargs) end +""" + IS_MINRESSolver(;kwargs...) + + Wrapper for the [MINRES solver](https://iterativesolvers.julialinearalgebra.org/dev/linear_systems/minres/). +""" function IS_MINRESSolver(;kwargs...) options = [:initially_zero,:skew_hermitian,:abstol,:reltol,:maxiter,:log,:verbose] @check all(map(opt -> opt ∈ options,keys(kwargs))) return IterativeLinearSolver(MINRESIterativeSolverType(),nothing,kwargs) end +""" + IS_SSORSolver(ω;kwargs...) + + Wrapper for the [SSOR solver](https://iterativesolvers.julialinearalgebra.org/dev/linear_systems/stationary/#SSOR). +""" function IS_SSORSolver(ω::Real;kwargs...) options = [:maxiter] @check all(map(opt -> opt ∈ options,keys(kwargs))) diff --git a/src/LinearSolvers/Krylov/CGSolvers.jl b/src/LinearSolvers/Krylov/CGSolvers.jl index 081da65d..0d3e681e 100644 --- a/src/LinearSolvers/Krylov/CGSolvers.jl +++ b/src/LinearSolvers/Krylov/CGSolvers.jl @@ -1,4 +1,12 @@ +""" + 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. +""" struct CGSolver <: Gridap.Algebra.LinearSolver Pl :: Gridap.Algebra.LinearSolver log :: ConvergenceLog{Float64} diff --git a/src/LinearSolvers/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl index 412f8c24..c4773890 100644 --- a/src/LinearSolvers/Krylov/FGMRESSolvers.jl +++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl @@ -1,5 +1,19 @@ -# FGMRES Solver +""" + 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. +""" struct FGMRESSolver <: Gridap.Algebra.LinearSolver m :: Int restart :: Bool diff --git a/src/LinearSolvers/Krylov/GMRESSolvers.jl b/src/LinearSolvers/Krylov/GMRESSolvers.jl index b7afa0bf..1c095923 100644 --- a/src/LinearSolvers/Krylov/GMRESSolvers.jl +++ b/src/LinearSolvers/Krylov/GMRESSolvers.jl @@ -1,4 +1,18 @@ -# GMRES Solver +""" + 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. +""" struct GMRESSolver <: Gridap.Algebra.LinearSolver m :: Int restart :: Bool diff --git a/src/LinearSolvers/Krylov/KrylovUtils.jl b/src/LinearSolvers/Krylov/KrylovUtils.jl index dcd76b11..0ffb0e19 100644 --- a/src/LinearSolvers/Krylov/KrylovUtils.jl +++ b/src/LinearSolvers/Krylov/KrylovUtils.jl @@ -1,10 +1,16 @@ """ - Computes the Krylov matrix-vector product y = Pl⁻¹⋅A⋅Pr⁻¹⋅x - by solving: + Computes the Krylov matrix-vector product + + `y = Pl⁻¹⋅A⋅Pr⁻¹⋅x` + + by solving + + ``` Pr⋅wr = x wl = A⋅wr Pl⋅y = wl + ``` """ function krylov_mul!(y,A,x,Pr,Pl,wr,wl) solve!(wr,Pr,x) @@ -24,10 +30,16 @@ function krylov_mul!(y,A,x,Pr::Nothing,Pl::Nothing,wr,wl) end """ - Computes the Krylov residual r = Pl⁻¹(A⋅x - b). - by solving: + Computes the Krylov residual + + `r = Pl⁻¹(A⋅x - b)` + + by solving + + ``` w = A⋅x - b Pl⋅r = w + ``` """ function krylov_residual!(r,x,A,b,Pl,w) mul!(w,A,x) diff --git a/src/LinearSolvers/Krylov/MINRESSolvers.jl b/src/LinearSolvers/Krylov/MINRESSolvers.jl index 8feb9fb1..dcc698c5 100644 --- a/src/LinearSolvers/Krylov/MINRESSolvers.jl +++ b/src/LinearSolvers/Krylov/MINRESSolvers.jl @@ -1,4 +1,12 @@ -# MINRES Solver +""" + struct MINRESSolver <: LinearSolver + ... + end + + MINRESSolver(m;Pr=nothing,Pl=nothing,maxiter=100,atol=1e-12,rtol=1.e-6,verbose=false,name="MINRES") + + MINRES solver, with optional right and left preconditioners `Pr` and `Pl`. +""" struct MINRESSolver <: Gridap.Algebra.LinearSolver Pr :: Union{Gridap.Algebra.LinearSolver,Nothing} Pl :: Union{Gridap.Algebra.LinearSolver,Nothing} diff --git a/src/LinearSolvers/PETSc/ElasticitySolvers.jl b/src/LinearSolvers/PETSc/ElasticitySolvers.jl index dc296f3a..42c01346 100644 --- a/src/LinearSolvers/PETSc/ElasticitySolvers.jl +++ b/src/LinearSolvers/PETSc/ElasticitySolvers.jl @@ -1,5 +1,12 @@ """ + struct ElasticitySolver <: LinearSolver + ... + end + + ElasticitySolver(space::FESpace; maxiter=500, atol=1.e-12, rtol=1.e-8) + GMRES + AMG solver, specifically designed for linear elasticity problems. + Follows PETSc's documentation for [PCAMG](https://petsc.org/release/manualpages/PC/PCGAMG.html) and [MatNullSpaceCreateRigidBody](https://petsc.org/release/manualpages/Mat/MatNullSpaceCreateRigidBody.html). """ diff --git a/src/LinearSolvers/PETSc/PETScCaches.jl b/src/LinearSolvers/PETSc/PETScCaches.jl index a2845149..b57f49fa 100644 --- a/src/LinearSolvers/PETSc/PETScCaches.jl +++ b/src/LinearSolvers/PETSc/PETScCaches.jl @@ -1,13 +1,15 @@ """ - Notes on this structure: + struct CachedPETScNS <: NumericalSetup - When converting julia vectors/PVectors to PETSc vectors, we purposely create aliasing - of the vector values. This means we can avoid copying data from one to another before solving, - but we need to be careful about it. + Wrapper around a PETSc NumericalSetup, providing highly efficiend reusable caches: - This structure takes care of this, and makes sure you do not attempt to solve the system - with julia vectors that are not the ones you used to create the solver cache. + When converting julia vectors/PVectors to PETSc vectors, we purposely create aliasing + of the vector values. This means we can avoid copying data from one to another before solving, + but we need to be careful about it. + + This structure takes care of this, and makes sure you do not attempt to solve the system + with julia vectors that are not the ones you used to create the solver cache. """ struct CachedPETScNS{TM,A} ns :: GridapPETSc.PETScLinearSolverNS{TM} diff --git a/src/LinearSolvers/PETSc/PETScUtils.jl b/src/LinearSolvers/PETSc/PETScUtils.jl index 1b8de887..24464808 100644 --- a/src/LinearSolvers/PETSc/PETScUtils.jl +++ b/src/LinearSolvers/PETSc/PETScUtils.jl @@ -2,6 +2,8 @@ # DoF coordinates """ + get_dof_coordinates(space::FESpace) + Given a lagrangian FESpace, returns the physical coordinates of the DoFs, as required by some PETSc solvers. See [PETSc documentation](https://petsc.org/release/manualpages/PC/PCSetCoordinates.html). """ diff --git a/src/LinearSolvers/RichardsonSmoothers.jl b/src/LinearSolvers/RichardsonSmoothers.jl index b343e80d..e72b29b7 100644 --- a/src/LinearSolvers/RichardsonSmoothers.jl +++ b/src/LinearSolvers/RichardsonSmoothers.jl @@ -1,22 +1,32 @@ +""" + struct RichardsonSmoother <: LinearSolver + ... + end + RichardsonSmoother(M::LinearSolver,niter::Int=1,ω::Float64=1.0) -struct RichardsonSmoother{A,B} <: Gridap.Algebra.LinearSolver - M :: Gridap.Algebra.LinearSolver - num_smooth_steps :: A - damping_factor :: B -end - -function RichardsonSmoother(M::Gridap.Algebra.LinearSolver, - num_smooth_steps::Integer=1, - damping_factor::Real=1.0) - A = typeof(num_smooth_steps) - B = typeof(damping_factor) - return RichardsonSmoother{A,B}(M,num_smooth_steps,damping_factor) + Performs `niter` Richardson iterations with relaxation parameter `ω` + using the linear solver `M`. + + Updates both the solution `x` and the residual `r` in place. +""" +struct RichardsonSmoother{A} <: Gridap.Algebra.LinearSolver + M :: A + niter :: Int64 + ω :: Float64 + function RichardsonSmoother( + M::Gridap.Algebra.LinearSolver, + niter::Integer=1, + ω::Real=1.0 + ) + A = typeof(M) + return RichardsonSmoother{A}(M,niter,ω) + end end -struct RichardsonSmootherSymbolicSetup{A} <: Gridap.Algebra.SymbolicSetup - smoother :: RichardsonSmoother - Mss :: A +struct RichardsonSmootherSymbolicSetup{A,B} <: Gridap.Algebra.SymbolicSetup + smoother :: RichardsonSmoother{A} + Mss :: B end function Gridap.Algebra.symbolic_setup(smoother::RichardsonSmoother,mat::AbstractMatrix) @@ -24,12 +34,12 @@ function Gridap.Algebra.symbolic_setup(smoother::RichardsonSmoother,mat::Abstrac return RichardsonSmootherSymbolicSetup(smoother,Mss) end -mutable struct RichardsonSmootherNumericalSetup{A,B,C,D} <: Gridap.Algebra.NumericalSetup - smoother :: RichardsonSmoother - A :: A - Adx :: B - dx :: C - Mns :: D +mutable struct RichardsonSmootherNumericalSetup{A,B,C,D,E} <: Gridap.Algebra.NumericalSetup + smoother :: RichardsonSmoother{A} + A :: B + Adx :: C + dx :: D + Mns :: E end function Gridap.Algebra.numerical_setup(ss::RichardsonSmootherSymbolicSetup, A::AbstractMatrix) @@ -45,11 +55,12 @@ end function Gridap.Algebra.solve!(x::AbstractVector,ns::RichardsonSmootherNumericalSetup,r::AbstractVector) Adx,dx,Mns = ns.Adx,ns.dx,ns.Mns + niter, ω = ns.smoother.niter, ns.smoother.ω iter = 1 - while iter <= ns.smoother.num_smooth_steps + while iter <= niter solve!(dx,Mns,r) - dx .= ns.smoother.damping_factor .* dx + dx .= ω .* dx x .= x .+ dx mul!(Adx, ns.A, dx) r .= r .- Adx diff --git a/src/LinearSolvers/SchurComplementSolvers.jl b/src/LinearSolvers/SchurComplementSolvers.jl index 38f94aa1..d38bdb21 100644 --- a/src/LinearSolvers/SchurComplementSolvers.jl +++ b/src/LinearSolvers/SchurComplementSolvers.jl @@ -7,6 +7,7 @@ where S = D - C A^-1 B """ + struct SchurComplementSolver{T1,T2,T3,T4} <: Gridap.Algebra.LinearSolver A :: T1 B :: T2