From e5bf1ee00f91a82bd28d06201fb91a36021456c2 Mon Sep 17 00:00:00 2001 From: Shreyas Date: Thu, 19 Dec 2024 16:55:45 +0100 Subject: [PATCH] Corrections and adding test and docs --- NEWS.md | 1 + docs/src/LinearSolvers.md | 12 +++ src/LinearSolvers/RichardsonLinearSolvers.jl | 46 +++++------ test/LinearSolvers/RichardsonLinearTests.jl | 76 +++++++++++++++++++ .../mpi/RichardsonLinearTests.jl | 10 +++ .../seq/RichardsonLinearTests.jl | 12 +++ test/LinearSolvers/seq/runtests.jl | 1 + 7 files changed, 135 insertions(+), 23 deletions(-) create mode 100644 test/LinearSolvers/RichardsonLinearTests.jl create mode 100644 test/LinearSolvers/mpi/RichardsonLinearTests.jl create mode 100644 test/LinearSolvers/seq/RichardsonLinearTests.jl diff --git a/NEWS.md b/NEWS.md index c21854d3..491f1404 100644 --- a/NEWS.md +++ b/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 diff --git a/docs/src/LinearSolvers.md b/docs/src/LinearSolvers.md index 6be63477..f9f98703 100644 --- a/docs/src/LinearSolvers.md +++ b/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** diff --git a/src/LinearSolvers/RichardsonLinearSolvers.jl b/src/LinearSolvers/RichardsonLinearSolvers.jl index 4922e0be..117dce7e 100644 --- a/src/LinearSolvers/RichardsonLinearSolvers.jl +++ b/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 \ No newline at end of file diff --git a/test/LinearSolvers/RichardsonLinearTests.jl b/test/LinearSolvers/RichardsonLinearTests.jl new file mode 100644 index 00000000..f8f5d242 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/LinearSolvers/mpi/RichardsonLinearTests.jl b/test/LinearSolvers/mpi/RichardsonLinearTests.jl new file mode 100644 index 00000000..7d93f9f1 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/LinearSolvers/seq/RichardsonLinearTests.jl b/test/LinearSolvers/seq/RichardsonLinearTests.jl new file mode 100644 index 00000000..23450b40 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/LinearSolvers/seq/runtests.jl b/test/LinearSolvers/seq/runtests.jl index 3bf20d6c..c925d630 100644 --- a/test/LinearSolvers/seq/runtests.jl +++ b/test/LinearSolvers/seq/runtests.jl @@ -4,3 +4,4 @@ include("KrylovTests.jl") include("IterativeSolversWrappersTests.jl") include("SmoothersTests.jl") include("GMGTests.jl") +include("RichardsonLinearTests.jl")