Extensions

Project.toml

@@ -9,9 +9,6 @@ BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 Gridap = "56d4f2e9-7ea1-5844-9cf6-b9c51ca7ce8e"
 GridapDistributed = "f9701e48-63b3-45aa-9a63-9bc6c271f355"
-GridapP4est = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9"
-GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
-IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
 LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195"
@@ -21,6 +18,16 @@ Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
 
+[weakdeps]
+GridapP4est = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9"
+GridapPETSc = "bcdc36c2-0c3e-11ea-095a-c9dadae499f1"
+IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"
+
+[extensions]
+GridapP4estExt = "GridapP4est"
+GridapPETScExt = "GridapPETSc"
+IterativeSolversExt = "IterativeSolvers"
+
 [compat]
 AbstractTrees = "0.4"
 BlockArrays = "0.16"
@@ -35,7 +42,7 @@ MPI = "0.20"
 NLsolve = "4.3.0"
 PartitionedArrays = "0.3"
 SparseMatricesCSR = "0.6.7"
-julia = "1.7"
+julia = "1.9"
 
 [extras]
 MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267"

ext/GridapP4estExt/GridapP4estExt.jl

@@ -0,0 +1,37 @@
+module GridapP4estExt
+
+using Gridap
+using GridapP4est
+
+export P4estCartesianModelHierarchy
+
+"""
+    P4estCartesianModelHierarchy(
+      ranks,np_per_level,domain,nc::NTuple{D,<:Integer};
+      num_refs_coarse::Integer = 0,
+      add_labels!::Function = (labels -> nothing),
+      map::Function = identity,
+      isperiodic::NTuple{D,Bool} = Tuple(fill(false,D))
+    ) where D
+
+  Returns a `ModelHierarchy` with a Cartesian model as coarsest level, using GridapP4est.jl. 
+  The i-th level will be distributed among `np_per_level[i]` processors. 
+  The seed model is given by `cmodel = CartesianDiscreteModel(domain,nc)`.
+"""
+function P4estCartesianModelHierarchy(
+  ranks,np_per_level,domain,nc::NTuple{D,<:Integer};
+  num_refs_coarse::Integer = 0,
+  add_labels!::Function = (labels -> nothing),
+  map::Function = identity,
+  isperiodic::NTuple{D,Bool} = Tuple(fill(false,D))
+) where D
+  cparts = generate_subparts(ranks,np_per_level[end])
+  cmodel = CartesianDiscreteModel(domain,nc;map,isperiodic)
+  add_labels!(get_face_labeling(cmodel))
+
+  coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse)
+  mh = ModelHierarchy(ranks,coarse_model,np_per_level)
+  return mh
+end
+
+end # module

ext/GridapPETScExt/GridapPETScExt.jl

@@ -0,0 +1,25 @@
+module GridapPETScExt
+
+using LinearAlgebra
+using SparseArrays
+using SparseMatricesCSR
+
+using Gridap
+using Gridap.Helpers, Gridap.Algebra, Gridap.CellData
+using Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField
+
+using GridapDistributed
+using PartitionedArrays
+using GridapPETSc
+
+using GridapSolvers
+using GridapSolvers.MultilevelTools
+using GridapSolvers.SolverInterfaces
+using GridapSolvers.PatchBasedSmoothers
+
+include("PETScUtils.jl")
+include("PETScCaches.jl")
+include("ElasticitySolvers.jl")
+include("HipmairXuSolvers.jl")
+
+end # module

src/LinearSolvers/PETSc/PETScCaches.jl → ext/GridapPETScExt/PETScCaches.jl

@@ -51,3 +51,38 @@ end
 function Algebra.numerical_setup!(ns::CachedPETScNS,mat::AbstractMatrix,x::AbstractVector)
   numerical_setup!(ns.ns,mat,x)
 end
+
+############################################################################################
+# Optimisations for GridapSolvers + PETSc
+
+function LinearSolvers.gmg_coarse_solver_caches(
+  solver::PETScLinearSolver,
+  smatrices::AbstractVector{<:AbstractMatrix},
+  work_vectors
+)
+  nlevs = num_levels(smatrices)
+  with_level(smatrices,nlevs) do AH
+    _, _, dxH, rH = work_vectors[nlevs-1]
+    cache = CachedPETScNS(
+      numerical_setup(symbolic_setup(solver, AH), AH), dxH, rH
+    )
+    return cache
+  end
+end
+
+function LinearSolvers.gmg_coarse_solver_caches(
+  solver::PETScLinearSolver,
+  smatrices::AbstractVector{<:AbstractMatrix},
+  svectors::AbstractVector{<:AbstractVector},
+  work_vectors
+)
+  nlevs = num_levels(smatrices)
+  with_level(smatrices,nlevs) do AH
+    _, _, dxH, rH = work_vectors[nlevs-1]
+    xH = svectors[nlevs]
+    cache = CachedPETScNS(
+      numerical_setup(symbolic_setup(solver, AH, xH), AH, xH), dxH,rH
+    )
+    return cache
+  end
+end

src/LinearSolvers/IterativeLinearSolvers.jl → ext/IterativeSolversExt.jl

@@ -1,3 +1,14 @@
+module IterativeSolversExt
+
+using LinearAlgebra, PartitionedArrays, SparseArrays
+
+using Gridap, GridapSolvers
+using IterativeSolvers
+
+export IS_ConjugateGradientSolver
+export IS_GMRESSolver
+export IS_MINRESSolver
+export IS_SSORSolver
 
 abstract type IterativeLinearSolverType end
 struct CGIterativeSolverType     <: IterativeLinearSolverType end
@@ -214,3 +225,5 @@ function Gridap.Algebra.solve!(::SSORIterativeSolverType,
   consistent!(x) |> fetch
   return x
 end
+
+end # module

src/LinearSolvers/CallbackSolver.jl

@@ -1,4 +1,18 @@
 
+"""
+    CallbackSolver(solver::LinearSolver,callback::Function)
+
+A linear solver that runs a callback function after solving the linear system. The callback
+function should take the solution vector as its only argument and return nothing, i.e 
+`callback(x::AbstractVector) -> nothing`.
+
+This structure is useful to add functionality to any linear solver, such as: 
+
+- Logging the solution, residuals, etc.
+- Monitoring properties of the solution, as it's divergence or mean. 
+- Modifying the solution in-place after solving the linear system, to apply a correction, 
+  for example.
+"""
 struct CallbackSolver{A,B} <: Algebra.LinearSolver
   solver :: A
   callback :: B

src/LinearSolvers/GMGLinearSolvers.jl

@@ -373,9 +373,6 @@ function gmg_coarse_solver_caches(
   with_level(smatrices,nlevs) do AH
     _, _, dxH, rH = work_vectors[nlevs-1]
     cache = numerical_setup(symbolic_setup(solver, AH), AH)
-    if isa(solver,PETScLinearSolver)
-      cache = CachedPETScNS(cache, dxH, rH)
-    end
     return cache
   end
 end
@@ -391,9 +388,6 @@ function gmg_coarse_solver_caches(
     _, _, dxH, rH = work_vectors[nlevs-1]
     xH = svectors[nlevs]
     cache = numerical_setup(symbolic_setup(solver, AH, xH), AH, xH)
-    if isa(solver,PETScLinearSolver)
-      cache = CachedPETScNS(cache, dxH, rH)
-    end
     return cache
   end
 end

src/LinearSolvers/LinearSolvers.jl

@@ -6,12 +6,10 @@ using LinearAlgebra
 using SparseArrays
 using SparseMatricesCSR
 using BlockArrays
-using IterativeSolvers
 
 using Gridap
 using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField
 using PartitionedArrays
-using GridapPETSc
 
 using GridapDistributed
 using GridapSolvers.MultilevelTools
@@ -46,17 +44,11 @@ include("Krylov/GMRESSolvers.jl")
 include("Krylov/FGMRESSolvers.jl")
 include("Krylov/MINRESSolvers.jl")
 
-include("PETSc/PETScUtils.jl")
-include("PETSc/PETScCaches.jl")
-include("PETSc/ElasticitySolvers.jl")
-include("PETSc/HipmairXuSolvers.jl")
-
 include("IdentityLinearSolvers.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")

src/MultilevelTools/FESpaceHierarchies.jl

@@ -1,10 +1,8 @@
-struct FESpaceHierarchyLevel{A,B,C,D,E}
+struct FESpaceHierarchyLevel{A,B,C}
   level               :: Int
   fe_space            :: A
   fe_space_red        :: B
-  cell_conformity     :: C
-  cell_conformity_red :: D
-  mh_level            :: E
+  mh_level            :: C
 end
 
 """
@@ -25,11 +23,9 @@ get_fe_space_before_redist(sh::FESpaceHierarchy,lev::Int) = get_fe_space_before_
 get_fe_space_before_redist(a::FESpaceHierarchyLevel) = a.fe_space
 
 get_cell_conformity(sh::FESpaceHierarchy,lev::Int) = get_cell_conformity(sh[lev])
-get_cell_conformity(a::FESpaceHierarchyLevel{A,Nothing}) where A = a.cell_conformity
-get_cell_conformity(a::FESpaceHierarchyLevel{A,B}) where {A,B} = a.cell_conformity_red
-
+get_cell_conformity(a::FESpaceHierarchyLevel) = get_cell_conformity(get_fe_space(a))
 get_cell_conformity_before_redist(sh::FESpaceHierarchy,lev::Int) = get_cell_conformity_before_redist(sh[lev])
-get_cell_conformity_before_redist(a::FESpaceHierarchyLevel) = a.cell_conformity
+get_cell_conformity_before_redist(a::FESpaceHierarchyLevel) = get_cell_conformity(get_fe_space_before_redist(a))
 
 get_model(sh::FESpaceHierarchy,level::Integer) = get_model(sh[level])
 get_model(a::FESpaceHierarchyLevel) = get_model(a.mh_level)
@@ -45,56 +41,22 @@ has_refinement(a::FESpaceHierarchyLevel) = has_refinement(a.mh_level)
 
 # Test/Trial FESpaces for ModelHierarchyLevels
 
-function _cell_conformity(
-  model::DiscreteModel,
-  reffe::Tuple{<:Gridap.FESpaces.ReferenceFEName,Any,Any}; 
-  conformity=nothing, kwargs...
-) :: CellConformity
-  basis, reffe_args, reffe_kwargs = reffe
-  cell_reffe = ReferenceFE(model,basis,reffe_args...;reffe_kwargs...)
-  conformity = Conformity(Gridap.Arrays.testitem(cell_reffe),conformity)
-  return CellConformity(cell_reffe,conformity)
-end
-
-function _cell_conformity(
-  model::DiscreteModel,
-  reffe::ReferenceFE; 
-  conformity=nothing, kwargs...
-) :: CellConformity
-  cell_reffe = CompressedArray([reffe],fill(one(Int8),num_cells(model)))
-  conformity = Conformity(reffe,conformity)
-  return CellConformity(cell_reffe,conformity)
-end
-
-function _cell_conformity(
-  model::GridapDistributed.DistributedDiscreteModel,args...;kwargs...
-) :: AbstractVector{<:CellConformity}
-  cell_conformities = map(local_views(model)) do model
-    _cell_conformity(model,args...;kwargs...)
-  end
-  return cell_conformities
-end
-
 function FESpaces.FESpace(mh::ModelHierarchyLevel,args...;kwargs...)
   if has_redistribution(mh)
     cparts, _ = get_old_and_new_parts(mh.red_glue,Val(false))
     Vh     = i_am_in(cparts) ? FESpace(get_model_before_redist(mh),args...;kwargs...) : nothing
     Vh_red = FESpace(get_model(mh),args...;kwargs...)
-    cell_conformity = i_am_in(cparts) ? _cell_conformity(get_model_before_redist(mh),args...;kwargs...) : nothing
-    cell_conformity_red = _cell_conformity(get_model(mh),args...;kwargs...)
   else
     Vh = FESpace(get_model(mh),args...;kwargs...)
     Vh_red = nothing
-    cell_conformity = _cell_conformity(get_model(mh),args...;kwargs...)
-    cell_conformity_red = nothing
   end
-  return FESpaceHierarchyLevel(mh.level,Vh,Vh_red,cell_conformity,cell_conformity_red,mh)
+  return FESpaceHierarchyLevel(mh.level,Vh,Vh_red,mh)
 end
 
 function FESpaces.TrialFESpace(a::FESpaceHierarchyLevel,args...;kwargs...)
   Uh     = !isa(a.fe_space,Nothing) ? TrialFESpace(a.fe_space,args...;kwargs...) : nothing
   Uh_red = !isa(a.fe_space_red,Nothing) ? TrialFESpace(a.fe_space_red,args...;kwargs...) : nothing
-  return FESpaceHierarchyLevel(a.level,Uh,Uh_red,a.cell_conformity,a.cell_conformity_red,a.mh_level)
+  return FESpaceHierarchyLevel(a.level,Uh,Uh_red,a.mh_level)
 end
 
 # Test/Trial FESpaces for ModelHierarchies/FESpaceHierarchy
@@ -132,9 +94,7 @@ function Gridap.MultiField.MultiFieldFESpace(spaces::Vector{<:FESpaceHierarchyLe
   level  = spaces[1].level
   Uh     = all(map(s -> !isa(s.fe_space,Nothing),spaces)) ? MultiFieldFESpace(map(s -> s.fe_space, spaces); kwargs...) : nothing
   Uh_red = all(map(s -> !isa(s.fe_space_red,Nothing),spaces)) ? MultiFieldFESpace(map(s -> s.fe_space_red, spaces); kwargs...) : nothing
-  cell_conformity = map(s -> s.cell_conformity, spaces)
-  cell_conformity_red = map(s -> s.cell_conformity_red, spaces)
-  return FESpaceHierarchyLevel(level,Uh,Uh_red,cell_conformity,cell_conformity_red,first(spaces).mh_level)
+  return FESpaceHierarchyLevel(level,Uh,Uh_red,first(spaces).mh_level)
 end
 
 function Gridap.MultiField.MultiFieldFESpace(spaces::Vector{<:HierarchicalArray};kwargs...)
@@ -158,15 +118,11 @@ function FESpaces.ConstantFESpace(mh::MultilevelTools.ModelHierarchyLevel)
     cparts, _ = get_old_and_new_parts(mh.red_glue,Val(false))
     Vh     = i_am_in(cparts) ? ConstantFESpace(get_model_before_redist(mh)) : nothing
     Vh_red = ConstantFESpace(get_model(mh))
-    cell_conformity = i_am_in(cparts) ? _cell_conformity(get_model_before_redist(mh),reffe) : nothing
-    cell_conformity_red = _cell_conformity(get_model(mh),reffe)
   else
     Vh = ConstantFESpace(get_model(mh))
     Vh_red = nothing
-    cell_conformity = _cell_conformity(get_model(mh),reffe)
-    cell_conformity_red = nothing
   end
-  return FESpaceHierarchyLevel(mh.level,Vh,Vh_red,cell_conformity,cell_conformity_red,mh)
+  return FESpaceHierarchyLevel(mh.level,Vh,Vh_red,mh)
 end
 
 # Computing system matrices

src/MultilevelTools/GridapFixes.jl

@@ -40,3 +40,48 @@ function Arrays.lazy_map(k::Broadcasting{<:Fields.VoidBasisMap},a::Arrays.LazyAr
   lazy_map(a.maps.value,args)
 end
 """
+
+############################################################################################
+# New API: get_cell_conformity
+
+using Gridap.FESpaces: CellConformity, NodeToDofGlue
+using Gridap.TensorValues: change_eltype
+
+function get_cell_polytopes(trian::Triangulation)
+  reffes = get_reffes(trian)
+  polys  = map(get_polytope,reffes)
+  ctypes = get_cell_type(trian)
+  return expand_cell_data(polys,ctypes)
+end
+
+function get_cell_conformity(space::UnstructuredFESpace{V,<:CellConformity}) where V
+  return space.metadata
+end
+
+function get_cell_conformity(space::UnstructuredFESpace{V,<:NodeToDofGlue{T}}) where {V,T}
+  cell_polys = get_cell_polytopes(get_triangulation(space))
+  polys, ctypes = compress_cell_data(cell_polys)
+  reffes = map(p -> LagrangianRefFE(change_eltype(T,eltype(V)),p,1), polys)
+  cell_reffe = expand_cell_data(reffes,ctypes)
+  return CellConformity(cell_reffe,H1Conformity())
+end
+
+for ST in [:TrialFESpace,ZeroMeanFESpace,FESpaceWithConstantFixed]
+  @eval begin 
+    function get_cell_conformity(space::$ST)
+      return get_cell_conformity(space.space)
+    end
+  end
+end
+
+function get_cell_conformity(space::MultiFieldFESpace)
+  map(get_cell_conformity,space)
+end
+
+function get_cell_conformity(space::GridapDistributed.DistributedFESpace)
+  map(get_cell_conformity,local_views(space))
+end
+
+function get_cell_conformity(space::GridapDistributed.DistributedMultiFieldFESpace)
+  map(get_cell_conformity,space)
+end