Added docs for LinearSolvers

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
 ```

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)))

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}

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

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

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)

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}

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).
 """

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}

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).
 """

src/LinearSolvers/RichardsonSmoothers.jl

@@ -1,35 +1,45 @@
+"""
+    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)
   Mss = symbolic_setup(smoother.M,mat)
   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

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