From e49618407fb3ba6ced293e1f6f826a9a6426c9e6 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 11 Dec 2024 10:05:32 +1100 Subject: [PATCH 1/8] Added StaggeredFEOperators in serial --- src/BlockSolvers/BlockFEOperators.jl | 156 ++++++++-------- src/BlockSolvers/BlockSolvers.jl | 12 +- src/BlockSolvers/StaggeredFEOperators.jl | 168 ++++++++++++++++++ src/MultilevelTools/GridapFixes.jl | 10 ++ test/BlockSolvers/BlockFEOperatorsTests.jl | 63 ++----- .../BlockSolvers/StaggeredFEOperatorsTests.jl | 55 ++++++ 6 files changed, 341 insertions(+), 123 deletions(-) create mode 100644 src/BlockSolvers/StaggeredFEOperators.jl create mode 100644 test/BlockSolvers/StaggeredFEOperatorsTests.jl diff --git a/src/BlockSolvers/BlockFEOperators.jl b/src/BlockSolvers/BlockFEOperators.jl index 9904586e..8c872baf 100644 --- a/src/BlockSolvers/BlockFEOperators.jl +++ b/src/BlockSolvers/BlockFEOperators.jl @@ -1,47 +1,62 @@ struct BlockFEOperator{NB,SB,P} <: FEOperator - global_op :: FEOperator - block_ops :: Matrix{<:Union{<:FEOperator,Missing,Nothing}} - is_nonlinear :: Matrix{Bool} + global_op :: FEOperator + block_ids :: Vector{CartesianIndex{2}} + block_ops :: Vector{FEOperator} + nonlinear :: Vector{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}}, + res::Vector{<:Union{<:Function,Missing,Nothing}}, jac::Matrix{<:Union{<:Function,Missing,Nothing}}, + args...; nonlinear = fill(true,size(res)) +) + keep(x) = !ismissing(x) && !isnothing(x) + ids = findall(keep, res) + @assert ids == findall(keep, jac) + _res = [res[I] for I in ids] + _jac = [jac[I] for I in ids] + return BlockFEOperator(_res,_jac,ids,args...; nonlinear = nonlinear[ids]) +end + +function BlockFEOperator( + contributions :: Vector{<:Tuple{Any,Function,Function}}, args...; kwargs... +) + ids = [CartesianIndex(c[1]) for c in contributions] + res = [c[2] for c in contributions] + jac = [c[3] for c in contributions] + return BlockFEOperator(res,jac,ids,args...;kwargs...) +end + +function BlockFEOperator( + res::Vector{<:Function}, + jac::Vector{<:Function}, + ids::Vector{CartesianIndex{2}}, trial::BlockFESpaceTypes, - test::BlockFESpaceTypes; + test ::BlockFESpaceTypes; kwargs... ) assem = SparseMatrixAssembler(test,trial) - return BlockFEOperator(res,jac,trial,test,assem) + return BlockFEOperator(res,jac,ids,trial,test,assem;kwargs...) end +# TODO: I think nonlinear should not be a kwarg, since its the whole point of this operator function BlockFEOperator( - res::Matrix{<:Union{<:Function,Missing,Nothing}}, - jac::Matrix{<:Union{<:Function,Missing,Nothing}}, + res::Vector{<:Function}, + jac::Vector{<:Function}, + ids::Vector{CartesianIndex{2}}, 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)) + nonlinear::Vector{Bool}=fill(true,length(res)) ) 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) + ranges = MultiField.get_block_ranges(NB,SB,P) + global_res = residual_from_blocks(ids,ranges,res) + global_jac = jacobian_from_blocks(ids,ranges,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) + block_ops = map(FEOperator,res,jac,blocks(trial),blocks(test),blocks(assem)) + return BlockFEOperator{NB,SB,P}(global_op,block_ids,block_ops,nonlinear) end # BlockArrays API @@ -58,10 +73,11 @@ Algebra.allocate_jacobian(op::BlockFEOperator,u) = allocate_jacobian(op.global_o 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) +function Algebra.jacobian!(A::AbstractBlockMatrix,op::BlockFEOperator,u) + A_blocks = blocks(A) + for (i,I) in enumerate(op.block_ids) + if op.nonlinear[i] + jacobian!(A_blocks[I],op.block_ops[i],u) end end return A @@ -69,47 +85,6 @@ 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 @@ -124,14 +99,13 @@ 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_sf_spaces = 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) @@ -140,3 +114,41 @@ function BlockArrays.blocks(f::GridapDistributed.DistributedMultiFieldFESpace{<: end return block_spaces end + +function liform_from_blocks(ids, ranges, liforms) + function biform(v) + c = DomainContribution() + for (I,lf) in zip(ids, liforms) + vk = v[ranges[I]] + c += lf(uk,vk) + end + return c + end + return biform +end + +function biform_from_blocks(ids, ranges, biforms) + function biform(u,v) + c = DomainContribution() + for (I,bf) in zip(ids, biforms) + uk = u[ranges[I[1]]] + vk = v[ranges[I[2]]] + c += bf(uk,vk) + end + return c + end + return biform +end + +function triform_from_blocks(ids, ranges, triforms) + function triform(u,du,dv) + c = DomainContribution() + for (I,tf) in zip(ids, triforms) + duk = du[ranges[I[1]]] + dvk = dv[ranges[I[2]]] + c += tf(u,duk,dvk) + end + return c + end + return triform +end diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl index 921f01b7..343759fe 100644 --- a/src/BlockSolvers/BlockSolvers.jl +++ b/src/BlockSolvers/BlockSolvers.jl @@ -13,16 +13,24 @@ module BlockSolvers using GridapSolvers.MultilevelTools using GridapSolvers.SolverInterfaces - include("BlockFEOperators.jl") + using GridapDistributed: to_parray_of_arrays, DistributedMultiFieldFESpace + + const MultiFieldFESpaceTypes = Union{<:MultiFieldFESpace,<:GridapDistributed.DistributedMultiFieldFESpace} + const BlockFESpaceTypes{NB,SB,P} = Union{<:MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}},<:GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}} include("BlockSolverInterfaces.jl") include("BlockDiagonalSolvers.jl") include("BlockTriangularSolvers.jl") - export BlockFEOperator + include("BlockFEOperators.jl") + include("StaggeredFEOperators.jl") export MatrixBlock, LinearSystemBlock, NonlinearSystemBlock, BiformBlock, TriformBlock export BlockDiagonalSolver export BlockTriangularSolver + + export BlockFEOperator + export StaggeredFEOperator, StaggeredAffineFEOperator, StaggeredFESolver + end diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl new file mode 100644 index 00000000..d460d227 --- /dev/null +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -0,0 +1,168 @@ + +abstract type StaggeredFEOperator{NB,SB} <: FESpaces.FEOperator end + +function MultiField.get_block_ranges(::StaggeredFEOperator{NB,SB}) where {NB,SB} + MultiField.get_block_ranges(NB,SB,Tuple(1:sum(SB))) +end + +function get_solution(op::StaggeredFEOperator{NB,SB}, xh::MultiFieldFEFunction, k) where {NB,SB} + r = MultiField.get_block_ranges(op)[k] + if length(r) == 1 + xh_k = xh[r[1]] + else + fv_k = blocks(xh.free_values)[k] + xh_k = MultiFieldFEFunction(fv_k, op.trials[k], xh.single_fe_functions[r]) + end + return xh_k +end + +# StaggeredFESolver + +struct StaggeredFESolver{NB} <: FESpaces.FESolver + solvers :: Vector{<:NonlinearSolver} + function StaggeredFESolver(solvers::Vector{<:NonlinearSolver}) + NB = length(solvers) + new{NB}(solvers) + end +end + +function Algebra.solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, ::Nothing) where NB + solvers = solver.solvers + xhs, caches, operators = (), (), () + for k in 1:NB + xh_k = get_solution(op,xh,k) + op_k = get_operator(op,xhs,k) + xh_k, cache_k = solve!(xh_k,solvers[k],op_k,nothing) + xhs, caches, operators = (xhs...,xh_k), (caches...,cache_k), (operators...,op_k) + end + return xh, (caches,operators) +end + +function Algebra.solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, cache) where NB + last_caches, last_operators = cache + solvers = solver.solvers + xhs, caches, operators = (), (), () + for k in 1:NB + xh_k = get_solution(op,xh,k) + op_k = get_operator!(last_operators[k],op,xhs,k) + xh_k, cache_k = solve!(xh_k,solvers[k],op_k,last_caches[k]) + xhs, caches, operators = (xhs...,xh_k), (caches...,cache_k), (operators...,op_k) + end + return xh, (caches,operators) +end + +# StaggeredAffineFEOperator + +struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} + biforms :: Vector{<:Function} + liforms :: Vector{<:Function} + trials :: Vector{<:FESpace} + tests :: Vector{<:FESpace} + assems :: Vector{<:Assembler} + trial :: BlockFESpaceTypes{NB,SB} + test :: BlockFESpaceTypes{NB,SB} + + function StaggeredAffineFEOperator( + biforms::Vector{<:Function}, + liforms::Vector{<:Function}, + trials::Vector{<:FESpace}, + tests::Vector{<:FESpace}, + assems::Vector{<:Assembler} + ) + @assert length(biforms) == length(liforms) == length(trials) == length(tests) == length(assems) + trial = combine_fespaces(trials) + test = combine_fespaces(tests) + NB, SB = length(trials), Tuple(map(num_fields,trials)) + new{NB,SB}(biforms,liforms,trials,tests,assems,trial,test) + end +end + +function StaggeredAffineFEOperator( + biforms::Vector{<:Function}, + liforms::Vector{<:Function}, + trials::Vector{<:FESpace}, + tests::Vector{<:FESpace} +) + assems = map(SparseMatrixAssembler,tests,trials) + return StaggeredAffineFEOperator(biforms,liforms,trials,tests,assems) +end + +FESpaces.get_trial(op::StaggeredAffineFEOperator) = op.trial +FESpaces.get_test(op::StaggeredAffineFEOperator) = op.test + +# Utils + +MultiField.num_fields(space::FESpace) = 1 + +# TODO: We could reuse gids in distributed +function combine_fespaces(spaces::Vector{<:FESpace}) + NB = length(spaces) + SB = Tuple(map(num_fields,spaces)) + + sf_spaces = vcat(map(split_fespace,spaces)...) + MultiFieldFESpace(sf_spaces; style = BlockMultiFieldStyle(NB,SB)) +end + +split_fespace(space::FESpace) = [space] +split_fespace(space::MultiFieldFESpaceTypes) = space.spaces + +function get_operator(op::StaggeredFEOperator{NB}, xhs, k) where NB + @assert NB >= k + a(uk,vk) = op.biforms[k](xhs,uk,vk) + l(vk) = op.liforms[k](xhs,vk) + # uhd = zero(op.trials[k]) + # A, b = assemble_matrix_and_vector(a,l,op.assems[k],op.trials[k],op.tests[k],uhd) + # return AffineOperator(A,b) + return AffineFEOperator(a,l,op.trials[k],op.tests[k],op.assems[k]) +end + +function get_operator!(op_k::AffineFEOperator, op::StaggeredFEOperator{NB}, xhs, k) where NB + @assert NB >= k + A, b = get_matrix(op_k), get_vector(op_k) + a(uk,vk) = op.biforms[k](xhs,uk,vk) + l(vk) = op.liforms[k](xhs,vk) + uhd = zero(op.trials[k]) + assemble_matrix_and_vector!(a,l,A,b,op.assems[k],op.trials[k],op.tests[k],uhd) + return AffineOperator(A,b) +end + +############################################################################################ +# StaggeredFEOperator + +struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} + jacobians :: Vector{<:Function} + residuals :: Vector{<:Function} + trials :: Vector{<:FESpace} + tests :: Vector{<:FESpace} + assems :: Vector{<:Assembler} + trial :: BlockFESpaceTypes{NB,SB} + test :: BlockFESpaceTypes{NB,SB} + + function StaggeredNonlinearFEOperator( + res::Vector{<:Function}, + jac::Vector{<:Function}, + trials::Vector{<:FESpace}, + tests::Vector{<:FESpace}, + assems::Vector{<:Assembler} + ) + @assert length(res) == length(jac) == length(trials) == length(tests) == length(assems) + trial = combine_fespaces(trials) + test = combine_fespaces(tests) + NB, SB = length(trials), Tuple(map(num_fields,trials)) + new{NB,SB}(res,jac,trials,tests,assems,trial,test) + end +end + +# TODO: Can be compute jacobians from residuals? +function StaggeredNonlinearFEOperator( + res::Vector{<:Function}, + jac::Vector{<:Function}, + trials::Vector{<:FESpace}, + tests::Vector{<:FESpace} +) + assems = map(SparseMatrixAssembler,tests,trials) + return StaggeredNonlinearFEOperator(res,jac,trials,tests,assems) +end + +FESpaces.get_trial(op::StaggeredNonlinearFEOperator) = op.trial +FESpaces.get_test(op::StaggeredNonlinearFEOperator) = op.test diff --git a/src/MultilevelTools/GridapFixes.jl b/src/MultilevelTools/GridapFixes.jl index 5db5568f..6e239f72 100644 --- a/src/MultilevelTools/GridapFixes.jl +++ b/src/MultilevelTools/GridapFixes.jl @@ -85,3 +85,13 @@ end function get_cell_conformity(space::GridapDistributed.DistributedMultiFieldFESpace) map(get_cell_conformity,space) end + +# Assembly + +# For some reason this signature is missing? + +function FESpaces.assemble_matrix_and_vector(f::Function,b::Function,a::Assembler,U::FESpace,V::FESpace,uhd) + v = get_fe_basis(V) + u = get_trial_fe_basis(U) + FESpaces.assemble_matrix_and_vector(a,FESpaces.collect_cell_matrix_and_vector(U,V,f(u,v),b(v),uhd)) +end diff --git a/test/BlockSolvers/BlockFEOperatorsTests.jl b/test/BlockSolvers/BlockFEOperatorsTests.jl index 71b94d68..e982f86c 100644 --- a/test/BlockSolvers/BlockFEOperatorsTests.jl +++ b/test/BlockSolvers/BlockFEOperatorsTests.jl @@ -4,54 +4,19 @@ 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) + +model = CartesianDiscreteModel((0,1,0,1),(4,4)) + +order = 1 +reffe = ReferenceFE(lagrangian,Float64,order) 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) +Ω = Triangulation(model) +dΩ = Measure(Ω,3*order) +jac(u,du,dv) = ∫(u * du * dv)dΩ +res(u,dv) = ∫(u * dv)dΩ + + + + + diff --git a/test/BlockSolvers/StaggeredFEOperatorsTests.jl b/test/BlockSolvers/StaggeredFEOperatorsTests.jl new file mode 100644 index 00000000..d3a7f826 --- /dev/null +++ b/test/BlockSolvers/StaggeredFEOperatorsTests.jl @@ -0,0 +1,55 @@ +module StaggeredFEOperatorsTests + +using Test +using Gridap, GridapDistributed, PartitionedArrays +using Gridap.MultiField +using BlockArrays +using GridapSolvers +using GridapSolvers.BlockSolvers + +model = CartesianDiscreteModel((0,1,0,1),(4,4)) + +order = 1 +reffe = ReferenceFE(lagrangian,Float64,order) +V = FESpace(model,reffe;dirichlet_tags="boundary") + +sol = [x -> x[1], x -> x[2], x -> x[1] + x[2], x -> x[1] - x[2]] +U1 = TrialFESpace(V,sol[1]) +U2 = TrialFESpace(V,sol[2]) +U3 = TrialFESpace(V,sol[3]) +U4 = TrialFESpace(V,sol[4]) + +mfs = BlockMultiFieldStyle(3,(1,2,1)) +X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs) +Y = MultiFieldFESpace([V,V,V,V];style=mfs) + +Ω = Triangulation(model) +dΩ = Measure(Ω,3*order) + +a1((),u1,v1) = ∫(u1 * v1)dΩ +l1((),v1) = ∫(sol[1] * v1)dΩ + +a2((u1,),(u2,u3),(v2,v3)) = ∫(u1 * u2 * v2)dΩ + ∫(u3 * v3)dΩ +l2((u1,),(v2,v3)) = ∫(sol[2] * u1 * v2)dΩ + ∫(sol[3] * v3)dΩ + +a3((u1,(u2,u3)),u4,v4) = ∫((u1 + u2) * u4 * v4)dΩ +l3((u1,(u2,u3)),v4) = ∫(sol[4] * (u1 + u2) * v4)dΩ + +UB1, VB1 = U1, V +UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V]) +UB3, VB3 = U4, V +op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],[UB1,UB2,UB3],[VB1,VB2,VB3]) + +solver = StaggeredFESolver([LUSolver(),LUSolver(),LUSolver()]) +xh = solve(solver,op) + +xh_exact = interpolate(sol,X) + +for k in 1:4 + eh_k = xh[k] - xh_exact[k] + e_k = sum(∫(eh_k * eh_k)dΩ) + println("Error in field $k: $e_k") + @test e_k < 1e-10 +end + +end # module \ No newline at end of file From a30f62c7a3baed98d46f0fa6cb63a05b15481228 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 11 Dec 2024 16:09:31 +1100 Subject: [PATCH 2/8] Added distributed StaggeredOperators --- src/BlockSolvers/BlockSolvers.jl | 50 ++++++++++--------- src/BlockSolvers/StaggeredFEOperators.jl | 24 +++++++-- src/SolverInterfaces/GridapExtras.jl | 44 ++++++++++++++++ src/SolverInterfaces/SolverInterfaces.jl | 4 ++ .../BlockSolvers/StaggeredFEOperatorsTests.jl | 12 +++-- 5 files changed, 101 insertions(+), 33 deletions(-) diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl index 343759fe..17ced70f 100644 --- a/src/BlockSolvers/BlockSolvers.jl +++ b/src/BlockSolvers/BlockSolvers.jl @@ -1,36 +1,38 @@ module BlockSolvers - 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 GridapDistributed +using LinearAlgebra +using SparseArrays +using SparseMatricesCSR +using BlockArrays +using IterativeSolvers - using GridapSolvers.MultilevelTools - using GridapSolvers.SolverInterfaces +using Gridap +using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField +using PartitionedArrays +using GridapDistributed - using GridapDistributed: to_parray_of_arrays, DistributedMultiFieldFESpace +using GridapSolvers.MultilevelTools +using GridapSolvers.SolverInterfaces - const MultiFieldFESpaceTypes = Union{<:MultiFieldFESpace,<:GridapDistributed.DistributedMultiFieldFESpace} - const BlockFESpaceTypes{NB,SB,P} = Union{<:MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}},<:GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}} +using GridapDistributed: to_parray_of_arrays +using GridapDistributed: DistributedMultiFieldFESpace, DistributedMultiFieldFEFunction - include("BlockSolverInterfaces.jl") - include("BlockDiagonalSolvers.jl") - include("BlockTriangularSolvers.jl") +const MultiFieldFESpaceTypes = Union{<:MultiFieldFESpace,<:GridapDistributed.DistributedMultiFieldFESpace} +const BlockFESpaceTypes{NB,SB,P} = Union{<:MultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}},<:GridapDistributed.DistributedMultiFieldFESpace{<:BlockMultiFieldStyle{NB,SB,P}}} - include("BlockFEOperators.jl") - include("StaggeredFEOperators.jl") +include("BlockSolverInterfaces.jl") +include("BlockDiagonalSolvers.jl") +include("BlockTriangularSolvers.jl") - export MatrixBlock, LinearSystemBlock, NonlinearSystemBlock, BiformBlock, TriformBlock +include("BlockFEOperators.jl") +include("StaggeredFEOperators.jl") - export BlockDiagonalSolver - export BlockTriangularSolver +export MatrixBlock, LinearSystemBlock, NonlinearSystemBlock, BiformBlock, TriformBlock - export BlockFEOperator - export StaggeredFEOperator, StaggeredAffineFEOperator, StaggeredFESolver +export BlockDiagonalSolver +export BlockTriangularSolver + +export BlockFEOperator +export StaggeredFEOperator, StaggeredAffineFEOperator, StaggeredFESolver end diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl index d460d227..3d38bdfb 100644 --- a/src/BlockSolvers/StaggeredFEOperators.jl +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -7,15 +7,30 @@ end function get_solution(op::StaggeredFEOperator{NB,SB}, xh::MultiFieldFEFunction, k) where {NB,SB} r = MultiField.get_block_ranges(op)[k] - if length(r) == 1 + if isone(length(r)) # SingleField xh_k = xh[r[1]] - else - fv_k = blocks(xh.free_values)[k] + else # MultiField + fv_k = blocks(get_free_dof_values(xh))[k] xh_k = MultiFieldFEFunction(fv_k, op.trials[k], xh.single_fe_functions[r]) end return xh_k end +function get_solution(op::StaggeredFEOperator{NB,SB}, xh::DistributedMultiFieldFEFunction, k) where {NB,SB} + r = MultiField.get_block_ranges(op)[k] + if isone(length(r)) # SingleField + xh_k = xh[r[1]] + else # MultiField + sf_k = xh.field_fe_fun[r] + fv_k = blocks(get_free_dof_values(xh))[k] + mf_k = map(local_views(op.trials[k]),partition(fv_k),map(local_views,sf_k)...) do Vk, fv_k, sf_k... + MultiFieldFEFunction(fv_k, Vk, [sf_k...]) + end + xh_k = DistributedMultiFieldFEFunction(sf_k, mf_k, fv_k) + end + return xh_k +end + # StaggeredFESolver struct StaggeredFESolver{NB} <: FESpaces.FESolver @@ -33,6 +48,7 @@ function Algebra.solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperat xh_k = get_solution(op,xh,k) op_k = get_operator(op,xhs,k) xh_k, cache_k = solve!(xh_k,solvers[k],op_k,nothing) + copy!(get_free_dof_values(get_solution(op,xh,k)),get_free_dof_values(xh_k)) xhs, caches, operators = (xhs...,xh_k), (caches...,cache_k), (operators...,op_k) end return xh, (caches,operators) @@ -104,7 +120,7 @@ function combine_fespaces(spaces::Vector{<:FESpace}) end split_fespace(space::FESpace) = [space] -split_fespace(space::MultiFieldFESpaceTypes) = space.spaces +split_fespace(space::MultiFieldFESpaceTypes) = [space...] function get_operator(op::StaggeredFEOperator{NB}, xhs, k) where NB @assert NB >= k diff --git a/src/SolverInterfaces/GridapExtras.jl b/src/SolverInterfaces/GridapExtras.jl index 8d867c55..d8784b63 100644 --- a/src/SolverInterfaces/GridapExtras.jl +++ b/src/SolverInterfaces/GridapExtras.jl @@ -12,3 +12,47 @@ end function Gridap.Algebra.numerical_setup!(ns::Gridap.Algebra.NumericalSetup,A::AbstractMatrix,x::AbstractVector) numerical_setup!(ns,A) end + +# Similar to PartitionedArrays.matching_local_indices, but cheaper since +# we do not try to match the indices. + +function same_local_indices(a::PRange,b::PRange) + partition(a) === partition(b) && return true + c = map(===,partition(a),partition(b)) + reduce(&,c,init=true) +end + +function same_local_indices(a::BlockPRange,b::BlockPRange) + c = map(same_local_indices,blocks(a),blocks(b)) + reduce(&,c,init=true) +end + +# The following is needed, otherwise the input vector `x` potentially does not match +# the domain of the operator `op`. + +function Algebra.solve!(x::PVector,ls::LinearSolver,op::AffineOperator,cache::Nothing) + A = op.matrix + b = op.vector + ss = symbolic_setup(ls,A) + ns = numerical_setup(ss,A) + y = allocate_in_domain(A) + copy!(y,x) + solve!(y,ns,b) + copy!(x,y) + consistent!(x) |> wait + return ns, y +end + +function Algebra.solve!(x::PVector,ls::LinearSolver,op::AffineOperator,cache,newmatrix::Bool) + A = op.matrix + b = op.vector + ns, y = cache + if newmatrix + numerical_setup!(ns,A) + end + copy!(y,x) + solve!(y,ns,b) + copy!(x,y) + consistent!(x) |> wait + return cache +end diff --git a/src/SolverInterfaces/SolverInterfaces.jl b/src/SolverInterfaces/SolverInterfaces.jl index 8b9aa724..4bcb3d32 100644 --- a/src/SolverInterfaces/SolverInterfaces.jl +++ b/src/SolverInterfaces/SolverInterfaces.jl @@ -4,6 +4,10 @@ using Gridap using Gridap.Helpers using Gridap.Algebra +using PartitionedArrays +using GridapDistributed +using GridapDistributed: BlockPRange + using AbstractTrees using Printf diff --git a/test/BlockSolvers/StaggeredFEOperatorsTests.jl b/test/BlockSolvers/StaggeredFEOperatorsTests.jl index d3a7f826..093d2c1a 100644 --- a/test/BlockSolvers/StaggeredFEOperatorsTests.jl +++ b/test/BlockSolvers/StaggeredFEOperatorsTests.jl @@ -5,9 +5,12 @@ using Gridap, GridapDistributed, PartitionedArrays using Gridap.MultiField using BlockArrays using GridapSolvers -using GridapSolvers.BlockSolvers +using GridapSolvers.BlockSolvers, GridapSolvers.LinearSolvers -model = CartesianDiscreteModel((0,1,0,1),(4,4)) +np = (2,2) +ranks = DebugArray(LinearIndices((prod(np),))) + +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4)) order = 1 reffe = ReferenceFE(lagrangian,Float64,order) @@ -40,14 +43,13 @@ UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V]) UB3, VB3 = U4, V op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],[UB1,UB2,UB3],[VB1,VB2,VB3]) -solver = StaggeredFESolver([LUSolver(),LUSolver(),LUSolver()]) +solver = StaggeredFESolver(fill(CGSolver(JacobiLinearSolver();rtol=1.e-10),3)) xh = solve(solver,op) xh_exact = interpolate(sol,X) - for k in 1:4 eh_k = xh[k] - xh_exact[k] - e_k = sum(∫(eh_k * eh_k)dΩ) + e_k = sqrt(sum(∫(eh_k * eh_k)dΩ)) println("Error in field $k: $e_k") @test e_k < 1e-10 end From 4e805601401d24c958bb7fec516a7777421d55cb Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 11 Dec 2024 20:19:05 +1100 Subject: [PATCH 3/8] Added documentation --- docs/src/BlockSolvers.md | 10 ++++ src/BlockSolvers/StaggeredFEOperators.jl | 63 ++++++++++++++++++++++++ 2 files changed, 73 insertions(+) diff --git a/docs/src/BlockSolvers.md b/docs/src/BlockSolvers.md index 49a85d93..ed746f20 100644 --- a/docs/src/BlockSolvers.md +++ b/docs/src/BlockSolvers.md @@ -59,3 +59,13 @@ BlockTriangularSolver BlockTriangularSolver(blocks::AbstractMatrix{<:SolverBlock},solvers ::AbstractVector{<:LinearSolver},) BlockTriangularSolver(solvers::AbstractVector{<:LinearSolver}) ``` + +## Staggered FE Operators + +```@docs +StaggeredFESolver +StaggeredFEOperator +AffineStaggeredFESolver +NonlinearStaggeredFESolver +solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, ::Nothing) +``` diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl index 3d38bdfb..53558ff3 100644 --- a/src/BlockSolvers/StaggeredFEOperators.jl +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -1,4 +1,22 @@ +""" + abstract type StaggeredFEOperator{NB,SB} <: FESpaces.FEOperator end + +Staggered operator, used to solve staggered problems. + +We define a staggered problem as a multi-variable non-linear problem where the equation +for the k-th variable `u_k` only depends on the previous variables `u_1,...,u_{k-1}` (and itself). + +Such a problem can then be solved by solving each variable sequentially, +using the previous variables as input. The most common examples of staggered problems +are one-directional coupling problems, where the variables are coupled in a chain-like manner. + +Two types of staggered operators are currently supported: + +- [`StaggeredAffineFEOperator`](@ref): when the k-th equation is linear in `u_k`. +- [`StaggeredNonlinearFEOperator`](@ref): when the k-th equation is non-linear in `u_k`. + +""" abstract type StaggeredFEOperator{NB,SB} <: FESpaces.FEOperator end function MultiField.get_block_ranges(::StaggeredFEOperator{NB,SB}) where {NB,SB} @@ -33,6 +51,13 @@ end # StaggeredFESolver +""" + struct StaggeredFESolver{NB} <: FESpaces.FESolver + solvers :: Vector{<:Union{LinearSolver,NonlinearSolver}} + end + +Solver for staggered problems. See [`StaggeredFEOperator`](@ref) for more details. +""" struct StaggeredFESolver{NB} <: FESpaces.FESolver solvers :: Vector{<:NonlinearSolver} function StaggeredFESolver(solvers::Vector{<:NonlinearSolver}) @@ -41,6 +66,11 @@ struct StaggeredFESolver{NB} <: FESpaces.FESolver end end +""" + solve(solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}) + solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, cache::Nothing) where NB + solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, cache) where NB +""" function Algebra.solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, ::Nothing) where NB solvers = solver.solvers xhs, caches, operators = (), (), () @@ -69,6 +99,25 @@ end # StaggeredAffineFEOperator +""" + struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} + ... + end + +Affine staggered operator, used to solve staggered problems +where the k-th equation is linear in `u_k`. + +Such a problem is formulated by a set of bilinear/linear form pairs: + + a_k((u_1,...,u_{k-1}),u_k,v_k) = ∫(...) + l_k((u_1,...,u_{k-1}),v_k) = ∫(...) + +than cam be assembled into a set of linear systems: + + A_k u_k = b_k + +where `A_k` and `b_k` only depend on the previous variables `u_1,...,u_{k-1}`. +""" struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} biforms :: Vector{<:Function} liforms :: Vector{<:Function} @@ -145,6 +194,20 @@ end ############################################################################################ # StaggeredFEOperator +""" + struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} + ... + end + +Nonlinear staggered operator, used to solve staggered problems +where the k-th equation is nonlinear in `u_k`. + +Such a problem is formulated by a set of residual/jacobian pairs: + + jac_k((u_1,...,u_{k-1}),u_k,du_k,dv_k) = ∫(...) + res_k((u_1,...,u_{k-1}),u_k,v_k) = ∫(...) + +""" struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} jacobians :: Vector{<:Function} residuals :: Vector{<:Function} From 4af05f8eb2fcfa2e59b2b699b8c23a89a108da91 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 11 Dec 2024 20:41:46 +1100 Subject: [PATCH 4/8] More documentation --- docs/src/BlockSolvers.md | 8 +++---- docs/src/index.md | 29 ++++++++++++++++++++++++ src/BlockSolvers/StaggeredFEOperators.jl | 26 +++++++++++++++++++++ 3 files changed, 59 insertions(+), 4 deletions(-) diff --git a/docs/src/BlockSolvers.md b/docs/src/BlockSolvers.md index ed746f20..bd3598e8 100644 --- a/docs/src/BlockSolvers.md +++ b/docs/src/BlockSolvers.md @@ -60,12 +60,12 @@ BlockTriangularSolver(blocks::AbstractMatrix{<:SolverBlock},solvers ::AbstractVe BlockTriangularSolver(solvers::AbstractVector{<:LinearSolver}) ``` -## Staggered FE Operators +## Staggered FEOperators ```@docs StaggeredFESolver StaggeredFEOperator -AffineStaggeredFESolver -NonlinearStaggeredFESolver -solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, ::Nothing) +StaggeredAffineFEOperator +StaggeredNonlinearFEOperator +solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperator{NB}, ::Nothing) where {NB} ``` diff --git a/docs/src/index.md b/docs/src/index.md index d80b3d16..d79076ed 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -10,6 +10,8 @@ GridapSolvers provides algebraic and non-algebraic solvers for the Gridap ecosys Solvers follow a modular design, where most blocks can be combined to produce PDE-taylored solvers for a wide range of problems. +## Table of contents + ```@contents Pages = [ "SolverInterfaces.md", @@ -20,3 +22,30 @@ Pages = [ "PatchBasedSmoothers.md" ] ``` + +## Documentation and examples + +A (hopefully) comprehensive documentation is available [here](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 + +```julia +pkg> add GridapSolvers +pkg> build +``` + +Building is required to link the external artifacts (e.g., PETSc, p4est) to the Julia environment. Restarting Julia is required after building in order to make the changes take effect. + +### Using custom binaries + +The previous installations steps will setup GridapSolvers to work using Julia's pre-compiled artifacts for MPI, PETSc and p4est. However, you can also link local copies of these libraries. This might be very desirable in clusters, where hardware-specific libraries might be faster/more stable than the ones provided by Julia. To do so, follow the next steps: + +- [MPI.jl](https://juliaparallel.org/MPI.jl/stable/configuration/) +- [GridapPETSc.jl](https://github.com/gridap/GridapPETSc.jl) +- [GridapP4est.jl](https://github.com/gridap/GridapP4est.jl), and [P4est_wrapper.jl](https://github.com/gridap/p4est_wrapper.jl) diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl index 53558ff3..2d3158e7 100644 --- a/src/BlockSolvers/StaggeredFEOperators.jl +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -142,6 +142,19 @@ struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} end end +""" + function StaggeredAffineFEOperator( + biforms::Vector{<:Function}, + liforms::Vector{<:Function}, + trials::Vector{<:FESpace}, + tests::Vector{<:FESpace}, + [assems::Vector{<:Assembler}] + ) + +Constructor for a `StaggeredAffineFEOperator` operator, taking in each +equation as a pair of bilinear/linear forms and the corresponding trial/test spaces. +The trial/test spaces can be single or multi-field spaces. +""" function StaggeredAffineFEOperator( biforms::Vector{<:Function}, liforms::Vector{<:Function}, @@ -233,6 +246,19 @@ struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} end # TODO: Can be compute jacobians from residuals? + +""" + function StaggeredNonlinearFEOperator( + res::Vector{<:Function}, + jac::Vector{<:Function}, + trials::Vector{<:FESpace}, + tests::Vector{<:FESpace} + ) + +Constructor for a `StaggeredNonlinearFEOperator` operator, taking in each +equation as a pair of residual/jacobian forms and the corresponding trial/test spaces. +The trial/test spaces can be single or multi-field spaces. +""" function StaggeredNonlinearFEOperator( res::Vector{<:Function}, jac::Vector{<:Function}, From d3a2e0678165fdb6e741f1a9392612da709bbc8c Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 11 Dec 2024 22:09:59 +1100 Subject: [PATCH 5/8] Updated NEWS.md --- CHANGELOG.md | 10 ---------- NEWS.md | 18 ++++++++++++++++++ 2 files changed, 18 insertions(+), 10 deletions(-) delete mode 100644 CHANGELOG.md create mode 100644 NEWS.md diff --git a/CHANGELOG.md b/CHANGELOG.md deleted file mode 100644 index 5aad4212..00000000 --- a/CHANGELOG.md +++ /dev/null @@ -1,10 +0,0 @@ -# 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. diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 00000000..c21854d3 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,18 @@ +# 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). + +## [Unreleased] + +### Added + +- Added support for GMG in serial. Since PR[#68](https://github.com/gridap/GridapSolvers.jl/pull/68). +- Added Vanka-like smoothers in serial. Since PR[#68](https://github.com/gridap/GridapSolvers.jl/pull/68). +- Added `StaggeredFEOperators` and `StaggeredFESolvers`. Since PR[#84](https://github.com/gridap/GridapSolvers.jl/pull/84). + +## Previous versions + +A changelog is not maintained for older versions than 0.4.0. From 6062c1d10fe0c25031b0db2c7fa2c4be9a292b19 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 12 Dec 2024 00:36:58 +1100 Subject: [PATCH 6/8] Added StaggeredNonlinearFEOpearators --- src/BlockSolvers/BlockFEOperators.jl | 6 +- src/BlockSolvers/BlockSolvers.jl | 4 +- src/BlockSolvers/StaggeredFEOperators.jl | 196 +++++++++++------- .../BlockSolvers/StaggeredFEOperatorsTests.jl | 126 ++++++++--- 4 files changed, 227 insertions(+), 105 deletions(-) diff --git a/src/BlockSolvers/BlockFEOperators.jl b/src/BlockSolvers/BlockFEOperators.jl index 8c872baf..ff5c98eb 100644 --- a/src/BlockSolvers/BlockFEOperators.jl +++ b/src/BlockSolvers/BlockFEOperators.jl @@ -92,7 +92,7 @@ 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]) + (length(range) == 1) ? f[range[1]] : MultiFieldFESpace(f.spaces[range]) end return block_spaces end @@ -104,11 +104,11 @@ function BlockArrays.blocks(f::GridapDistributed.DistributedMultiFieldFESpace{<: if (length(range) == 1) space = f[range[1]] else - global_sf_spaces = f[range] + global_sf_spaces = f.field_fe_space[range] local_sf_spaces = 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) + space = DistributedMultiFieldFESpace(global_sf_spaces,local_mf_spaces,gids,vector_type) end space end diff --git a/src/BlockSolvers/BlockSolvers.jl b/src/BlockSolvers/BlockSolvers.jl index 17ced70f..23d29c89 100644 --- a/src/BlockSolvers/BlockSolvers.jl +++ b/src/BlockSolvers/BlockSolvers.jl @@ -14,6 +14,8 @@ using GridapDistributed using GridapSolvers.MultilevelTools using GridapSolvers.SolverInterfaces +using Gridap.MultiField: BlockSparseMatrixAssembler + using GridapDistributed: to_parray_of_arrays using GridapDistributed: DistributedMultiFieldFESpace, DistributedMultiFieldFEFunction @@ -33,6 +35,6 @@ export BlockDiagonalSolver export BlockTriangularSolver export BlockFEOperator -export StaggeredFEOperator, StaggeredAffineFEOperator, StaggeredFESolver +export StaggeredFEOperator, StaggeredAffineFEOperator, StaggeredNonlinearFEOperator, StaggeredFESolver end diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl index 2d3158e7..46765e4b 100644 --- a/src/BlockSolvers/StaggeredFEOperators.jl +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -23,6 +23,20 @@ function MultiField.get_block_ranges(::StaggeredFEOperator{NB,SB}) where {NB,SB} MultiField.get_block_ranges(NB,SB,Tuple(1:sum(SB))) end +# TODO: This is type piracy -> move to Gridap +MultiField.num_fields(space::FESpace) = 1 + +# TODO: We could reuse gids in distributed +function combine_fespaces(spaces::Vector{<:FESpace}) + NB = length(spaces) + SB = Tuple(map(num_fields,spaces)) + sf_spaces = vcat(map(split_fespace,spaces)...) + MultiFieldFESpace(sf_spaces; style = BlockMultiFieldStyle(NB,SB)) +end + +split_fespace(space::FESpace) = [space] +split_fespace(space::MultiFieldFESpaceTypes) = [space...] + function get_solution(op::StaggeredFEOperator{NB,SB}, xh::MultiFieldFEFunction, k) where {NB,SB} r = MultiField.get_block_ranges(op)[k] if isone(length(r)) # SingleField @@ -78,7 +92,6 @@ function Algebra.solve!(xh, solver::StaggeredFESolver{NB}, op::StaggeredFEOperat xh_k = get_solution(op,xh,k) op_k = get_operator(op,xhs,k) xh_k, cache_k = solve!(xh_k,solvers[k],op_k,nothing) - copy!(get_free_dof_values(get_solution(op,xh,k)),get_free_dof_values(xh_k)) xhs, caches, operators = (xhs...,xh_k), (caches...,cache_k), (operators...,op_k) end return xh, (caches,operators) @@ -127,12 +140,25 @@ struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} trial :: BlockFESpaceTypes{NB,SB} test :: BlockFESpaceTypes{NB,SB} + @doc """ + function StaggeredAffineFEOperator( + biforms :: Vector{<:Function}, + liforms :: Vector{<:Function}, + trials :: Vector{<:FESpace}, + tests :: Vector{<:FESpace}, + [assems :: Vector{<:Assembler}] + ) + + Constructor for a `StaggeredAffineFEOperator` operator, taking in each + equation as a pair of bilinear/linear forms and the corresponding trial/test spaces. + The trial/test spaces can be single or multi-field spaces. + """ function StaggeredAffineFEOperator( - biforms::Vector{<:Function}, - liforms::Vector{<:Function}, - trials::Vector{<:FESpace}, - tests::Vector{<:FESpace}, - assems::Vector{<:Assembler} + biforms :: Vector{<:Function}, + liforms :: Vector{<:Function}, + trials :: Vector{<:FESpace}, + tests :: Vector{<:FESpace}, + assems :: Vector{<:Assembler} = map(SparseMatrixAssembler,tests,trials) ) @assert length(biforms) == length(liforms) == length(trials) == length(tests) == length(assems) trial = combine_fespaces(trials) @@ -140,68 +166,53 @@ struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} NB, SB = length(trials), Tuple(map(num_fields,trials)) new{NB,SB}(biforms,liforms,trials,tests,assems,trial,test) end -end -""" + @doc """ function StaggeredAffineFEOperator( - biforms::Vector{<:Function}, - liforms::Vector{<:Function}, - trials::Vector{<:FESpace}, - tests::Vector{<:FESpace}, - [assems::Vector{<:Assembler}] - ) - -Constructor for a `StaggeredAffineFEOperator` operator, taking in each -equation as a pair of bilinear/linear forms and the corresponding trial/test spaces. -The trial/test spaces can be single or multi-field spaces. -""" -function StaggeredAffineFEOperator( - biforms::Vector{<:Function}, - liforms::Vector{<:Function}, - trials::Vector{<:FESpace}, - tests::Vector{<:FESpace} -) - assems = map(SparseMatrixAssembler,tests,trials) - return StaggeredAffineFEOperator(biforms,liforms,trials,tests,assems) + biforms :: Vector{<:Function}, + liforms :: Vector{<:Function}, + trial :: BlockFESpaceTypes{NB,SB,P}, + test :: BlockFESpaceTypes{NB,SB,P}, + [assem :: BlockSparseMatrixAssembler{NB,NV,SB,P}] + ) where {NB,NV,SB,P} + + Constructor for a `StaggeredAffineFEOperator` operator, taking in each + equation as a pair of bilinear/linear forms and the global trial/test spaces. + """ + function StaggeredAffineFEOperator( + biforms :: Vector{<:Function}, + liforms :: Vector{<:Function}, + trial :: BlockFESpaceTypes{NB,SB,P}, + test :: BlockFESpaceTypes{NB,SB,P}, + assem :: BlockSparseMatrixAssembler{NB,NV,SB,P} = SparseMatrixAssembler(trial,test) + ) where {NB,NV,SB,P} + @assert length(biforms) == length(liforms) == NB + @assert P == Tuple(1:sum(SB)) "Permutations not supported" + trials = blocks(trial) + tests = blocks(test) + assems = diag(blocks(assem)) + new{NB,SB}(biforms,liforms,trials,tests,assems,trial,test) + end end FESpaces.get_trial(op::StaggeredAffineFEOperator) = op.trial FESpaces.get_test(op::StaggeredAffineFEOperator) = op.test -# Utils - -MultiField.num_fields(space::FESpace) = 1 - -# TODO: We could reuse gids in distributed -function combine_fespaces(spaces::Vector{<:FESpace}) - NB = length(spaces) - SB = Tuple(map(num_fields,spaces)) - - sf_spaces = vcat(map(split_fespace,spaces)...) - MultiFieldFESpace(sf_spaces; style = BlockMultiFieldStyle(NB,SB)) -end - -split_fespace(space::FESpace) = [space] -split_fespace(space::MultiFieldFESpaceTypes) = [space...] - -function get_operator(op::StaggeredFEOperator{NB}, xhs, k) where NB +function get_operator(op::StaggeredAffineFEOperator{NB}, xhs, k) where NB @assert NB >= k a(uk,vk) = op.biforms[k](xhs,uk,vk) l(vk) = op.liforms[k](xhs,vk) - # uhd = zero(op.trials[k]) - # A, b = assemble_matrix_and_vector(a,l,op.assems[k],op.trials[k],op.tests[k],uhd) - # return AffineOperator(A,b) return AffineFEOperator(a,l,op.trials[k],op.tests[k],op.assems[k]) end -function get_operator!(op_k::AffineFEOperator, op::StaggeredFEOperator{NB}, xhs, k) where NB +function get_operator!(op_k::AffineFEOperator, op::StaggeredAffineFEOperator{NB}, xhs, k) where NB @assert NB >= k A, b = get_matrix(op_k), get_vector(op_k) a(uk,vk) = op.biforms[k](xhs,uk,vk) l(vk) = op.liforms[k](xhs,vk) uhd = zero(op.trials[k]) assemble_matrix_and_vector!(a,l,A,b,op.assems[k],op.trials[k],op.tests[k],uhd) - return AffineOperator(A,b) + return op_k end ############################################################################################ @@ -222,20 +233,32 @@ Such a problem is formulated by a set of residual/jacobian pairs: """ struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} - jacobians :: Vector{<:Function} residuals :: Vector{<:Function} + jacobians :: Vector{<:Function} trials :: Vector{<:FESpace} tests :: Vector{<:FESpace} assems :: Vector{<:Assembler} trial :: BlockFESpaceTypes{NB,SB} test :: BlockFESpaceTypes{NB,SB} + @doc """ + function StaggeredNonlinearFEOperator( + res :: Vector{<:Function}, + jac :: Vector{<:Function}, + trials :: Vector{<:FESpace}, + tests :: Vector{<:FESpace} + ) + + Constructor for a `StaggeredNonlinearFEOperator` operator, taking in each + equation as a pair of residual/jacobian forms and the corresponding trial/test spaces. + The trial/test spaces can be single or multi-field spaces. + """ function StaggeredNonlinearFEOperator( - res::Vector{<:Function}, - jac::Vector{<:Function}, - trials::Vector{<:FESpace}, - tests::Vector{<:FESpace}, - assems::Vector{<:Assembler} + res :: Vector{<:Function}, + jac :: Vector{<:Function}, + trials :: Vector{<:FESpace}, + tests :: Vector{<:FESpace}, + assems :: Vector{<:Assembler} = map(SparseMatrixAssembler,tests,trials) ) @assert length(res) == length(jac) == length(trials) == length(tests) == length(assems) trial = combine_fespaces(trials) @@ -243,31 +266,50 @@ struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} NB, SB = length(trials), Tuple(map(num_fields,trials)) new{NB,SB}(res,jac,trials,tests,assems,trial,test) end -end -# TODO: Can be compute jacobians from residuals? - -""" + @doc """ function StaggeredNonlinearFEOperator( - res::Vector{<:Function}, - jac::Vector{<:Function}, - trials::Vector{<:FESpace}, - tests::Vector{<:FESpace} - ) - -Constructor for a `StaggeredNonlinearFEOperator` operator, taking in each -equation as a pair of residual/jacobian forms and the corresponding trial/test spaces. -The trial/test spaces can be single or multi-field spaces. -""" -function StaggeredNonlinearFEOperator( - res::Vector{<:Function}, - jac::Vector{<:Function}, - trials::Vector{<:FESpace}, - tests::Vector{<:FESpace} -) - assems = map(SparseMatrixAssembler,tests,trials) - return StaggeredNonlinearFEOperator(res,jac,trials,tests,assems) + res :: Vector{<:Function}, + jac :: Vector{<:Function}, + trial :: BlockFESpaceTypes{NB,SB,P}, + test :: BlockFESpaceTypes{NB,SB,P}, + [assem :: BlockSparseMatrixAssembler{NB,NV,SB,P}] + ) where {NB,NV,SB,P} + + Constructor for a `StaggeredNonlinearFEOperator` operator, taking in each + equation as a pair of bilinear/linear forms and the global trial/test spaces. + """ + function StaggeredNonlinearFEOperator( + res :: Vector{<:Function}, + jac :: Vector{<:Function}, + trial :: BlockFESpaceTypes{NB,SB,P}, + test :: BlockFESpaceTypes{NB,SB,P}, + assem :: BlockSparseMatrixAssembler{NB,NV,SB,P} = SparseMatrixAssembler(trial,test) + ) where {NB,NV,SB,P} + @assert length(res) == length(jac) == NB + @assert P == Tuple(1:sum(SB)) "Permutations not supported" + trials = blocks(trial) + tests = blocks(test) + assems = diag(blocks(assem)) + new{NB,SB}(res,jac,trials,tests,assems,trial,test) + end end FESpaces.get_trial(op::StaggeredNonlinearFEOperator) = op.trial FESpaces.get_test(op::StaggeredNonlinearFEOperator) = op.test + +function get_operator(op::StaggeredNonlinearFEOperator{NB}, xhs, k) where NB + @assert NB >= k + jac(uk,duk,dvk) = op.jacobians[k](xhs,uk,duk,dvk) + res(uk,vk) = op.residuals[k](xhs,uk,vk) + return FESpaces.FEOperatorFromWeakForm(res,jac,op.trials[k],op.tests[k],op.assems[k]) +end + +function get_operator!( + op_k::FESpaces.FEOperatorFromWeakForm, op::StaggeredNonlinearFEOperator{NB}, xhs, k +) where NB + @assert NB >= k + jac(uk,duk,dvk) = op.jacobians[k](xhs,uk,duk,dvk) + res(uk,vk) = op.residuals[k](xhs,uk,vk) + return FESpaces.FEOperatorFromWeakForm(res,jac,op.trials[k],op.tests[k],op.assems[k]) +end diff --git a/test/BlockSolvers/StaggeredFEOperatorsTests.jl b/test/BlockSolvers/StaggeredFEOperatorsTests.jl index 093d2c1a..c4a8597d 100644 --- a/test/BlockSolvers/StaggeredFEOperatorsTests.jl +++ b/test/BlockSolvers/StaggeredFEOperatorsTests.jl @@ -5,13 +5,83 @@ using Gridap, GridapDistributed, PartitionedArrays using Gridap.MultiField using BlockArrays using GridapSolvers -using GridapSolvers.BlockSolvers, GridapSolvers.LinearSolvers +using GridapSolvers.BlockSolvers, GridapSolvers.LinearSolvers, GridapSolvers.NonlinearSolvers + +function test_solution(xh,sol,X,dΩ) + N = length(sol) + xh_exact = interpolate(sol,X) + for k in 1:N + eh_k = xh[k] - xh_exact[k] + e_k = sqrt(sum(∫(eh_k * eh_k)dΩ)) + println("Error in field $k: $e_k") + @test e_k < 1e-10 + end +end + +function driver_affine(model,verbose) + order = 1 + reffe = ReferenceFE(lagrangian,Float64,order) + V = FESpace(model,reffe;dirichlet_tags="boundary") + + sol = [x -> x[1], x -> x[2], x -> x[1] + x[2], x -> x[1] - x[2]] + U1 = TrialFESpace(V,sol[1]) + U2 = TrialFESpace(V,sol[2]) + U3 = TrialFESpace(V,sol[3]) + U4 = TrialFESpace(V,sol[4]) + + # Define weakforms + Ω = Triangulation(model) + dΩ = Measure(Ω,3*order) + + a1((),u1,v1) = ∫(u1 * v1)dΩ + l1((),v1) = ∫(sol[1] * v1)dΩ + + a2((u1,),(u2,u3),(v2,v3)) = ∫(u1 * u2 * v2)dΩ + ∫(u3 * v3)dΩ + l2((u1,),(v2,v3)) = ∫(sol[2] * u1 * v2)dΩ + ∫(sol[3] * v3)dΩ + + a3((u1,(u2,u3)),u4,v4) = ∫((u1 + u2) * u4 * v4)dΩ + l3((u1,(u2,u3)),v4) = ∫(sol[4] * (u1 + u2) * v4)dΩ + + # Define solver: Each block will be solved with a CG solver + lsolver = CGSolver(JacobiLinearSolver();rtol=1.e-10,verbose=verbose) + solver = StaggeredFESolver(fill(lsolver,3)) + + # Create operator from full spaces + mfs = BlockMultiFieldStyle(3,(1,2,1)) + X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs) + Y = MultiFieldFESpace([V,V,V,V];style=mfs) + op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],X,Y) + xh = solve(solver,op) + test_solution(xh,sol,X,dΩ) + + # Create operator from components + UB1, VB1 = U1, V + UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V]) + UB3, VB3 = U4, V + op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],[UB1,UB2,UB3],[VB1,VB2,VB3]) + + # Solve and keep caches for reuse + xh = zero(X) + xh, cache = solve!(xh,solver,op); + xh, cache = solve!(xh,solver,op,cache); + test_solution(xh,sol,X,dΩ) + + return true +end + + + +############################################################################################ np = (2,2) ranks = DebugArray(LinearIndices((prod(np),))) - +verbose = i_am_main(ranks) model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4)) +@testset "StaggeredAffineFEOperators" driver_affine(model,verbose) + +model = CartesianDiscreteModel((0,1,0,1),(4,4)) + order = 1 reffe = ReferenceFE(lagrangian,Float64,order) V = FESpace(model,reffe;dirichlet_tags="boundary") @@ -22,36 +92,44 @@ U2 = TrialFESpace(V,sol[2]) U3 = TrialFESpace(V,sol[3]) U4 = TrialFESpace(V,sol[4]) -mfs = BlockMultiFieldStyle(3,(1,2,1)) -X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs) -Y = MultiFieldFESpace([V,V,V,V];style=mfs) - +# Define weakforms Ω = Triangulation(model) -dΩ = Measure(Ω,3*order) +dΩ = Measure(Ω,4*order) + +F(u::Function) = x -> (u(x) + 1) * u(x) +F(u) = (u + 1) * u +dF(u,du) = 2.0 * u * du + du + +j1((),u1,du1,dv1) = ∫(dF(u1,du1) * dv1)dΩ +r1((),u1,v1) = ∫((F(u1) - F(sol[1])) * v1)dΩ -a1((),u1,v1) = ∫(u1 * v1)dΩ -l1((),v1) = ∫(sol[1] * v1)dΩ +j2((u1,),(u2,u3),(du2,du3),(dv2,dv3)) = ∫(u1 * dF(u2,du2) * dv2)dΩ + ∫(dF(u3,du3) * dv3)dΩ +r2((u1,),(u2,u3),(v2,v3)) = ∫(u1 * (F(u2) - F(sol[2])) * v2)dΩ + ∫((F(u3) - F(sol[3])) * v3)dΩ -a2((u1,),(u2,u3),(v2,v3)) = ∫(u1 * u2 * v2)dΩ + ∫(u3 * v3)dΩ -l2((u1,),(v2,v3)) = ∫(sol[2] * u1 * v2)dΩ + ∫(sol[3] * v3)dΩ +j3((u1,(u2,u3)),u4,du4,dv4) = ∫(u3 * dF(u4,du4) * dv4)dΩ +r3((u1,(u2,u3)),u4,v4) = ∫(u3 * (F(u4) - F(sol[4])) * v4)dΩ -a3((u1,(u2,u3)),u4,v4) = ∫((u1 + u2) * u4 * v4)dΩ -l3((u1,(u2,u3)),v4) = ∫(sol[4] * (u1 + u2) * v4)dΩ +# Define solver: Each block will be solved with a LU solver +lsolver = LUSolver() +nlsolver = NewtonSolver(lsolver;rtol=1.e-10,verbose=verbose) +solver = StaggeredFESolver(fill(nlsolver,3)) +# Create operator from full spaces +mfs = BlockMultiFieldStyle(3,(1,2,1)) +X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs) +Y = MultiFieldFESpace([V,V,V,V];style=mfs) +op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],X,Y) +xh = solve(solver,op) +test_solution(xh,sol,X,dΩ) + +# Create operator from components UB1, VB1 = U1, V UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V]) UB3, VB3 = U4, V -op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],[UB1,UB2,UB3],[VB1,VB2,VB3]) - -solver = StaggeredFESolver(fill(CGSolver(JacobiLinearSolver();rtol=1.e-10),3)) -xh = solve(solver,op) +op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],[UB1,UB2,UB3],[VB1,VB2,VB3]) -xh_exact = interpolate(sol,X) -for k in 1:4 - eh_k = xh[k] - xh_exact[k] - e_k = sqrt(sum(∫(eh_k * eh_k)dΩ)) - println("Error in field $k: $e_k") - @test e_k < 1e-10 -end +xh = zero(X) +xh, cache = solve!(xh,solver,op); +test_solution(xh,sol,X,dΩ) end # module \ No newline at end of file From 456387681131b93aeb3d63b866d4c1c2073e2758 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 12 Dec 2024 08:54:53 +1100 Subject: [PATCH 7/8] Activated tests --- src/BlockSolvers/StaggeredFEOperators.jl | 3 +- .../BlockSolvers/StaggeredFEOperatorsTests.jl | 136 ++++++++++-------- .../mpi/StaggeredFEOperatorsTests.jl | 9 ++ .../seq/StaggeredFEOperatorsTests.jl | 7 + test/BlockSolvers/seq/runtests.jl | 1 + 5 files changed, 91 insertions(+), 65 deletions(-) create mode 100644 test/BlockSolvers/mpi/StaggeredFEOperatorsTests.jl create mode 100644 test/BlockSolvers/seq/StaggeredFEOperatorsTests.jl diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl index 46765e4b..0ba7676e 100644 --- a/src/BlockSolvers/StaggeredFEOperators.jl +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -210,8 +210,7 @@ function get_operator!(op_k::AffineFEOperator, op::StaggeredAffineFEOperator{NB} A, b = get_matrix(op_k), get_vector(op_k) a(uk,vk) = op.biforms[k](xhs,uk,vk) l(vk) = op.liforms[k](xhs,vk) - uhd = zero(op.trials[k]) - assemble_matrix_and_vector!(a,l,A,b,op.assems[k],op.trials[k],op.tests[k],uhd) + assemble_matrix_and_vector!(a,l,A,b,op.assems[k],op.trials[k],op.tests[k],zero(op.trials[k])) return op_k end diff --git a/test/BlockSolvers/StaggeredFEOperatorsTests.jl b/test/BlockSolvers/StaggeredFEOperatorsTests.jl index c4a8597d..01771498 100644 --- a/test/BlockSolvers/StaggeredFEOperatorsTests.jl +++ b/test/BlockSolvers/StaggeredFEOperatorsTests.jl @@ -7,13 +7,13 @@ using BlockArrays using GridapSolvers using GridapSolvers.BlockSolvers, GridapSolvers.LinearSolvers, GridapSolvers.NonlinearSolvers -function test_solution(xh,sol,X,dΩ) +function test_solution(xh,sol,X,dΩ,verbose) N = length(sol) xh_exact = interpolate(sol,X) for k in 1:N eh_k = xh[k] - xh_exact[k] e_k = sqrt(sum(∫(eh_k * eh_k)dΩ)) - println("Error in field $k: $e_k") + verbose && println("Error in field $k: $e_k") @test e_k < 1e-10 end end @@ -43,7 +43,7 @@ function driver_affine(model,verbose) l3((u1,(u2,u3)),v4) = ∫(sol[4] * (u1 + u2) * v4)dΩ # Define solver: Each block will be solved with a CG solver - lsolver = CGSolver(JacobiLinearSolver();rtol=1.e-10,verbose=verbose) + lsolver = CGSolver(JacobiLinearSolver();rtol=1.e-12,verbose=verbose) solver = StaggeredFESolver(fill(lsolver,3)) # Create operator from full spaces @@ -52,7 +52,7 @@ function driver_affine(model,verbose) Y = MultiFieldFESpace([V,V,V,V];style=mfs) op = StaggeredAffineFEOperator([a1,a2,a3],[l1,l2,l3],X,Y) xh = solve(solver,op) - test_solution(xh,sol,X,dΩ) + test_solution(xh,sol,X,dΩ,verbose) # Create operator from components UB1, VB1 = U1, V @@ -64,72 +64,82 @@ function driver_affine(model,verbose) xh = zero(X) xh, cache = solve!(xh,solver,op); xh, cache = solve!(xh,solver,op,cache); - test_solution(xh,sol,X,dΩ) + test_solution(xh,sol,X,dΩ,verbose) return true end +function test_nonlinear(model,verbose) + order = 1 + reffe = ReferenceFE(lagrangian,Float64,order) + V = FESpace(model,reffe;dirichlet_tags="boundary") + + sol = [x -> x[1], x -> x[2], x -> x[1] + x[2], x -> 2.0*x[1]] + U1 = TrialFESpace(V,sol[1]) + U2 = TrialFESpace(V,sol[2]) + U3 = TrialFESpace(V,sol[3]) + U4 = TrialFESpace(V,sol[4]) + + # Define weakforms + Ω = Triangulation(model) + dΩ = Measure(Ω,4*order) + + F(u::Function) = x -> (u(x) + 1) * u(x) + F(u) = (u + 1) * u + dF(u,du) = 2.0 * u * du + du + + j1((),u1,du1,dv1) = ∫(dF(u1,du1) * dv1)dΩ + r1((),u1,v1) = ∫((F(u1) - F(sol[1])) * v1)dΩ + + j2((u1,),(u2,u3),(du2,du3),(dv2,dv3)) = ∫(u1 * dF(u2,du2) * dv2)dΩ + ∫(dF(u3,du3) * dv3)dΩ + r2((u1,),(u2,u3),(v2,v3)) = ∫(u1 * (F(u2) - F(sol[2])) * v2)dΩ + ∫((F(u3) - F(sol[3])) * v3)dΩ + + j3((u1,(u2,u3)),u4,du4,dv4) = ∫(u3 * dF(u4,du4) * dv4)dΩ + r3((u1,(u2,u3)),u4,v4) = ∫(u3 * (F(u4) - F(sol[4])) * v4)dΩ + + # Define solver: Each block will be solved with a LU solver + lsolver = LUSolver() + nlsolver = NewtonSolver(lsolver;rtol=1.e-10,verbose=verbose) + solver = StaggeredFESolver(fill(nlsolver,3)) + + # Create operator from full spaces + mfs = BlockMultiFieldStyle(3,(1,2,1)) + X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs) + Y = MultiFieldFESpace([V,V,V,V];style=mfs) + op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],X,Y) + xh = zero(X) + xh, cache = solve!(xh,solver,op); + test_solution(xh,sol,X,dΩ,verbose) + + # Create operator from components + UB1, VB1 = U1, V + UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V]) + UB3, VB3 = U4, V + op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],[UB1,UB2,UB3],[VB1,VB2,VB3]) + + xh = zero(X) + xh, cache = solve!(xh,solver,op,cache); + test_solution(xh,sol,X,dΩ,verbose) + return true +end -############################################################################################ - -np = (2,2) -ranks = DebugArray(LinearIndices((prod(np),))) -verbose = i_am_main(ranks) -model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4)) - -@testset "StaggeredAffineFEOperators" driver_affine(model,verbose) - -model = CartesianDiscreteModel((0,1,0,1),(4,4)) - -order = 1 -reffe = ReferenceFE(lagrangian,Float64,order) -V = FESpace(model,reffe;dirichlet_tags="boundary") - -sol = [x -> x[1], x -> x[2], x -> x[1] + x[2], x -> x[1] - x[2]] -U1 = TrialFESpace(V,sol[1]) -U2 = TrialFESpace(V,sol[2]) -U3 = TrialFESpace(V,sol[3]) -U4 = TrialFESpace(V,sol[4]) - -# Define weakforms -Ω = Triangulation(model) -dΩ = Measure(Ω,4*order) - -F(u::Function) = x -> (u(x) + 1) * u(x) -F(u) = (u + 1) * u -dF(u,du) = 2.0 * u * du + du - -j1((),u1,du1,dv1) = ∫(dF(u1,du1) * dv1)dΩ -r1((),u1,v1) = ∫((F(u1) - F(sol[1])) * v1)dΩ - -j2((u1,),(u2,u3),(du2,du3),(dv2,dv3)) = ∫(u1 * dF(u2,du2) * dv2)dΩ + ∫(dF(u3,du3) * dv3)dΩ -r2((u1,),(u2,u3),(v2,v3)) = ∫(u1 * (F(u2) - F(sol[2])) * v2)dΩ + ∫((F(u3) - F(sol[3])) * v3)dΩ - -j3((u1,(u2,u3)),u4,du4,dv4) = ∫(u3 * dF(u4,du4) * dv4)dΩ -r3((u1,(u2,u3)),u4,v4) = ∫(u3 * (F(u4) - F(sol[4])) * v4)dΩ - -# Define solver: Each block will be solved with a LU solver -lsolver = LUSolver() -nlsolver = NewtonSolver(lsolver;rtol=1.e-10,verbose=verbose) -solver = StaggeredFESolver(fill(nlsolver,3)) - -# Create operator from full spaces -mfs = BlockMultiFieldStyle(3,(1,2,1)) -X = MultiFieldFESpace([U1,U2,U3,U4];style=mfs) -Y = MultiFieldFESpace([V,V,V,V];style=mfs) -op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],X,Y) -xh = solve(solver,op) -test_solution(xh,sol,X,dΩ) +function driver(model,verbose) + @testset "StaggeredAffineFEOperators" driver_affine(model,verbose) + @testset "StaggeredNonlinearFEOperators" driver_affine(model,verbose) +end -# Create operator from components -UB1, VB1 = U1, V -UB2, VB2 = MultiFieldFESpace([U2,U3]), MultiFieldFESpace([V,V]) -UB3, VB3 = U4, V -op = StaggeredNonlinearFEOperator([r1,r2,r3],[j1,j2,j3],[UB1,UB2,UB3],[VB1,VB2,VB3]) +# Distributed +function main(distribute,parts) + ranks = distribute(LinearIndices((prod(parts),))) + model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(8,8)) + driver(model,false) +end -xh = zero(X) -xh, cache = solve!(xh,solver,op); -test_solution(xh,sol,X,dΩ) +# Serial +function main() + model = CartesianDiscreteModel((0,1,0,1),(8,8)) + driver(model,false) +end end # module \ No newline at end of file diff --git a/test/BlockSolvers/mpi/StaggeredFEOperatorsTests.jl b/test/BlockSolvers/mpi/StaggeredFEOperatorsTests.jl new file mode 100644 index 00000000..97ac0bce --- /dev/null +++ b/test/BlockSolvers/mpi/StaggeredFEOperatorsTests.jl @@ -0,0 +1,9 @@ +module StaggeredFEOperatorsMPITests +using MPI, PartitionedArrays +include("../StaggeredFEOperatorsTests.jl") + +with_mpi() do distribute + StaggeredFEOperatorsTests.main(distribute,(2,2)) +end + +end \ No newline at end of file diff --git a/test/BlockSolvers/seq/StaggeredFEOperatorsTests.jl b/test/BlockSolvers/seq/StaggeredFEOperatorsTests.jl new file mode 100644 index 00000000..df49be90 --- /dev/null +++ b/test/BlockSolvers/seq/StaggeredFEOperatorsTests.jl @@ -0,0 +1,7 @@ +module StaggeredFEOperatorsSequentialTests +using PartitionedArrays +include("../StaggeredFEOperatorsTests.jl") + +StaggeredFEOperatorsTests.main() + +end \ No newline at end of file diff --git a/test/BlockSolvers/seq/runtests.jl b/test/BlockSolvers/seq/runtests.jl index 593d9768..30bb9c78 100644 --- a/test/BlockSolvers/seq/runtests.jl +++ b/test/BlockSolvers/seq/runtests.jl @@ -2,3 +2,4 @@ using Test @testset "BlockDiagonalSolvers" begin include("BlockDiagonalSolversTests.jl") end @testset "BlockTriangularSolvers" begin include("BlockTriangularSolversTests.jl") end +@testset "StaggeredFEOperators" begin include("StaggeredFEOperatorsTests.jl") end From cc769ac4d79c892f2b06499ccb9122f95b4ef3e2 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Sat, 14 Dec 2024 09:50:36 +1100 Subject: [PATCH 8/8] Minor --- src/BlockSolvers/StaggeredFEOperators.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/BlockSolvers/StaggeredFEOperators.jl b/src/BlockSolvers/StaggeredFEOperators.jl index 0ba7676e..c1f38510 100644 --- a/src/BlockSolvers/StaggeredFEOperators.jl +++ b/src/BlockSolvers/StaggeredFEOperators.jl @@ -158,7 +158,7 @@ struct StaggeredAffineFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} liforms :: Vector{<:Function}, trials :: Vector{<:FESpace}, tests :: Vector{<:FESpace}, - assems :: Vector{<:Assembler} = map(SparseMatrixAssembler,tests,trials) + assems :: Vector{<:Assembler} = map(SparseMatrixAssembler,trials,tests) ) @assert length(biforms) == length(liforms) == length(trials) == length(tests) == length(assems) trial = combine_fespaces(trials) @@ -257,7 +257,7 @@ struct StaggeredNonlinearFEOperator{NB,SB} <: StaggeredFEOperator{NB,SB} jac :: Vector{<:Function}, trials :: Vector{<:FESpace}, tests :: Vector{<:FESpace}, - assems :: Vector{<:Assembler} = map(SparseMatrixAssembler,tests,trials) + assems :: Vector{<:Assembler} = map(SparseMatrixAssembler,trials,tests) ) @assert length(res) == length(jac) == length(trials) == length(tests) == length(assems) trial = combine_fespaces(trials)