Skip to content

Commit

Permalink
Added BlockFEOperators
Browse files Browse the repository at this point in the history
  • Loading branch information
JordiManyer committed Dec 27, 2023
1 parent 3e97b3b commit b1fd9db
Show file tree
Hide file tree
Showing 8 changed files with 283 additions and 5 deletions.
4 changes: 2 additions & 2 deletions src/BlockSolvers/BlockDiagonalSolvers.jl
Original file line number Diff line number Diff line change
@@ -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}
Expand Down
142 changes: 142 additions & 0 deletions src/BlockSolvers/BlockFEOperators.jl
Original file line number Diff line number Diff line change
@@ -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
16 changes: 14 additions & 2 deletions src/BlockSolvers/BlockSolverInterfaces.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,13 +19,25 @@ 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
f :: Function
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
Expand All @@ -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)
Expand All @@ -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)
Expand Down
4 changes: 4 additions & 0 deletions src/BlockSolvers/BlockSolvers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
8 changes: 8 additions & 0 deletions src/GridapSolvers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -25,6 +27,9 @@ module GridapSolvers
export RestrictionOperator, ProlongationOperator
export setup_transfer_operators

# BlockSolvers
export BlockDiagonalSolver

# LinearSolvers
export JacobiLinearSolver
export RichardsonSmoother
Expand All @@ -37,7 +42,10 @@ module GridapSolvers
export IS_MINRESSolver
export IS_SSORSolver

export CGSolver
export MINRESSolver
export GMRESSolver
export FGMRESSolver

# PatchBasedSmoothers
export PatchDecomposition
Expand Down
1 change: 0 additions & 1 deletion src/LinearSolvers/Krylov/FGMRESSolvers.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
56 changes: 56 additions & 0 deletions test/BlockSolvers/BlockDiagonalSolversTests.jl
Original file line number Diff line number Diff line change
@@ -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)
= Measure(Ω,4)

sol(x) = sum(x)
a((u1,u2),(v1,v2)) = (u1v1 + u2v2)*
l((v1,v2)) = (solv1 - solv2)*

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) -> (uv)*
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)
57 changes: 57 additions & 0 deletions test/BlockSolvers/BlockFEOperatorsTests.jl
Original file line number Diff line number Diff line change
@@ -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)
= Measure(Ω,4)

u0 = zero(Y)
sol(x) = sum(x)

# Reference operator
a_ref((u1,u2),(v1,v2)) = (u1v1 + u2v2)*
l_ref((v1,v2)) = (solv1 + solv2)*
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) = (uv)*
l(v) = (solv)*
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)

0 comments on commit b1fd9db

Please sign in to comment.