Merge pull request #87 from shreyas02/main

Adding RichardsonLinearSolver

NEWS.md

@@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Added support for GMG in serial. Since PR[#68](https://github.com/gridap/GridapSolvers.jl/pull/68).
 - Added Vanka-like smoothers in serial. Since PR[#68](https://github.com/gridap/GridapSolvers.jl/pull/68).
 - Added `StaggeredFEOperators` and `StaggeredFESolvers`. Since PR[#84](https://github.com/gridap/GridapSolvers.jl/pull/84).
+- Added `RichardsonLinearSolver`. Since PR[#87](https://github.com/gridap/GridapSolvers.jl/pull/87).
 
 ## Previous versions
 

docs/src/LinearSolvers.md

@@ -16,6 +16,18 @@ CurrentModule = GridapSolvers.LinearSolvers
   krylov_residual!
 ```
 
+## 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.
+
+```math
+  x_{k+1} = x_k + \omega Pl^{-1} (b - A x_k)
+```
+
+```@docs
+  RichardsonLinearSolver
+```
+
 ## 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**

src/LinearSolvers/LinearSolvers.jl

@@ -26,6 +26,7 @@ export GMGLinearSolver
 export BlockDiagonalSmoother
 export SchurComplementSolver
 export SchwarzLinearSolver
+export RichardsonLinearSolver
 
 export CallbackSolver
 
@@ -58,6 +59,7 @@ include("LinearSolverFromSmoothers.jl")
 include("JacobiLinearSolvers.jl")
 include("RichardsonSmoothers.jl")
 include("SymGaussSeidelSmoothers.jl")
+include("RichardsonLinearSolvers.jl")
 
 include("GMGLinearSolvers.jl")
 include("IterativeLinearSolvers.jl")

src/LinearSolvers/RichardsonLinearSolvers.jl

@@ -0,0 +1,110 @@
+"""
+    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.  
+"""
+struct RichardsonLinearSolver<:Gridap.Algebra.LinearSolver
+    ω::Union{Vector{Float64},Float64}
+    Pl::Union{Gridap.Algebra.LinearSolver,Nothing}
+    log::ConvergenceLog{Float64}
+end
+
+function RichardsonLinearSolver(ω,maxiter;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")
+    tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
+    log  = ConvergenceLog(name,tols,verbose=verbose)
+    return RichardsonLinearSolver(ω,Pl,log)
+end
+
+struct RichardsonLinearSymbolicSetup<:Gridap.Algebra.SymbolicSetup
+    solver
+end
+
+function Gridap.Algebra.symbolic_setup(solver::RichardsonLinearSolver,A::AbstractMatrix)
+    return RichardsonLinearSymbolicSetup(solver)
+end
+
+function get_solver_caches(solver::RichardsonLinearSolver, A::AbstractMatrix)
+    ω = solver.ω
+    z = allocate_in_domain(A)
+    r = allocate_in_domain(A)
+    α = allocate_in_domain(A)
+    fill!(z,0.0)
+    fill!(r,0.0)
+    fill!(α,1.0)
+    return ω, z, r, α
+end
+
+mutable struct RichardsonLinearNumericalSetup<:Gridap.Algebra.NumericalSetup
+    solver
+    A
+    Pl_ns
+    caches
+end
+
+function Gridap.Algebra.numerical_setup(ss::RichardsonLinearSymbolicSetup, A::AbstractMatrix)
+    solver = ss.solver
+    Pl_ns = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A),A) : nothing 
+    caches = get_solver_caches(solver,A)
+    return RichardsonLinearNumericalSetup(solver,A,Pl_ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup(ss::RichardsonLinearSymbolicSetup, A::AbstractMatrix, x::AbstractVector)
+    solver = ss.solver
+    Pl_ns = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A,x),A,x) : nothing 
+    caches = get_solver_caches(solver,A)
+    return RichardsonLinearNumericalSetup(solver,A,Pl_ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::RichardsonLinearNumericalSetup, A::AbstractMatrix)
+    if !isa(ns.Pl_ns,Nothing)
+        numerical_setup!(ns.Pl_ns,A)
+    end
+    ns.A = A 
+    return ns 
+end
+
+function Gridap.Algebra.numerical_setup!(ns::RichardsonLinearNumericalSetup, A::AbstractMatrix, x::AbstractVector)
+    if !isa(ns.Pl_ns,Nothing)
+        numerical_setup!(ns.Pl_ns,A,x)
+    end
+    ns.A = A 
+    return ns 
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector, ns:: RichardsonLinearNumericalSetup, b::AbstractVector)
+    solver,A,Pl,caches = ns.solver,ns.A,ns.Pl_ns,ns.caches
+    ω, z, r, α = caches
+    log = solver.log
+    # Relaxation parameters 
+    α .*= ω
+    # residual
+    r .= b
+    mul!(r, A, x, -1, 1)
+    done = init!(log,norm(r))
+    if !isa(ns.Pl_ns,Nothing) # Case when a preconditioner is applied
+        while !done
+            solve!(z, Pl, r) # Apply preconditioner r = PZ 
+            x .+= α.* z
+            r .= b
+            mul!(r, A, x, -1, 1)
+            done = update!(log,norm(r))
+        end
+        finalize!(log,norm(r))
+    else    # Case when no preconditioner is applied
+        while !done
+            x .+= α.* r
+            r .= b
+            mul!(r, A, x, -1, 1)
+            done = update!(log,norm(r))
+        end
+        finalize!(log,norm(r))
+    end
+    return x
+end

test/LinearSolvers/RichardsonLinearTests.jl

@@ -0,0 +1,76 @@
+module RichardsonLinearTests
+
+using Test
+using Gridap, Gridap.Algebra
+using GridapDistributed
+using PartitionedArrays
+using IterativeSolvers
+
+using GridapSolvers
+using GridapSolvers.LinearSolvers
+
+sol(x) = x[1] + x[2]
+f(x)   = -Δ(sol)(x)
+
+function test_solver(solver,op,Uh,dΩ)
+  A, b = get_matrix(op), get_vector(op);
+  ns = numerical_setup(symbolic_setup(solver,A),A)
+
+  x = allocate_in_domain(A); fill!(x,0.0)
+  solve!(x,ns,b)
+
+  u  = interpolate(sol,Uh)
+  uh = FEFunction(Uh,x)
+  eh = uh - u
+  E  = sum(∫(eh*eh)*dΩ)
+  @test E < 1.e-6
+end
+
+function get_mesh(parts,np)
+  Dc = length(np)
+  if Dc == 2
+    domain = (0,1,0,1)
+    nc = (8,8)
+  else
+    @assert Dc == 3
+    domain = (0,1,0,1,0,1)
+    nc = (8,8,8)
+  end
+  if prod(np) == 1
+    model = CartesianDiscreteModel(domain,nc)
+  else
+    model = CartesianDiscreteModel(parts,np,domain,nc)
+  end
+  return model
+end
+
+function main(distribute,np)
+  parts = distribute(LinearIndices((prod(np),)))
+  model = get_mesh(parts,np)
+
+  order  = 1
+  qorder = order*2 + 1
+  reffe  = ReferenceFE(lagrangian,Float64,order)
+  Vh     = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")
+  Uh     = TrialFESpace(Vh,sol)
+  u      = interpolate(sol,Uh)
+
+  Ω      = Triangulation(model)
+  dΩ     = Measure(Ω,qorder)
+  a(u,v) = ∫(∇(v)⋅∇(u))*dΩ
+  l(v)   = ∫(v⋅f)*dΩ
+  op = AffineFEOperator(a,l,Uh,Vh)
+
+  P = JacobiLinearSolver()
+  verbose = i_am_main(parts)
+
+  # RichardsonLinearSolver with a left preconditioner
+  solver = RichardsonLinearSolver(0.5,1000; Pl = P, rtol = 1e-8, verbose = verbose)
+  test_solver(solver,op,Uh,dΩ)
+
+  # RichardsonLinearSolver without a preconditioner
+  solver = RichardsonLinearSolver(0.5,1000; rtol = 1e-8, verbose = verbose)
+  test_solver(solver,op,Uh,dΩ)
+end
+
+end

test/LinearSolvers/mpi/RichardsonLinearTests.jl

@@ -0,0 +1,10 @@
+module RichardsonLinearTestsMPI
+using MPI, PartitionedArrays
+include("../RichardsonLinearTests.jl")
+
+with_mpi() do distribute 
+  RichardsonLinearTests.main(distribute,(2,2))   # 2D
+  RichardsonLinearTests.main(distribute,(2,2,1)) # 3D
+end
+
+end

test/LinearSolvers/seq/RichardsonLinearTests.jl

@@ -0,0 +1,12 @@
+module RichardsonLinearTestsSequential
+using PartitionedArrays
+include("../RichardsonLinearTests.jl")
+
+with_debug() do distribute
+  RichardsonLinearTests.main(distribute,(1,1))   # 2D - serial
+  RichardsonLinearTests.main(distribute,(2,2))   # 2D
+  RichardsonLinearTests.main(distribute,(1,1,1)) # 3D - serial
+  RichardsonLinearTests.main(distribute,(2,2,1)) # 3D 
+end
+
+end

test/LinearSolvers/seq/runtests.jl

@@ -4,3 +4,4 @@ include("KrylovTests.jl")
 include("IterativeSolversWrappersTests.jl")
 include("SmoothersTests.jl")
 include("GMGTests.jl")
+include("RichardsonLinearTests.jl")