Skip to content

Commit 32c0e5c

Browse files
JordiManyerJordiManyer
committed
Merge branch 'main' of github.com:gridap/GridapSolvers.jl into extensions
2 parents 18c6d3a + d4fb231 commit 32c0e5c

File tree

9 files changed

+225
-1
lines changed

9 files changed

+225
-1
lines changed

‎.github/workflows/CI.yml

+1-1
Original file line numberDiff line numberDiff line change
@@ -88,7 +88,7 @@ jobs:
8888
cd $SOURCES_DIR
8989
./configure --prefix=$PETSC_INSTALL --with-cc=mpicc --with-cxx=mpicxx --with-fc=mpif90 \
9090
--download-mumps --download-scalapack --download-parmetis --download-metis \
91-
--download-ptscotch --with-debugging --with-x=0 --with-shared=1 \
91+
--download-fblaslapack --download-ptscotch --with-debugging --with-x=0 --with-shared=1 \
9292
--with-mpi=1 --with-64-bit-indices
9393
make
9494
make install

‎NEWS.md

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

1617
## Previous versions
1718

‎docs/src/LinearSolvers.md

+12
Original file line numberDiff line numberDiff line change
@@ -16,6 +16,18 @@ CurrentModule = GridapSolvers.LinearSolvers
1616
krylov_residual!
1717
```
1818

19+
## Richardson iterative solver
20+
21+
Given a linear system ``Ax = b`` and a left **preconditioner** ``Pl``, this iterative solver performs the following iteration until the solution converges.
22+
23+
```math
24+
x_{k+1} = x_k + \omega Pl^{-1} (b - A x_k)
25+
```
26+
27+
```@docs
28+
RichardsonLinearSolver
29+
```
30+
1931
## Smoothers
2032

2133
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

+2
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,7 @@ export GMGLinearSolver
2424
export BlockDiagonalSmoother
2525
export SchurComplementSolver
2626
export SchwarzLinearSolver
27+
export RichardsonLinearSolver
2728

2829
export CallbackSolver
2930

@@ -51,6 +52,7 @@ include("LinearSolverFromSmoothers.jl")
5152
include("JacobiLinearSolvers.jl")
5253
include("RichardsonSmoothers.jl")
5354
include("SymGaussSeidelSmoothers.jl")
55+
include("RichardsonLinearSolvers.jl")
5456

5557
include("GMGLinearSolvers.jl")
5658
include("SchurComplementSolvers.jl")

‎src/LinearSolvers/RichardsonLinearSolvers.jl

+110
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,110 @@
1+
"""
2+
struct RichardsonLinearSolver <: LinearSolver
3+
...
4+
end
5+
6+
RichardsonLinearSolver(ω,maxiter;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")
7+
8+
Richardson Iteration, with an optional left preconditioners `Pl`.
9+
10+
The relaxation parameter (ω) can either be of type Float64 or Vector{Float64}.
11+
This gives flexiblity in relaxation.
12+
"""
13+
struct RichardsonLinearSolver<:Gridap.Algebra.LinearSolver
14+
ω::Union{Vector{Float64},Float64}
15+
Pl::Union{Gridap.Algebra.LinearSolver,Nothing}
16+
log::ConvergenceLog{Float64}
17+
end
18+
19+
function RichardsonLinearSolver(ω,maxiter;Pl=nothing,rtol=1e-10,atol=1e-6,verbose=true,name = "RichardsonLinearSolver")
20+
tols = SolverTolerances{Float64}(maxiter=maxiter,atol=atol,rtol=rtol)
21+
log = ConvergenceLog(name,tols,verbose=verbose)
22+
return RichardsonLinearSolver(ω,Pl,log)
23+
end
24+
25+
struct RichardsonLinearSymbolicSetup<:Gridap.Algebra.SymbolicSetup
26+
solver
27+
end
28+
29+
function Gridap.Algebra.symbolic_setup(solver::RichardsonLinearSolver,A::AbstractMatrix)
30+
return RichardsonLinearSymbolicSetup(solver)
31+
end
32+
33+
function get_solver_caches(solver::RichardsonLinearSolver, A::AbstractMatrix)
34+
ω = solver.ω
35+
z = allocate_in_domain(A)
36+
r = allocate_in_domain(A)
37+
α = allocate_in_domain(A)
38+
fill!(z,0.0)
39+
fill!(r,0.0)
40+
fill!(α,1.0)
41+
return ω, z, r, α
42+
end
43+
44+
mutable struct RichardsonLinearNumericalSetup<:Gridap.Algebra.NumericalSetup
45+
solver
46+
A
47+
Pl_ns
48+
caches
49+
end
50+
51+
function Gridap.Algebra.numerical_setup(ss::RichardsonLinearSymbolicSetup, A::AbstractMatrix)
52+
solver = ss.solver
53+
Pl_ns = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A),A) : nothing
54+
caches = get_solver_caches(solver,A)
55+
return RichardsonLinearNumericalSetup(solver,A,Pl_ns,caches)
56+
end
57+
58+
function Gridap.Algebra.numerical_setup(ss::RichardsonLinearSymbolicSetup, A::AbstractMatrix, x::AbstractVector)
59+
solver = ss.solver
60+
Pl_ns = !isnothing(solver.Pl) ? numerical_setup(symbolic_setup(solver.Pl,A,x),A,x) : nothing
61+
caches = get_solver_caches(solver,A)
62+
return RichardsonLinearNumericalSetup(solver,A,Pl_ns,caches)
63+
end
64+
65+
function Gridap.Algebra.numerical_setup!(ns::RichardsonLinearNumericalSetup, A::AbstractMatrix)
66+
if !isa(ns.Pl_ns,Nothing)
67+
numerical_setup!(ns.Pl_ns,A)
68+
end
69+
ns.A = A
70+
return ns
71+
end
72+
73+
function Gridap.Algebra.numerical_setup!(ns::RichardsonLinearNumericalSetup, A::AbstractMatrix, x::AbstractVector)
74+
if !isa(ns.Pl_ns,Nothing)
75+
numerical_setup!(ns.Pl_ns,A,x)
76+
end
77+
ns.A = A
78+
return ns
79+
end
80+
81+
function Gridap.Algebra.solve!(x::AbstractVector, ns:: RichardsonLinearNumericalSetup, b::AbstractVector)
82+
solver,A,Pl,caches = ns.solver,ns.A,ns.Pl_ns,ns.caches
83+
ω, z, r, α = caches
84+
log = solver.log
85+
# Relaxation parameters
86+
α .*= ω
87+
# residual
88+
r .= b
89+
mul!(r, A, x, -1, 1)
90+
done = init!(log,norm(r))
91+
if !isa(ns.Pl_ns,Nothing) # Case when a preconditioner is applied
92+
while !done
93+
solve!(z, Pl, r) # Apply preconditioner r = PZ
94+
x .+= α.* z
95+
r .= b
96+
mul!(r, A, x, -1, 1)
97+
done = update!(log,norm(r))
98+
end
99+
finalize!(log,norm(r))
100+
else # Case when no preconditioner is applied
101+
while !done
102+
x .+= α.* r
103+
r .= b
104+
mul!(r, A, x, -1, 1)
105+
done = update!(log,norm(r))
106+
end
107+
finalize!(log,norm(r))
108+
end
109+
return x
110+
end

‎test/LinearSolvers/RichardsonLinearTests.jl

+76
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
module RichardsonLinearTests
2+
3+
using Test
4+
using Gridap, Gridap.Algebra
5+
using GridapDistributed
6+
using PartitionedArrays
7+
using IterativeSolvers
8+
9+
using GridapSolvers
10+
using GridapSolvers.LinearSolvers
11+
12+
sol(x) = x[1] + x[2]
13+
f(x) = -Δ(sol)(x)
14+
15+
function test_solver(solver,op,Uh,dΩ)
16+
A, b = get_matrix(op), get_vector(op);
17+
ns = numerical_setup(symbolic_setup(solver,A),A)
18+
19+
x = allocate_in_domain(A); fill!(x,0.0)
20+
solve!(x,ns,b)
21+
22+
u = interpolate(sol,Uh)
23+
uh = FEFunction(Uh,x)
24+
eh = uh - u
25+
E = sum(∫(eh*eh)*dΩ)
26+
@test E < 1.e-6
27+
end
28+
29+
function get_mesh(parts,np)
30+
Dc = length(np)
31+
if Dc == 2
32+
domain = (0,1,0,1)
33+
nc = (8,8)
34+
else
35+
@assert Dc == 3
36+
domain = (0,1,0,1,0,1)
37+
nc = (8,8,8)
38+
end
39+
if prod(np) == 1
40+
model = CartesianDiscreteModel(domain,nc)
41+
else
42+
model = CartesianDiscreteModel(parts,np,domain,nc)
43+
end
44+
return model
45+
end
46+
47+
function main(distribute,np)
48+
parts = distribute(LinearIndices((prod(np),)))
49+
model = get_mesh(parts,np)
50+
51+
order = 1
52+
qorder = order*2 + 1
53+
reffe = ReferenceFE(lagrangian,Float64,order)
54+
Vh = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary")
55+
Uh = TrialFESpace(Vh,sol)
56+
u = interpolate(sol,Uh)
57+
58+
Ω = Triangulation(model)
59+
dΩ = Measure(Ω,qorder)
60+
a(u,v) = ∫(∇(v)⋅∇(u))*dΩ
61+
l(v) = ∫(v⋅f)*dΩ
62+
op = AffineFEOperator(a,l,Uh,Vh)
63+
64+
P = JacobiLinearSolver()
65+
verbose = i_am_main(parts)
66+
67+
# RichardsonLinearSolver with a left preconditioner
68+
solver = RichardsonLinearSolver(0.5,1000; Pl = P, rtol = 1e-8, verbose = verbose)
69+
test_solver(solver,op,Uh,dΩ)
70+
71+
# RichardsonLinearSolver without a preconditioner
72+
solver = RichardsonLinearSolver(0.5,1000; rtol = 1e-8, verbose = verbose)
73+
test_solver(solver,op,Uh,dΩ)
74+
end
75+
76+
end

‎test/LinearSolvers/mpi/RichardsonLinearTests.jl

+10
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,10 @@
1+
module RichardsonLinearTestsMPI
2+
using MPI, PartitionedArrays
3+
include("../RichardsonLinearTests.jl")
4+
5+
with_mpi() do distribute
6+
RichardsonLinearTests.main(distribute,(2,2)) # 2D
7+
RichardsonLinearTests.main(distribute,(2,2,1)) # 3D
8+
end
9+
10+
end

‎test/LinearSolvers/seq/RichardsonLinearTests.jl

+12
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,12 @@
1+
module RichardsonLinearTestsSequential
2+
using PartitionedArrays
3+
include("../RichardsonLinearTests.jl")
4+
5+
with_debug() do distribute
6+
RichardsonLinearTests.main(distribute,(1,1)) # 2D - serial
7+
RichardsonLinearTests.main(distribute,(2,2)) # 2D
8+
RichardsonLinearTests.main(distribute,(1,1,1)) # 3D - serial
9+
RichardsonLinearTests.main(distribute,(2,2,1)) # 3D
10+
end
11+
12+
end

‎test/LinearSolvers/seq/runtests.jl

+1
Original file line numberDiff line numberDiff line change
@@ -4,3 +4,4 @@ include("KrylovTests.jl")
44
include("IterativeSolversWrappersTests.jl")
55
include("SmoothersTests.jl")
66
include("GMGTests.jl")
7+
include("RichardsonLinearTests.jl")

0 commit comments

Comments
 (0)