Vanka smoothers

CHANGELOG.md

@@ -0,0 +1,10 @@
+# Changelog
+
+All notable changes to this project will be documented in this file.
+
+The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
+and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
+
+## Previous versions
+
+A changelog is not maintained for older versions than 0.5.0.

README.md

@@ -19,6 +19,14 @@ Solvers follow a modular design, where most blocks can be combined to produce PD
 - **Geometric Multigrid**: We provide a full-fledged geometric multigrid solver. Highly scalable adaptivity and redistribution of meshes, provided by `p4est` through `GridapP4est.jl`.
 - **PETSc interface**: Full access to PETSc algebraic solvers, through `GridapPETSc.jl`, with full interoperability with the rest of the aforementioned solvers.
 
+## Documentation and examples
+
+A (hopefully) comprehensive documentation is available at [https://gridap.github.io/GridapSolvers.jl/stable/](https://gridap.github.io/GridapSolvers.jl/stable/).
+
+A list of examples is available in the documentation. These include some very well known examples such as the Stokes, Incompressible Navier-Stokes and Darcy problems. The featured scripts are available in `test/Applications`.
+
+An example on how to use the library within an HPC cluster is available in `joss_paper/scalability`. The included example and drivers are used to generate the scalability results in our [JOSS paper](https://doi.org/10.21105/joss.07162).
+
 ## Installation
 
 GridapSolvers is a registered package in the official [Julia package registry](https://github.com/JuliaRegistries/General).  Thus, the installation of GridapSolvers is straight forward using the [Julia's package manager](https://julialang.github.io/Pkg.jl/v1/). Open the Julia REPL, type `]` to enter package mode, and install as follows
@@ -40,4 +48,29 @@ The previous installations steps will setup GridapSolvers to work using Julia's
 
 ## Citation
 
-In order to give credit to the `GridapSolvers` contributors, we ask that you please reference our [JOSS paper](https://joss.theoj.org/papers/10.21105/joss.07162) along with the required citations for [Gridap](https://github.com/gridap/Gridap.jl?tab=readme-ov-file#how-to-cite-gridap).
+In order to give credit to the `GridapSolvers` contributors, we simply ask you to cite the `Gridap` main project as indicated [here](https://github.com/gridap/Gridap.jl#how-to-cite-gridap) and the sub-packages you use as indicated in the corresponding repositories. Please, use the reference below in any publication in which you have made use of `GridapSolvers`:
+
+```
+@article{Manyer2024, 
+  doi = {10.21105/joss.07162}, 
+  url = {https://doi.org/10.21105/joss.07162}, 
+  year = {2024}, 
+  publisher = {The Open Journal}, 
+  volume = {9}, 
+  number = {102}, 
+  pages = {7162}, 
+  author = {Jordi Manyer and Alberto F. Martín and Santiago Badia}, 
+  title = {GridapSolvers.jl: Scalable multiphysics finite element solvers in Julia}, 
+  journal = {Journal of Open Source Software} 
+} 
+```
+
+## Contributing
+
+GridapSolvers is a collaborative project open to contributions. If you want to contribute, please take into account:
+
+  - Before opening a PR with a significant contribution, contact the project administrators by [opening an issue](https://github.com/gridap/GridapSolvers.jl/issues/new) describing what you are willing to implement. Wait for feedback from other community members.
+  - We adhere to the contribution and code-of-conduct instructions of the Gridap.jl project, available [here](https://github.com/gridap/Gridap.jl/blob/master/CONTRIBUTING.md) and [here](https://github.com/gridap/Gridap.jl/blob/master/CODE_OF_CONDUCT.md), resp.  Please, carefully read and follow the instructions in these files.
+  - Open a PR with your contribution.
+
+Want to help? We have [issues waiting for help](https://github.com/gridap/GridapSolvers.jl/labels/help%20wanted). You can start contributing to the GridapSolvers project by solving some of those issues.

docs/src/SolverInterfaces.md

@@ -26,5 +26,6 @@ CurrentModule = GridapSolvers.SolverInterfaces
   init!
   update!
   finalize!
+  set_depth!
   print_message
 ```

src/GridapSolvers.jl

@@ -1,57 +1,58 @@
 module GridapSolvers
 
-  include("SolverInterfaces/SolverInterfaces.jl")
-  include("MultilevelTools/MultilevelTools.jl")
-  include("BlockSolvers/BlockSolvers.jl")
-  include("PatchBasedSmoothers/PatchBasedSmoothers.jl")
-  include("LinearSolvers/LinearSolvers.jl")
-  include("NonlinearSolvers/NonlinearSolvers.jl")
-
-  using GridapSolvers.SolverInterfaces
-  using GridapSolvers.MultilevelTools
-  using GridapSolvers.BlockSolvers
-  using GridapSolvers.LinearSolvers
-  using GridapSolvers.PatchBasedSmoothers
-  using GridapSolvers.NonlinearSolvers
-
-  # MultilevelTools
-  export get_parts, generate_level_parts, generate_subparts
-
-  export ModelHierarchy, CartesianModelHierarchy, P4estCartesianModelHierarchy
-  export num_levels, get_level, get_level_parts
-  export get_model, get_model_before_redist
-
-  export FESpaceHierarchy
-  export get_fe_space, get_fe_space_before_redist
-  export compute_hierarchy_matrices
-
-  export DistributedGridTransferOperator
-  export RestrictionOperator, ProlongationOperator
-  export setup_transfer_operators
-
-  # BlockSolvers
-  export BlockDiagonalSolver
-
-  # LinearSolvers
-  export JacobiLinearSolver
-  export RichardsonSmoother
-  export SymGaussSeidelSmoother
-  export GMGLinearSolver
-  export BlockDiagonalSmoother
-
-  export ConjugateGradientSolver
-  export IS_GMRESSolver
-  export IS_MINRESSolver
-  export IS_SSORSolver
-
-  export CGSolver
-  export MINRESSolver
-  export GMRESSolver
-  export FGMRESSolver
-
-  # PatchBasedSmoothers
-  export PatchDecomposition
-  export PatchFESpace
-  export PatchBasedLinearSolver
+include("SolverInterfaces/SolverInterfaces.jl")
+include("MultilevelTools/MultilevelTools.jl")
+include("BlockSolvers/BlockSolvers.jl")
+include("PatchBasedSmoothers/PatchBasedSmoothers.jl")
+include("LinearSolvers/LinearSolvers.jl")
+include("NonlinearSolvers/NonlinearSolvers.jl")
+
+using GridapSolvers.SolverInterfaces
+using GridapSolvers.MultilevelTools
+using GridapSolvers.BlockSolvers
+using GridapSolvers.LinearSolvers
+using GridapSolvers.PatchBasedSmoothers
+using GridapSolvers.NonlinearSolvers
+
+# MultilevelTools
+export get_parts, generate_level_parts, generate_subparts
+
+export ModelHierarchy, CartesianModelHierarchy, P4estCartesianModelHierarchy
+export num_levels, get_level, get_level_parts
+export get_model, get_model_before_redist
+
+export FESpaceHierarchy
+export get_fe_space, get_fe_space_before_redist
+export compute_hierarchy_matrices
+
+export DistributedGridTransferOperator
+export RestrictionOperator, ProlongationOperator
+export setup_transfer_operators
+
+# BlockSolvers
+export BlockDiagonalSolver
+
+# LinearSolvers
+export JacobiLinearSolver
+export RichardsonSmoother
+export SymGaussSeidelSmoother
+export GMGLinearSolver
+export BlockDiagonalSmoother
+
+export ConjugateGradientSolver
+export IS_GMRESSolver
+export IS_MINRESSolver
+export IS_SSORSolver
+
+export CGSolver
+export MINRESSolver
+export GMRESSolver
+export FGMRESSolver
+
+# PatchBasedSmoothers
+export PatchDecomposition
+export PatchFESpace
+export PatchBasedLinearSolver
+export Closure
 
 end

src/LinearSolvers/CallbackSolver.jl

@@ -0,0 +1,52 @@
+
+struct CallbackSolver{A,B} <: Algebra.LinearSolver
+  solver :: A
+  callback :: B
+
+  function CallbackSolver(solver::LinearSolver,callback::Function)
+    A = typeof(solver)
+    B = typeof(callback)
+    new{A,B}(solver,callback)
+  end
+end
+
+struct CallbackSolverSS{A,B} <: Algebra.SymbolicSetup
+  solver :: A
+  ss :: B
+end
+
+function Algebra.symbolic_setup(solver::CallbackSolver,mat::AbstractMatrix)
+  ss = Algebra.symbolic_setup(solver.solver,mat)
+  return CallbackSolverSS(solver,ss)
+end
+
+struct CallbackSolverNS{A,B} <: Algebra.NumericalSetup
+  solver :: A
+  ns :: B
+end
+
+function Algebra.numerical_setup(ss::CallbackSolverSS,mat::AbstractMatrix)
+  ns = Algebra.numerical_setup(ss.ss,mat)
+  return CallbackSolverNS(ss.solver,ns)
+end
+
+function Algebra.numerical_setup(ss::CallbackSolverSS,mat::AbstractMatrix,x::AbstractVector)
+  ns = Algebra.numerical_setup(ss.ss,mat,x)
+  return CallbackSolverNS(ss.solver,ns)
+end
+
+function Algebra.numerical_setup!(ns::CallbackSolverNS,mat::AbstractMatrix)
+  Algebra.numerical_setup!(ns.ns,mat)
+  return ns
+end
+
+function Algebra.numerical_setup!(ns::CallbackSolverNS,mat::AbstractMatrix,x::AbstractVector)
+  Algebra.numerical_setup!(ns.ns,mat,x)
+  return ns
+end
+
+function Algebra.solve!(x::AbstractVector,ns::CallbackSolverNS,b::AbstractVector)
+  solve!(x,ns.ns,b)
+  ns.solver.callback(x)
+  return x
+end

src/LinearSolvers/GMGLinearSolvers.jl

@@ -415,7 +415,9 @@ function gmg_work_vectors(smatrices::AbstractVector{<:AbstractMatrix})
   return work_vectors
 end
 
-function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVector,Nothing},ns::GMGNumericalSetup)
+function apply_GMG_level!(
+  lev::Integer,xh::Union{<:AbstractVector,Nothing},rh::Union{<:AbstractVector,Nothing},ns::GMGNumericalSetup
+)
   mh = ns.solver.mh
   parts = get_level_parts(mh,lev)
   if i_am_in(parts)

src/LinearSolvers/LinearSolverFromSmoothers.jl

@@ -0,0 +1,50 @@
+
+struct LinearSolverFromSmoother{A} <: Algebra.LinearSolver
+  smoother :: A
+end
+
+struct LinearSolverFromSmootherSS{A,B} <: Algebra.SymbolicSetup
+  smoother :: A
+  smoother_ss :: B
+end
+
+struct LinearSolverFromSmootherNS{A,B,C} <: Algebra.NumericalSetup
+  smoother    :: A
+  smoother_ns :: B
+  caches      :: C
+end
+
+function Gridap.Algebra.symbolic_setup(solver::LinearSolverFromSmoother, mat::AbstractMatrix)
+  ss = symbolic_setup(solver.smoother,mat)
+  return LinearSolverFromSmootherSS(solver.smoother,ss)
+end
+
+function Gridap.Algebra.numerical_setup(ss::LinearSolverFromSmootherSS, mat::AbstractMatrix)
+  ns = numerical_setup(ss.smoother_ss, mat)
+  caches = allocate_in_domain(mat)
+  return LinearSolverFromSmootherNS(ss.smoother,ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup(ss::LinearSolverFromSmootherSS, mat::AbstractMatrix, vec::AbstractVector)
+  ns = numerical_setup(ss.smoother_ss, mat, vec)
+  caches = allocate_in_domain(mat)
+  return LinearSolverFromSmootherNS(ss.smoother,ns,caches)
+end
+
+function Gridap.Algebra.numerical_setup!(ns::LinearSolverFromSmootherNS, mat::AbstractMatrix)
+  numerical_setup!(ns.smoother_ns, mat)
+  return ns
+end
+
+function Gridap.Algebra.numerical_setup!(ns::LinearSolverFromSmootherNS, mat::AbstractMatrix, vec::AbstractVector)
+  numerical_setup!(ns.smoother_ns, mat, vec)
+  return ns
+end
+
+function Gridap.Algebra.solve!(x::AbstractVector,ns::LinearSolverFromSmootherNS, b::AbstractVector)
+  r = ns.caches
+  fill!(x,zero(eltype(x)))
+  copy!(r,b)
+  solve!(x,ns.smoother_ns,r)
+  return x
+end

src/LinearSolvers/LinearSolvers.jl

@@ -18,12 +18,16 @@ using GridapSolvers.MultilevelTools
 using GridapSolvers.SolverInterfaces
 using GridapSolvers.PatchBasedSmoothers
 
+export LinearSolverFromSmoother
 export JacobiLinearSolver
 export RichardsonSmoother
 export SymGaussSeidelSmoother
 export GMGLinearSolver
 export BlockDiagonalSmoother
 export SchurComplementSolver
+export SchwarzLinearSolver
+
+export CallbackSolver
 
 # Wrappers for IterativeSolvers.jl
 export IS_ConjugateGradientSolver
@@ -49,12 +53,18 @@ include("PETSc/ElasticitySolvers.jl")
 include("PETSc/HipmairXuSolvers.jl")
 
 include("IdentityLinearSolvers.jl")
+
+include("LinearSolverFromSmoothers.jl")
 include("JacobiLinearSolvers.jl")
 include("RichardsonSmoothers.jl")
 include("SymGaussSeidelSmoothers.jl")
+
 include("GMGLinearSolvers.jl")
 include("IterativeLinearSolvers.jl")
 include("SchurComplementSolvers.jl")
 include("MatrixSolvers.jl")
+include("SchwarzLinearSolvers.jl")
+
+include("CallbackSolver.jl")
 
 end

src/LinearSolvers/SchwarzLinearSolvers.jl

@@ -0,0 +1,49 @@
+
+# TODO: 
+# - Implement the multiplicative case
+# - Add support for weights/averaging when aggregating the additive case?
+
+struct SchwarzLinearSolver{T,S,A} <: Algebra.LinearSolver
+  local_solvers::A
+  function SchwarzLinearSolver(
+    solver::Union{S,AbstractVector{<:S}};
+    type = :additive
+  ) where S <: Algebra.LinearSolver
+    @check type in (:additive,:multiplicative)
+    @notimplementedif type == :multiplicative # TODO
+    A = typeof(solver)
+    new{type,S,A}(solver)
+  end
+end
+
+struct SchwarzSymbolicSetup{T,S,A,B} <: Algebra.SymbolicSetup
+  solver::SchwarzLinearSolver{T,S,A}
+  local_ss::B
+end
+
+function Algebra.symbolic_setup(s::SchwarzLinearSolver,mat::AbstractMatrix)
+  # TODO: This is where we should compute the comm coloring for the multiplicative case
+  expand(s) = map(m -> s,partition(mat))
+  expand(s::AbstractVector) = s
+
+  local_solvers = expand(s.local_solvers)
+  local_ss = map(symbolic_setup,local_solvers,partition(mat))
+  return SchwarzSymbolicSetup(s,local_ss)
+end
+
+struct SchwarzNumericalSetup{T,S,A,B} <: Algebra.NumericalSetup
+  solver::SchwarzLinearSolver{T,S,A}
+  local_ns::B
+end
+
+function Algebra.numerical_setup(ss::SchwarzSymbolicSetup,mat::PSparseMatrix)
+  local_ns = map(numerical_setup,ss.local_ss,partition(mat))
+  return SchwarzNumericalSetup(ss.solver,local_ns)
+end
+
+function Algebra.solve!(x::PVector,ns::SchwarzNumericalSetup{:additive},b::PVector)
+  map(solve!,partition(x),ns.local_ns,partition(b))
+  assemble!(x) |> wait
+  consistent!(x) |> wait
+  return x
+end