From b1fd9db9a4d26fdad33c2d81fd2105cf0946b4dc Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 28 Dec 2023 10:37:31 +1100 Subject: [PATCH] Added BlockFEOperators --- src/BlockSolvers/BlockDiagonalSolvers.jl | 4 +- src/BlockSolvers/BlockFEOperators.jl | 142 ++++++++++++++++++ src/BlockSolvers/BlockSolverInterfaces.jl | 16 +- src/BlockSolvers/BlockSolvers.jl | 4 + src/GridapSolvers.jl | 8 + src/LinearSolvers/Krylov/FGMRESSolvers.jl | 1 - .../BlockSolvers/BlockDiagonalSolversTests.jl | 56 +++++++ test/BlockSolvers/BlockFEOperatorsTests.jl | 57 +++++++ 8 files changed, 283 insertions(+), 5 deletions(-) create mode 100644 src/BlockSolvers/BlockFEOperators.jl create mode 100644 test/BlockSolvers/BlockDiagonalSolversTests.jl create mode 100644 test/BlockSolvers/BlockFEOperatorsTests.jl diff --git a/src/BlockSolvers/BlockDiagonalSolvers.jl b/src/BlockSolvers/BlockDiagonalSolvers.jl index 9b3a30ed..30787edb 100644 --- a/src/BlockSolvers/BlockDiagonalSolvers.jl +++ b/src/BlockSolvers/BlockDiagonalSolvers.jl @@ -1,7 +1,7 @@ struct BlockDiagonalSolver{N,A,B} <: Gridap.Algebra.LinearSolver - blocks :: B - solvers :: C + blocks :: A + solvers :: B function BlockDiagonalSolver( blocks :: AbstractVector{<:SolverBlock}, solvers :: AbstractVector{<:Gridap.Algebra.LinearSolver} diff --git a/src/BlockSolvers/BlockFEOperators.jl b/src/BlockSolvers/BlockFEOperators.jl new file mode 100644 index 00000000..9904586e --- /dev/null +++ b/src/BlockSolvers/BlockFEOperators.jl @@ -0,0 +1,142 @@ + +struct BlockFEOperator{NB,SB,P} <: FEOperator + global_op :: FEOperator + block_ops :: Matrix{<:Union{<:FEOperator,Missing,Nothing}} + is_nonlinear :: Matrix{Bool} +end + +const BlockFESpaceTypes{NB,SB,P} = Union{<:MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}},<:GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}} + +function BlockFEOperator( + res::Matrix{<:Union{<:Function,Missing,Nothing}}, + jac::Matrix{<:Union{<:Function,Missing,Nothing}}, + trial::BlockFESpaceTypes, + test::BlockFESpaceTypes; + kwargs... +) + assem = SparseMatrixAssembler(test,trial) + return BlockFEOperator(res,jac,trial,test,assem) +end + +function BlockFEOperator( + res::Matrix{<:Union{<:Function,Missing,Nothing}}, + jac::Matrix{<:Union{<:Function,Missing,Nothing}}, + trial::BlockFESpaceTypes{NB,SB,P}, + test::BlockFESpaceTypes{NB,SB,P}, + assem::MultiField.BlockSparseMatrixAssembler{NB,NV,SB,P}; + is_nonlinear::Matrix{Bool}=fill(true,(NB,NB)) +) where {NB,NV,SB,P} + @check size(res,1) == size(jac,1) == NB + @check size(res,2) == size(jac,2) == NB + + global_res = residual_from_blocks(NB,SB,P,res) + global_jac = jacobian_from_blocks(NB,SB,P,jac) + global_op = FEOperator(global_res,global_jac,trial,test,assem) + + trial_blocks = blocks(trial) + test_blocks = blocks(test) + assem_blocks = blocks(assem) + block_ops = map(CartesianIndices(res)) do I + if !ismissing(res[I]) && !isnothing(res[I]) + FEOperator(res[I],jac[I],test_blocks[I[1]],trial_blocks[I[2]],assem_blocks[I]) + end + end + return BlockFEOperator{NB,SB,P}(global_op,block_ops,is_nonlinear) +end + +# BlockArrays API + +BlockArrays.blocks(op::BlockFEOperator) = op.block_ops + +# FEOperator API + +FESpaces.get_test(op::BlockFEOperator) = get_test(op.global_op) +FESpaces.get_trial(op::BlockFEOperator) = get_trial(op.global_op) +Algebra.allocate_residual(op::BlockFEOperator,u) = allocate_residual(op.global_op,u) +Algebra.residual(op::BlockFEOperator,u) = residual(op.global_op,u) +Algebra.allocate_jacobian(op::BlockFEOperator,u) = allocate_jacobian(op.global_op,u) +Algebra.jacobian(op::BlockFEOperator,u) = jacobian(op.global_op,u) +Algebra.residual!(b::AbstractVector,op::BlockFEOperator,u) = residual!(b,op.global_op,u) + +function Algebra.jacobian!(A::AbstractBlockMatrix,op::BlockFEOperator{NB},u) where NB + map(blocks(A),blocks(op),op.is_nonlinear) do A,op,nl + if nl + residual!(A,op,u) + end + end + return A +end + +# Private methods + +function residual_from_blocks(NB,SB,P,block_residuals) + function res(u,v) + block_ranges = MultiField.get_block_ranges(NB,SB,P) + block_u = map(r -> (length(r) == 1) ? u[r[1]] : Tuple(u[r]), block_ranges) + block_v = map(r -> (length(r) == 1) ? v[r[1]] : Tuple(v[r]), block_ranges) + block_contrs = map(CartesianIndices(block_residuals)) do I + if !ismissing(block_residuals[I]) && !isnothing(block_residuals[I]) + block_residuals[I](block_u[I[2]],block_v[I[1]]) + end + end + return add_block_contribs(block_contrs) + end + return res +end + +function jacobian_from_blocks(NB,SB,P,block_jacobians) + function jac(u,du,dv) + block_ranges = MultiField.get_block_ranges(NB,SB,P) + block_u = map(r -> (length(r) == 1) ? u[r[1]] : Tuple(u[r]) , block_ranges) + block_du = map(r -> (length(r) == 1) ? du[r[1]] : Tuple(du[r]), block_ranges) + block_dv = map(r -> (length(r) == 1) ? dv[r[1]] : Tuple(dv[r]), block_ranges) + block_contrs = map(CartesianIndices(block_jacobians)) do I + if !ismissing(block_jacobians[I]) && !isnothing(block_jacobians[I]) + block_jacobians[I](block_u[I[2]],block_du[I[2]],block_dv[I[1]]) + end + end + return add_block_contribs(block_contrs) + end + return jac +end + +function add_block_contribs(contrs) + c = contrs[1] + for ci in contrs[2:end] + if !ismissing(ci) && !isnothing(ci) + c = c + ci + end + end + return c +end + +function BlockArrays.blocks(a::MultiField.BlockSparseMatrixAssembler) + return a.block_assemblers +end + +function BlockArrays.blocks(f::MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}) where {NB,SB,P} + block_ranges = MultiField.get_block_ranges(NB,SB,P) + block_spaces = map(block_ranges) do range + (length(range) == 1) ? f[range[1]] : MultiFieldFESpace(f[range]) + end + return block_spaces +end + +function BlockArrays.blocks(f::GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}) where {NB,SB,P} + block_gids = blocks(get_free_dof_ids(f)) + + block_ranges = MultiField.get_block_ranges(NB,SB,P) + block_spaces = map(block_ranges,block_gids) do range, gids + if (length(range) == 1) + space = f[range[1]] + else + global_sf_spaces = f[range] + local_sf_spaces = GridapDistributed.to_parray_of_arrays(map(local_views,global_sf_spaces)) + local_mf_spaces = map(MultiFieldFESpace,local_sf_spaces) + vector_type = GridapDistributed._find_vector_type(local_mf_spaces,gids) + space = MultiFieldFESpace(global_sf_spaces,local_mf_spaces,gids,vector_type) + end + space + end + return block_spaces +end diff --git a/src/BlockSolvers/BlockSolverInterfaces.jl b/src/BlockSolvers/BlockSolverInterfaces.jl index 82f86a90..38bc28ca 100644 --- a/src/BlockSolvers/BlockSolverInterfaces.jl +++ b/src/BlockSolvers/BlockSolverInterfaces.jl @@ -19,6 +19,12 @@ struct BiformBlock <: LinearSolverBlock trial :: FESpace test :: FESpace assem :: Assembler + function BiformBlock(f::Function, + trial::FESpace, + test::FESpace, + assem=SparseMatrixAssembler(trial,test)) + return new(f,trial,test,assem) + end end struct TriformBlock <: NonlinearSolverBlock @@ -26,6 +32,12 @@ struct TriformBlock <: NonlinearSolverBlock trial :: FESpace test :: FESpace assem :: Assembler + function TriformBlock(f::Function, + trial::FESpace, + test::FESpace, + assem=SparseMatrixAssembler(trial,test)) + return new(f,trial,test,assem) + end end # Instantiate blocks @@ -48,7 +60,7 @@ function instantiate_block_cache(block::BiformBlock,mat::AbstractMatrix) end instantiate_block_cache(block::LinearSystemBlock,mat::AbstractMatrix) = mat -function instantiate_block_cache(block::TriformSolverBlock,mat::AbstractMatrix,x::AbstractVector) +function instantiate_block_cache(block::TriformBlock,mat::AbstractMatrix,x::AbstractVector) uh = FEFunction(block.trial,x) f(u,v) = block.f(uh,u,v) return assemble_matrix(f,block.assem,block.trial,block.test) @@ -64,7 +76,7 @@ function update_block_cache!(cache,block::NonlinearSolverBlock,mat::AbstractMatr @abstractmethod end function update_block_cache!(cache,block::LinearSolverBlock,mat::AbstractMatrix,x::AbstractVector) - update_block!(cache,block,mat) + update_block_cache!(cache,block,mat) end function update_block_cache!(cache,block::TriformBlock,mat::AbstractMatrix,x::AbstractVector) diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl index f3b78d65..2a7d0fa9 100644 --- a/src/BlockSolvers/BlockSolvers.jl +++ b/src/BlockSolvers/BlockSolvers.jl @@ -13,10 +13,14 @@ module BlockSolvers using GridapSolvers.MultilevelTools using GridapSolvers.SolverInterfaces + include("BlockFEOperators.jl") + include("BlockSolverInterfaces.jl") include("BlockDiagonalSolvers.jl") include("BlockTriangularSolvers.jl") + export BlockFEOperator + export BlockDiagonalSolver export BlockTriangularSolver end diff --git a/src/GridapSolvers.jl b/src/GridapSolvers.jl index 12870f53..5316d068 100644 --- a/src/GridapSolvers.jl +++ b/src/GridapSolvers.jl @@ -2,11 +2,13 @@ module GridapSolvers include("SolverInterfaces/SolverInterfaces.jl") include("MultilevelTools/MultilevelTools.jl") + include("BlockSolvers/BlockSolvers.jl") include("LinearSolvers/LinearSolvers.jl") include("PatchBasedSmoothers/PatchBasedSmoothers.jl") using GridapSolvers.SolverInterfaces using GridapSolvers.MultilevelTools + using GridapSolvers.BlockSolvers using GridapSolvers.LinearSolvers using GridapSolvers.PatchBasedSmoothers @@ -25,6 +27,9 @@ module GridapSolvers export RestrictionOperator, ProlongationOperator export setup_transfer_operators + # BlockSolvers + export BlockDiagonalSolver + # LinearSolvers export JacobiLinearSolver export RichardsonSmoother @@ -37,7 +42,10 @@ module GridapSolvers export IS_MINRESSolver export IS_SSORSolver + export CGSolver + export MINRESSolver export GMRESSolver + export FGMRESSolver # PatchBasedSmoothers export PatchDecomposition diff --git a/src/LinearSolvers/Krylov/FGMRESSolvers.jl b/src/LinearSolvers/Krylov/FGMRESSolvers.jl index 34da4ae4..31adf264 100644 --- a/src/LinearSolvers/Krylov/FGMRESSolvers.jl +++ b/src/LinearSolvers/Krylov/FGMRESSolvers.jl @@ -70,7 +70,6 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::FGMRESNumericalSetup,b::Abs V, Z, zl, H, g, c, s = caches fill!(V[1],zero(eltype(V[1]))) - fill!(zr,zero(eltype(zr))) fill!(zl,zero(eltype(zl))) # Initial residual diff --git a/test/BlockSolvers/BlockDiagonalSolversTests.jl b/test/BlockSolvers/BlockDiagonalSolversTests.jl new file mode 100644 index 00000000..ad1efac6 --- /dev/null +++ b/test/BlockSolvers/BlockDiagonalSolversTests.jl @@ -0,0 +1,56 @@ +using BlockArrays, LinearAlgebra +using Gridap, Gridap.MultiField, Gridap.Algebra +using PartitionedArrays, GridapDistributed +using GridapSolvers + +np = (2,2) +ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) +end + +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(8,8)) + +reffe = ReferenceFE(lagrangian,Float64,1) +V = FESpace(model,reffe) + +mfs = BlockMultiFieldStyle() +Y = MultiFieldFESpace([V,V];style=mfs) + +Ω = Triangulation(model) +dΩ = Measure(Ω,4) + +sol(x) = sum(x) +a((u1,u2),(v1,v2)) = ∫(u1⋅v1 + u2⋅v2)*dΩ +l((v1,v2)) = ∫(sol⋅v1 - sol⋅v2)*dΩ + +op = AffineFEOperator(a,l,Y,Y) +A, b = get_matrix(op), get_vector(op); + + +# 1) From system blocks +s1 = BlockDiagonalSolver([LUSolver(),LUSolver()]) +ss1 = symbolic_setup(s1,A) +ns1 = numerical_setup(ss1,A) +numerical_setup!(ns1,A) + +x1 = allocate_in_domain(A); fill!(x1,0.0) +solve!(x1,ns1,b) + +# 2) From matrix blocks +s2 = BlockDiagonalSolver([A[Block(1,1)],A[Block(2,2)]],[LUSolver(),LUSolver()]) +ss2 = symbolic_setup(s2,A) +ns2 = numerical_setup(ss2,A) +numerical_setup!(ns2,A) + +x2 = allocate_in_domain(A); fill!(x2,0.0) +solve!(x2,ns2,b) + +# 3) From weakform blocks +aii = (u,v) -> ∫(u⋅v)*dΩ +s3 = BlockDiagonalSolver([aii,aii],[V,V],[V,V],[LUSolver(),LUSolver()]) +ss3 = symbolic_setup(s3,A) +ns3 = numerical_setup(ss3,A) +numerical_setup!(ns3,A) + +x3 = allocate_in_domain(A); fill!(x3,0.0) +solve!(x3,ns3,b) diff --git a/test/BlockSolvers/BlockFEOperatorsTests.jl b/test/BlockSolvers/BlockFEOperatorsTests.jl new file mode 100644 index 00000000..71b94d68 --- /dev/null +++ b/test/BlockSolvers/BlockFEOperatorsTests.jl @@ -0,0 +1,57 @@ +using Test +using BlockArrays, LinearAlgebra +using Gridap, Gridap.MultiField, Gridap.Algebra +using PartitionedArrays, GridapDistributed +using GridapSolvers, GridapSolvers.BlockSolvers + +function same_block_array(A,B) + map(blocks(A),blocks(B)) do A, B + t = map(partition(A),partition(B)) do A, B + A ≈ B + end + reduce(&,t) + end |> all +end + +np = (2,2) +ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) +end + +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(8,8)) + +reffe = ReferenceFE(lagrangian,Float64,1) +V = FESpace(model,reffe) + +mfs = BlockMultiFieldStyle() +Y = MultiFieldFESpace([V,V];style=mfs) + +Ω = Triangulation(model) +dΩ = Measure(Ω,4) + +u0 = zero(Y) +sol(x) = sum(x) + +# Reference operator +a_ref((u1,u2),(v1,v2)) = ∫(u1⋅v1 + u2⋅v2)*dΩ +l_ref((v1,v2)) = ∫(sol⋅v1 + sol⋅v2)*dΩ +res_ref(u,v) = a_ref(u,v) - l_ref(v) +jac_ref(u,du,dv) = a_ref(du,dv) +op_ref = FEOperator(res_ref,jac_ref,Y,Y) +A_ref = jacobian(op_ref,u0) +b_ref = residual(op_ref,u0) + +# Block operator +a(u,v) = ∫(u⋅v)*dΩ +l(v) = ∫(sol⋅v)*dΩ +res(u,v) = a(u,v) - l(v) +jac(u,du,dv) = a(du,dv) + +res_blocks = collect(reshape([res,missing,missing,res],(2,2))) +jac_blocks = collect(reshape([jac,missing,missing,jac],(2,2))) +op = BlockFEOperator(res_blocks,jac_blocks,Y,Y) +A = jacobian(op,u0) +b = residual(op,u0) + +@test same_block_array(A,A_ref) +@test same_block_array(b,b_ref)