Corrections and adding test and docs

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/RichardsonLinearSolvers.jl

@@ -3,7 +3,7 @@
       ...
     end
 
-    RichardsonLinearSolver(ω,IterMax;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")
+    RichardsonLinearSolver(ω,maxiter;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")
 
   Richardson Iteration, with an optional left preconditioners `Pl`.
 
@@ -14,14 +14,13 @@
 struct RichardsonLinearSolver<:Gridap.Algebra.LinearSolver
     ω::Union{Vector{Float64},Float64}
     Pl::Union{Gridap.Algebra.LinearSolver,Nothing}
-    IterMax::Int64
     log::ConvergenceLog{Float64}
 end
 
-function RichardsonLinearSolver(ω,IterMax;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")
-    tols = SolverTolerances{Float64}(maxiter=IterMax,atol=atol,rtol=rtol)
+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,IterMax,log)
+    return RichardsonLinearSolver(ω,Pl,log)
 end
 
 struct RichardsonLinearSymbolicSetup<:Gridap.Algebra.SymbolicSetup
@@ -34,7 +33,13 @@ end
 
 function get_solver_caches(solver::RichardsonLinearSolver, A::AbstractMatrix)
     ω = solver.ω
-    return ω
+    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
@@ -76,36 +81,31 @@ end
 
 function Gridap.Algebra.solve!(x::AbstractVector, ns:: RichardsonLinearNumericalSetup, b::AbstractVector)
     solver,A,Pl,caches = ns.solver,ns.A,ns.Pl_ns,ns.caches
-    ω = caches
+    ω, z, r, α = caches
     log = solver.log
-    # Relaxation parameters if ω is of type Float64
-    if isa(ω,Float64)
-        temp = ω
-        ω = temp.*ones(eltype(x), size(x))
-    end
-    # Initialization 
-    z = 0*similar(x)
-    r = 0*similar(x)
+    # Relaxation parameters 
+    α .*= ω
     # residual
-    copyto!(r,b - mul!(r,A,x))
+    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 
-            copyto!(x,x + ω.* z)
-            copyto!(r,b - mul!(r,A,x))
+            x .+= α.* z
+            r .= b
+            mul!(r, A, x, -1, 1)
             done = update!(log,norm(r))
         end
         finalize!(log,norm(r))
-        return x
     else    # Case when no preconditioner is applied
         while !done
-            copyto!(z,r)
-            copyto!(x,x + ω.* z)
-            copyto!(r,b - mul!(r,A,x))
+            x .+= α.* r
+            r .= b
+            mul!(r, A, x, -1, 1)
             done = update!(log,norm(r))
         end
         finalize!(log,norm(r))
-        return x
     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")