From b33cddb62549d3dee211978276980e5ea41e1e24 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 30 Jul 2025 09:58:30 +1000 Subject: [PATCH 1/7] implement split multifield grad ops --- src/Autodiff.jl | 28 ++++++++++++++ src/MultiField.jl | 14 ++++--- test.jl | 93 +++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 129 insertions(+), 6 deletions(-) create mode 100644 test.jl diff --git a/src/Autodiff.jl b/src/Autodiff.jl index 9651648..8d3a55d 100644 --- a/src/Autodiff.jl +++ b/src/Autodiff.jl @@ -80,6 +80,34 @@ function FESpaces._change_argument(op,f,local_trians,uh::DistributedADTypes) g end +# Distributed counterpart of: src/MultiField/MultiFieldAutodiff.jl + +for (op,_op) in ((:gradient,:_gradient),(:jacobian,:_jacobian)) + @eval begin + function FESpaces.$(op)(f::Function,uh::DistributedMultiFieldFEFunction;ad_type=:split) + fuh = f(uh) + if ad_type == :split + MultiField.multifield_autodiff_split($op,f,uh,fuh) + elseif ad_type == :monolithic + FESpaces.$(_op)(f,uh,fuh) + else + @notimplemented """Unknown ad_type = $ad_type + Options: + - :split -- compute the gradient for each field separately, then merge + - :monolithic -- compute the gradient for all fields together + """ + end + end + end +end + +function MultiField._combine_contributions(op::Function,terms,fuh::DistributedDomainContribution) + local_terms = map(local_views(fuh),local_views.(terms)...) do fuh,terms... + MultiField._combine_contributions(op,terms,fuh) + end + DistributedDomainContribution(local_terms) +end + # Distributed counterpart of: src/Arrays/Autodiff.jl # autodiff_array_xxx diff --git a/src/MultiField.jl b/src/MultiField.jl index c20a2a5..510b667 100644 --- a/src/MultiField.jl +++ b/src/MultiField.jl @@ -33,6 +33,8 @@ MultiField.num_fields(m::DistributedMultiFieldCellField) = length(m.field_fe_fun Base.iterate(m::DistributedMultiFieldCellField) = iterate(m.field_fe_fun) Base.iterate(m::DistributedMultiFieldCellField,state) = iterate(m.field_fe_fun,state) Base.getindex(m::DistributedMultiFieldCellField,field_id::Integer) = m.field_fe_fun[field_id] +Base.getindex(m::DistributedMultiFieldCellField,field_id::AbstractUnitRange) = m.field_fe_fun[field_id] +Base.lastindex(m::DistributedMultiFieldCellField) = num_fields(m) Base.length(m::DistributedMultiFieldCellField) = num_fields(m) function LinearAlgebra.dot(a::DistributedMultiFieldCellField,b::DistributedMultiFieldCellField) @@ -161,7 +163,7 @@ function FESpaces.EvaluationFunction( isconsistent=false ) free_values = change_ghost(_free_values,f.gids;is_consistent=isconsistent,make_consistent=true) - + # Create distributed single field functions field_fe_fun = DistributedSingleFieldFEFunction[] for i in 1:num_fields(f) @@ -221,7 +223,7 @@ function FESpaces.interpolate_everywhere!( ) msg = "free_values and FESpace have incompatible index partitions." @check PartitionedArrays.matching_local_indices(axes(free_values,1),get_free_dof_ids(space)) msg - + # Interpolate each field field_fe_fun = DistributedSingleFieldFEFunction[] for i in 1:num_fields(space) @@ -396,16 +398,16 @@ function generate_multi_field_gids( end collect(keys(dict)) end - + f_p_parts_snd, f_p_parts_rcv = map(x->assembly_neighbors(partition(x)),f_frange) |> tuple_of_arrays p_f_parts_snd = map(v,f_p_parts_snd...) p_f_parts_rcv = map(v,f_p_parts_rcv...) p_neigs_snd = map(merge_neigs,p_f_parts_snd) p_neigs_rcv = map(merge_neigs,p_f_parts_rcv) - + exchange_graph = ExchangeGraph(p_neigs_snd,p_neigs_rcv) assembly_neighbors(p_iset;neighbors=exchange_graph) - + PRange(p_iset) end @@ -446,7 +448,7 @@ end # BlockSparseMatrixAssemblers -const DistributedBlockSparseMatrixAssembler{R,C} = +const DistributedBlockSparseMatrixAssembler{R,C} = MultiField.BlockSparseMatrixAssembler{R,C,<:AbstractMatrix{<:DistributedSparseMatrixAssembler}} function FESpaces.SparseMatrixAssembler( diff --git a/test.jl b/test.jl new file mode 100644 index 0000000..a971327 --- /dev/null +++ b/test.jl @@ -0,0 +1,93 @@ +using Pkg; Pkg.activate("./GridapDistributed"); + +using Test + +using Gridap +using Gridap.Algebra +using Gridap.Arrays +using Gridap.Geometry +using Gridap.CellData +using Gridap.ReferenceFEs +using Gridap.FESpaces +using Gridap.MultiField + +using GridapDistributed +using PartitionedArrays + +function half_empty_trian(ranks,model) + cell_ids = get_cell_gids(model) + trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids + cell_mask = zeros(Bool, num_cells(model)) + if rank ∈ (3,4) + cell_mask[own_to_local(ids)] .= true + end + Triangulation(model,cell_mask) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + +np = (2,2) +ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) +end + +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(5,5)) +Ω = Triangulation(model) +# Γ = BoundaryTriangulation(model) +Γ = half_empty_trian(ranks,model) +dΩ = Measure(Ω,2) +dΓ = Measure(Γ,2) +V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1)) +V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1)) +V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1)) +X = MultiFieldFESpace([V1,V2,V3]) +uh = zero(X); + +f(xh) = ∫(xh[1]+xh[2]⋅xh[2]+xh[1]*xh[3])dΓ +df(v,xh) = ∫(v[1]+2*v[2]⋅xh[2]+v[1]*xh[3]+xh[1]*v[3])dΓ +du = gradient(f,uh) +du_vec = assemble_vector(du,X) +df_vec = assemble_vector(v->df(v,uh),X) + +@test df_vec ≈ du_vec + +f2(xh,yh) = ∫(xh[1]⋅yh[1]+xh[2]⋅yh[2]+xh[1]⋅xh[2]⋅yh[2]+xh[1]*xh[3]*yh[3])dΓ +dv = get_fe_basis(X); +j = jacobian(uh->f2(uh,dv),uh) +J = assemble_matrix(j,X,X) + +f2_jac(xh,dxh,yh) = ∫(dxh[1]⋅yh[1]+dxh[2]⋅yh[2]+dxh[1]⋅xh[2]⋅yh[2]+xh[1]⋅dxh[2]⋅yh[2]+dxh[1]*xh[3]*yh[3]+xh[1]*dxh[3]*yh[3])dΓ +op = FEOperator(f2,f2_jac,X,X) +J_fwd = jacobian(op,uh) + +@test reduce(&,map(≈,partition(J),partition(J_fwd))) + +# Skel +V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) +V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) +V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1),conformity=:L2) +X = MultiFieldFESpace([V1,V2,V3]) +uh = zero(X); +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ,2) + +f(xh) = ∫(mean(xh[1])+mean(xh[2])⋅mean(xh[2])+mean(xh[1])*mean(xh[3]))dΛ +df(v,xh) = ∫(mean(v[1])+2*mean(v[2])⋅mean(xh[2])+mean(v[1])*mean(xh[3])+mean(xh[1])*mean(v[3]))dΛ +du = gradient(f,uh) +du_vec = assemble_vector(du,X) +df_vec = assemble_vector(v->df(v,uh),X) + +@test df_vec ≈ du_vec + +# Skel jac +f2(xh,yh) = ∫(mean(xh[1])⋅mean(yh[1])+mean(xh[2])⋅mean(yh[2])+mean(xh[1])⋅mean(xh[2])⋅mean(yh[2])+mean(xh[1])*mean(xh[3])*mean(yh[3]))dΛ +dv = get_fe_basis(X); +j = jacobian(uh->f2(uh,dv),uh); +J = assemble_matrix(j,X,X) + +f2_jac(xh,dxh,yh) = ∫(mean(dxh[1])⋅mean(yh[1])+mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])⋅mean(xh[2])⋅mean(yh[2]) + + mean(xh[1])⋅mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])*mean(xh[3])*mean(yh[3])+mean(xh[1])*mean(dxh[3])*mean(yh[3]))dΛ +op = FEOperator(f2,f2_jac,X,X) +J_fwd = jacobian(op,uh) + +@test reduce(&,map(≈,partition(J),partition(J_fwd))) \ No newline at end of file From 721721912f90ea8469b1203eef16b1e96d8b3e72 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Wed, 6 Aug 2025 12:16:02 +1000 Subject: [PATCH 2/7] Save work --- src/Geometry.jl | 118 ++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 104 insertions(+), 14 deletions(-) diff --git a/src/Geometry.jl b/src/Geometry.jl index 3e25f38..8bc57ff 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -648,6 +648,8 @@ function filter_cells_when_needed( remove_ghost_cells(trian,cell_gids) end +# Removing ghost cells + function remove_ghost_cells(trian::DistributedTriangulation,gids) trians = map(remove_ghost_cells,local_views(trian),partition(gids)) model = get_background_model(trian) @@ -671,10 +673,8 @@ function remove_ghost_cells( remove_ghost_cells(glue,trian,gids) end -function remove_ghost_cells( - trian::AdaptedTriangulation{Dc,Dp,<:Union{SkeletonTriangulation,BoundaryTriangulation}},gids -) where {Dc,Dp} - remove_ghost_cells(trian.trian,gids) +function remove_ghost_cells(trian::AdaptedTriangulation,gids) + AdaptedTriangulation(remove_ghost_cells(trian.trian,gids),trian.adapted_model) end function remove_ghost_cells(glue::FaceToFaceGlue,trian,gids) @@ -720,24 +720,20 @@ function _find_owned_skeleton_facets(glue,gids) findall(part->part==part_id(gids),tface_to_part) end +# Adding ghost cells + function add_ghost_cells(dtrian::DistributedTriangulation) dmodel = get_background_model(dtrian) add_ghost_cells(dmodel,dtrian) end -function _covers_all_faces( - dmodel::DistributedDiscreteModel{Dm}, - dtrian::DistributedTriangulation{Dt} -) where {Dm,Dt} - covers_all_faces = map(local_views(dmodel),local_views(dtrian)) do model, trian - glue = get_glue(trian,Val(Dt)) - @assert isa(glue,FaceToFaceGlue) - isa(glue.tface_to_mface,IdentityVector) - end - reduce(&,covers_all_faces,init=true) +function add_ghost_cells(dmodel::DistributedDiscreteModel,dtrian::DistributedTriangulation) + T = eltype(local_views(dtrian)) + add_ghost_cells(T,dmodel,dtrian) end function add_ghost_cells( + ::Type{<:Triangulation}, dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} ) where {Dm,Dt} @@ -776,6 +772,100 @@ function add_ghost_cells( return DistributedTriangulation(trians,dmodel) end +function add_ghost_cells( + ::Type{<:BoundaryTriangulation}, + dmodel::DistributedDiscreteModel, + dtrian::DistributedTriangulation{Dt} +) where Dt + println("Flag 0") + bf_dtrian = add_ghost_cells(Triangulation, dmodel, dtrian) + (bf_trian === dtrian) && return dtrian # Case 1 above + + # Otherwise, we have created the underlying BodyFittedTriangulation but still have to + # create the expanded FaceToCellGlue... + + face_gids = get_face_gids(dmodel, Dt) + bgface_to_lcell = map(local_views(dtrian), partition(face_gids)) do trian, face_gids + glue = trian.glue + bgface_to_lcell = fill(Int8(1), length(face_gids)) + bgface_to_lcell[glue.face_to_bgface] .= glue.face_to_lcell + end + consistent!(PVector(bgface_to_lcell, partition(face_gids))) |> wait + + trians = map(local_views(bf_dtrian), local_views(dmodel), bgface_to_lcell) do bf_trian, model, bgface_to_lcell + topo = get_grid_topology(model) + face_grid = get_grid(bf_trian) + cell_grid = get_grid(model) + face_to_bgface = bf_trian.tface_to_mface + glue = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell) + BoundaryTriangulation(bf_trian,glue) + end + + return DistributedTriangulation(trians,dmodel) +end + +function add_ghost_cells( + ::Type{<:SkeletonTriangulation}, + dmodel::DistributedDiscreteModel, + dtrian::DistributedTriangulation{Dt} +) where Dt + bf_dtrian = add_ghost_cells(Triangulation, dmodel, dtrian) + (bf_trian === dtrian) && return dtrian # Case 1 above + + # Otherwise, we have created the underlying BodyFittedTriangulation but still have to + # create the expanded FaceToCellGlue... + + face_gids = get_face_gids(dmodel, Dt) + bgface_to_lcell_plus = map(local_views(dtrian), partition(face_gids)) do trian, face_gids + glue = trian.plus.glue + bgface_to_lcell = fill(Int8(1), length(face_gids)) + bgface_to_lcell[glue.face_to_bgface] .= glue.face_to_lcell + end + t_plus = consistent!(PVector(bgface_to_lcell, partition(face_gids))) + bgface_to_lcell_minus = map(local_views(dtrian), partition(face_gids)) do trian, face_gids + glue = trian.minus.glue + bgface_to_lcell = fill(Int8(1), length(face_gids)) + bgface_to_lcell[glue.face_to_bgface] .= glue.face_to_lcell + end + t_minus = consistent!(PVector(bgface_to_lcell, partition(face_gids))) + + wait(t_plus) + wait(t_minus) + + trians = map( + local_views(bf_dtrian), local_views(dmodel), bgface_to_lcell_plus, bgface_to_lcell_minus + ) do bf_trian, model, bgface_to_lcell_plus, bgface_to_lcell_minus + topo = get_grid_topology(model) + face_grid = get_grid(bf_trian) + cell_grid = get_grid(model) + face_to_bgface = bf_trian.tface_to_mface + + glue_plus = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell_plus) + plus = BoundaryTriangulation(bf_trian,glue_plus) + + glue_minus = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell_minus) + minus = BoundaryTriangulation(bf_trian,glue_minus) + + SkeletonTriangulation(plus,minus) + end + + return DistributedTriangulation(trians,dmodel) +end + +function _covers_all_faces( + dmodel::DistributedDiscreteModel{Dm}, + dtrian::DistributedTriangulation{Dt} +) where {Dm,Dt} + covers_all_faces = map(local_views(dmodel),local_views(dtrian)) do model, trian + glue = get_glue(trian,Val(Dt)) + @assert isa(glue,FaceToFaceGlue) + isa(glue.tface_to_mface,IdentityVector) + end + reduce(&,covers_all_faces,init=true) +end + +# Triangulation gids + function generate_cell_gids(dtrian::DistributedTriangulation) dmodel = get_background_model(dtrian) generate_cell_gids(dmodel,dtrian) From 6cfeb873646513943f757e9ccfa76ce76aa4ad55 Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Sat, 16 Aug 2025 11:41:59 +1000 Subject: [PATCH 3/7] Changes to fix issue 177 --- src/Geometry.jl | 234 +++++++++++++++--------------------------------- 1 file changed, 71 insertions(+), 163 deletions(-) diff --git a/src/Geometry.jl b/src/Geometry.jl index 8bc57ff..d04520f 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -451,15 +451,19 @@ end # This object cannot implement the Triangulation interface in a strict sense """ """ -struct DistributedTriangulation{Dc,Dp,A,B} <: GridapType - trians::A - model::B +struct DistributedTriangulation{Dc,Dp,A,B,C} <: GridapType + trians ::A + model ::B + metadata::C function DistributedTriangulation( trians::AbstractArray{<:Triangulation{Dc,Dp}}, - model::DistributedDiscreteModel) where {Dc,Dp} + model::DistributedDiscreteModel; + metadata = nothing + ) where {Dc,Dp} A = typeof(trians) B = typeof(model) - new{Dc,Dp,A,B}(trians,model) + C = typeof(metadata) + new{Dc,Dp,A,B,C}(trians,model,metadata) end end @@ -494,58 +498,65 @@ end # Triangulation constructors -function Geometry.Triangulation( - model::DistributedDiscreteModel;kwargs...) - D=num_cell_dims(model) +function Geometry.Triangulation(model::DistributedDiscreteModel;kwargs...) + D = num_cell_dims(model) Triangulation(no_ghost,ReferenceFE{D},model;kwargs...) end -function Geometry.BoundaryTriangulation( - model::DistributedDiscreteModel;kwargs...) +function Geometry.Triangulation(::Type{ReferenceFE{D}}, model::DistributedDiscreteModel;kwargs...) where D + Triangulation(no_ghost, ReferenceFE{D}, model; kwargs...) +end + +function Geometry.BoundaryTriangulation(model::DistributedDiscreteModel;kwargs...) BoundaryTriangulation(no_ghost,model;kwargs...) end -function Geometry.BoundaryTriangulation( - trian::DistributedTriangulation;kwargs...) +function Geometry.BoundaryTriangulation(trian::DistributedTriangulation;kwargs...) BoundaryTriangulation(no_ghost,trian;kwargs...) end -function Geometry.SkeletonTriangulation( - model::DistributedDiscreteModel;kwargs...) +function Geometry.SkeletonTriangulation(model::DistributedDiscreteModel;kwargs...) SkeletonTriangulation(no_ghost,model;kwargs...) end -function Geometry.SkeletonTriangulation( - trian::DistributedTriangulation;kwargs...) +function Geometry.SkeletonTriangulation(trian::DistributedTriangulation;kwargs...) SkeletonTriangulation(no_ghost,trian;kwargs...) end +function Geometry.Triangulation(portion, model::DistributedDiscreteModel;kwargs...) + D = num_cell_dims(model) + Triangulation(portion,ReferenceFE{D},model;kwargs...) +end + function Geometry.Triangulation( portion,::Type{ReferenceFE{Dt}},model::DistributedDiscreteModel{Dm};kwargs...) where {Dt,Dm} # Generate global ordering for the faces of dimension Dt (if needed) gids = get_face_gids(model,Dt) - trians = map(local_views(model),partition(gids)) do model, gids - Triangulation(portion,gids,ReferenceFE{Dt},model;kwargs...) + trians = map(local_views(model)) do model + Triangulation(ReferenceFE{Dt},model;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end function Geometry.BoundaryTriangulation( portion,model::DistributedDiscreteModel{Dc};kwargs...) where Dc gids = get_face_gids(model,Dc) - trians = map(local_views(model),partition(gids)) do model, gids - BoundaryTriangulation(portion,gids,model;kwargs...) + trians = map(local_views(model)) do model + BoundaryTriangulation(model;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end function Geometry.SkeletonTriangulation( portion,model::DistributedDiscreteModel{Dc};kwargs...) where Dc gids = get_face_gids(model,Dc) - trians = map(local_views(model),partition(gids)) do model, gids - SkeletonTriangulation(portion,gids,model;kwargs...) + trians = map(local_views(model)) do model + SkeletonTriangulation(model;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end # NOTE: The following constructors require adding back the ghost cells: @@ -557,10 +568,11 @@ function Geometry.BoundaryTriangulation( model = get_background_model(trian) gids = get_cell_gids(model) ghosted_trian = add_ghost_cells(trian) - trians = map(local_views(ghosted_trian),partition(gids)) do trian, gids - BoundaryTriangulation(portion,gids,trian;kwargs...) + trians = map(local_views(ghosted_trian)) do trian + BoundaryTriangulation(trian;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end function Geometry.SkeletonTriangulation( @@ -569,32 +581,29 @@ function Geometry.SkeletonTriangulation( model = get_background_model(trian) gids = get_cell_gids(model) ghosted_trian = add_ghost_cells(trian) - trians = map(local_views(ghosted_trian),partition(gids)) do trian, gids - SkeletonTriangulation(portion,gids,trian;kwargs...) + trians = map(local_views(ghosted_trian)) do trian + SkeletonTriangulation(trian;kwargs...) end - DistributedTriangulation(trians,model) + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end -function Geometry.Triangulation( - portion,gids::AbstractLocalIndices, args...;kwargs...) +function Geometry.Triangulation(portion,gids::AbstractLocalIndices, args...;kwargs...) trian = Triangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end -function Geometry.BoundaryTriangulation( - portion,gids::AbstractLocalIndices,args...;kwargs...) +function Geometry.BoundaryTriangulation(portion,gids::AbstractLocalIndices,args...;kwargs...) trian = BoundaryTriangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end -function Geometry.SkeletonTriangulation( - portion,gids::AbstractLocalIndices,args...;kwargs...) +function Geometry.SkeletonTriangulation(portion,gids::AbstractLocalIndices,args...;kwargs...) trian = SkeletonTriangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end -function Geometry.InterfaceTriangulation( - portion,gids::AbstractLocalIndices,args...;kwargs...) +function Geometry.InterfaceTriangulation(portion,gids::AbstractLocalIndices,args...;kwargs...) trian = InterfaceTriangulation(args...;kwargs...) filter_cells_when_needed(portion,gids,trian) end @@ -605,55 +614,29 @@ function Geometry.InterfaceTriangulation(a::DistributedTriangulation,b::Distribu DistributedTriangulation(trians,a.model) end -function Geometry.Triangulation( - portion, model::DistributedDiscreteModel;kwargs...) - D = num_cell_dims(model) - Triangulation(portion,ReferenceFE{D},model;kwargs...) -end - -function Geometry.Triangulation( - ::Type{ReferenceFE{D}}, model::DistributedDiscreteModel;kwargs...) where D - Triangulation(no_ghost, ReferenceFE{D}, model; kwargs...) -end - -function filter_cells_when_needed( - portion::WithGhost, - cell_gids::AbstractLocalIndices, - trian::Triangulation) +# Filtering cells dispatch - trian +@inline function filter_cells_when_needed( + portion::Union{WithGhost,FullyAssembledRows},cell_gids,trian) + return trian end -function filter_cells_when_needed( - portion::NoGhost, - cell_gids::AbstractLocalIndices, - trian::Triangulation) - - remove_ghost_cells(trian,cell_gids) -end - -function filter_cells_when_needed( - portion::FullyAssembledRows, - cell_gids::AbstractLocalIndices, - trian::Triangulation) - - trian +@inline function filter_cells_when_needed( + portion::Union{NoGhost,SubAssembledRows},cell_gids,trian) + return remove_ghost_cells(trian,cell_gids) end -function filter_cells_when_needed( - portion::SubAssembledRows, - cell_gids::AbstractLocalIndices, - trian::Triangulation) +# Removing ghost cells - remove_ghost_cells(trian,cell_gids) +struct RemoveGhostsMetadata{A} + parents::A end -# Removing ghost cells - function remove_ghost_cells(trian::DistributedTriangulation,gids) trians = map(remove_ghost_cells,local_views(trian),partition(gids)) model = get_background_model(trian) - return DistributedTriangulation(trians,model) + metadata = RemoveGhostsMetadata(local_views(trian)) + return DistributedTriangulation(trians,model;metadata) end function remove_ghost_cells(trian::Triangulation,gids) @@ -717,7 +700,7 @@ function _find_owned_skeleton_facets(glue,gids) end tface_to_part[tface] = part end - findall(part->part==part_id(gids),tface_to_part) + return findall(isequal(part_id(gids)),tface_to_part) end # Adding ghost cells @@ -728,14 +711,19 @@ function add_ghost_cells(dtrian::DistributedTriangulation) end function add_ghost_cells(dmodel::DistributedDiscreteModel,dtrian::DistributedTriangulation) - T = eltype(local_views(dtrian)) - add_ghost_cells(T,dmodel,dtrian) + add_ghost_cells(dtrian.metadata,dmodel,dtrian) end +# We already have the parents saved up function add_ghost_cells( - ::Type{<:Triangulation}, - dmodel::DistributedDiscreteModel{Dm}, - dtrian::DistributedTriangulation{Dt} + metadata::RemoveGhostsMetadata, dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} +) where {Dm,Dt} + DistributedTriangulation(metadata.parents,dmodel) +end + +# We have to reconstruct the ghosted triangulation +function add_ghost_cells( + metadata, dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} ) where {Dm,Dt} tface_to_mface = map(local_views(dtrian)) do trian @@ -772,86 +760,6 @@ function add_ghost_cells( return DistributedTriangulation(trians,dmodel) end -function add_ghost_cells( - ::Type{<:BoundaryTriangulation}, - dmodel::DistributedDiscreteModel, - dtrian::DistributedTriangulation{Dt} -) where Dt - println("Flag 0") - bf_dtrian = add_ghost_cells(Triangulation, dmodel, dtrian) - (bf_trian === dtrian) && return dtrian # Case 1 above - - # Otherwise, we have created the underlying BodyFittedTriangulation but still have to - # create the expanded FaceToCellGlue... - - face_gids = get_face_gids(dmodel, Dt) - bgface_to_lcell = map(local_views(dtrian), partition(face_gids)) do trian, face_gids - glue = trian.glue - bgface_to_lcell = fill(Int8(1), length(face_gids)) - bgface_to_lcell[glue.face_to_bgface] .= glue.face_to_lcell - end - consistent!(PVector(bgface_to_lcell, partition(face_gids))) |> wait - - trians = map(local_views(bf_dtrian), local_views(dmodel), bgface_to_lcell) do bf_trian, model, bgface_to_lcell - topo = get_grid_topology(model) - face_grid = get_grid(bf_trian) - cell_grid = get_grid(model) - face_to_bgface = bf_trian.tface_to_mface - glue = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell) - BoundaryTriangulation(bf_trian,glue) - end - - return DistributedTriangulation(trians,dmodel) -end - -function add_ghost_cells( - ::Type{<:SkeletonTriangulation}, - dmodel::DistributedDiscreteModel, - dtrian::DistributedTriangulation{Dt} -) where Dt - bf_dtrian = add_ghost_cells(Triangulation, dmodel, dtrian) - (bf_trian === dtrian) && return dtrian # Case 1 above - - # Otherwise, we have created the underlying BodyFittedTriangulation but still have to - # create the expanded FaceToCellGlue... - - face_gids = get_face_gids(dmodel, Dt) - bgface_to_lcell_plus = map(local_views(dtrian), partition(face_gids)) do trian, face_gids - glue = trian.plus.glue - bgface_to_lcell = fill(Int8(1), length(face_gids)) - bgface_to_lcell[glue.face_to_bgface] .= glue.face_to_lcell - end - t_plus = consistent!(PVector(bgface_to_lcell, partition(face_gids))) - bgface_to_lcell_minus = map(local_views(dtrian), partition(face_gids)) do trian, face_gids - glue = trian.minus.glue - bgface_to_lcell = fill(Int8(1), length(face_gids)) - bgface_to_lcell[glue.face_to_bgface] .= glue.face_to_lcell - end - t_minus = consistent!(PVector(bgface_to_lcell, partition(face_gids))) - - wait(t_plus) - wait(t_minus) - - trians = map( - local_views(bf_dtrian), local_views(dmodel), bgface_to_lcell_plus, bgface_to_lcell_minus - ) do bf_trian, model, bgface_to_lcell_plus, bgface_to_lcell_minus - topo = get_grid_topology(model) - face_grid = get_grid(bf_trian) - cell_grid = get_grid(model) - face_to_bgface = bf_trian.tface_to_mface - - glue_plus = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell_plus) - plus = BoundaryTriangulation(bf_trian,glue_plus) - - glue_minus = FaceToCellGlue(topo,cell_grid,face_grid,face_to_bgface,bgface_to_lcell_minus) - minus = BoundaryTriangulation(bf_trian,glue_minus) - - SkeletonTriangulation(plus,minus) - end - - return DistributedTriangulation(trians,dmodel) -end - function _covers_all_faces( dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} From d04d19a3217f620f1b293c67c049c28d9435685a Mon Sep 17 00:00:00 2001 From: Jordi Manyer Date: Sat, 16 Aug 2025 12:15:14 +1000 Subject: [PATCH 4/7] Updated NEWS --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index cad1599..d9b9a1e 100644 --- a/NEWS.md +++ b/NEWS.md @@ -5,6 +5,12 @@ 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] + +### Fixed + +- Fixed issue [#177](https://github.com/gridap/GridapDistributed.jl/issues/177) and [#170](https://github.com/gridap/GridapDistributed.jl/issues/170). Since PR[#180](https://github.com/gridap/GridapDistributed.jl/pull/180). + ## [0.4.8] 2025-06-11 ### Added From 0db17a551b1aaba971a3a634d703e740405d1568 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Wed, 17 Sep 2025 12:19:58 +1000 Subject: [PATCH 5/7] added mwe --- test.jl | 37 +++++++++++++++++++++++++++++++---- test_mwe.jl | 55 +++++++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 88 insertions(+), 4 deletions(-) create mode 100644 test_mwe.jl diff --git a/test.jl b/test.jl index a971327..b896f3f 100644 --- a/test.jl +++ b/test.jl @@ -14,7 +14,22 @@ using Gridap.MultiField using GridapDistributed using PartitionedArrays -function half_empty_trian(ranks,model) +function partial_trian(ranks,model) + cell_ids = get_cell_gids(model) + trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids + cell_mask = zeros(Bool, num_cells(model)) + if rank ∈ (1,2) + cell_mask[own_to_local(ids)] .= true + else + t = own_to_local(ids) + cell_mask[t[1:floor(Int,length(t)/2)]] .= true + end + Triangulation(model,cell_mask) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + +function half_empty_trian(ranks,model) # edge case cell_ids = get_cell_gids(model) trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids cell_mask = zeros(Bool, num_cells(model)) @@ -26,15 +41,27 @@ function half_empty_trian(ranks,model) GridapDistributed.DistributedTriangulation(trians,model) end +function trian_with_empty_procs(ranks,model) + cell_ids = get_cell_gids(model) + trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids + cell_mask = zeros(Bool, num_cells(model)) + if rank ∈ (1,2) + t = own_to_local(ids) + cell_mask[t[1:floor(Int,length(t)/2)]] .= true + end + Triangulation(model,cell_mask) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + np = (2,2) ranks = with_debug() do distribute distribute(LinearIndices((prod(np),))) end -model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(5,5)) +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(10,10)) Ω = Triangulation(model) -# Γ = BoundaryTriangulation(model) -Γ = half_empty_trian(ranks,model) +Γ = BoundaryTriangulation(model) dΩ = Measure(Ω,2) dΓ = Measure(Γ,2) V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1)) @@ -63,6 +90,8 @@ J_fwd = jacobian(op,uh) @test reduce(&,map(≈,partition(J),partition(J_fwd))) # Skel +Γ = partial_trian(ranks,model) +# Γ = half_empty_trian(ranks,model) V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1),conformity=:L2) diff --git a/test_mwe.jl b/test_mwe.jl new file mode 100644 index 0000000..1825f4d --- /dev/null +++ b/test_mwe.jl @@ -0,0 +1,55 @@ +using Pkg; Pkg.activate("./GridapDistributed"); + +using Test + +using Gridap +using Gridap.Algebra +using Gridap.Arrays +using Gridap.Geometry +using Gridap.CellData +using Gridap.ReferenceFEs +using Gridap.FESpaces +using Gridap.MultiField + +using GridapDistributed +using PartitionedArrays + +function half_empty_trian(ranks,model) # edge case + cell_ids = get_cell_gids(model) + trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids + cell_mask = zeros(Bool, num_cells(model)) + if rank ∈ (1,2) + cell_mask[own_to_local(ids)] .= true + end + Triangulation(model,cell_mask) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + +np = (2,2) +ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) +end + +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(5,5)) +Γ = half_empty_trian(ranks,model) +V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) +V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) +V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1),conformity=:L2) +X = MultiFieldFESpace([V1,V2,V3]) +uh = zero(X); +Λ = SkeletonTriangulation(model) +dΛ = Measure(Λ,2) + +# Skel jac +f2(xh,yh) = ∫(mean(xh[1])⋅mean(yh[1])+mean(xh[2])⋅mean(yh[2])+mean(xh[1])⋅mean(xh[2])⋅mean(yh[2])+mean(xh[1])*mean(xh[3])*mean(yh[3]))dΛ +dv = get_fe_basis(X); +j = jacobian(uh->f2(uh,dv),uh); +J = assemble_matrix(j,X,X) + +f2_jac(xh,dxh,yh) = ∫(mean(dxh[1])⋅mean(yh[1])+mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])⋅mean(xh[2])⋅mean(yh[2]) + + mean(xh[1])⋅mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])*mean(xh[3])*mean(yh[3])+mean(xh[1])*mean(dxh[3])*mean(yh[3]))dΛ +op = FEOperator(f2,f2_jac,X,X) +J_fwd = jacobian(op,uh) + +@test reduce(&,map(≈,partition(J),partition(J_fwd))) \ No newline at end of file From 8c72f0afcf48476213b3456fbf26e05f574f4e1a Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 19 Sep 2025 19:56:49 +1000 Subject: [PATCH 6/7] testing skel --- src/Geometry.jl | 169 ++++++++++++++++++++++++++++++++---------- test.jl | 122 ------------------------------ test/AutodiffTests.jl | 108 ++++++++++++++++++++++++++- test_mwe.jl | 4 +- 4 files changed, 239 insertions(+), 164 deletions(-) delete mode 100644 test.jl diff --git a/src/Geometry.jl b/src/Geometry.jl index a96dabb..80c58f5 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -38,14 +38,27 @@ Geometry.num_point_dims(::Type{<:DistributedGrid{Dc,Dp}}) where {Dc,Dp} = Dp # This object cannot implement the GridTopology interface in a strict sense """ """ -struct DistributedGridTopology{Dc,Dp,A} <: GridapType +struct DistributedGridTopology{Dc,Dp,A,B} <: GridapType topos::A - function DistributedGridTopology(topos::AbstractArray{<:GridTopology{Dc,Dp}}) where {Dc,Dp} + face_gids::B + function DistributedGridTopology( + topos::AbstractArray{<:GridTopology{Dc,Dp}}, + face_gids::AbstractArray{<:PRange} + ) where {Dc,Dp} A = typeof(topos) - new{Dc,Dp,A}(topos) + B = typeof(face_gids) + new{Dc,Dp,A,B}(topos,face_gids) end end +function DistributedGridTopology( + topos::AbstractArray{<:GridTopology{Dc,Dp}}, cell_gids::PRange +) where {Dc,Dp} + face_gids = Vector{PRange}(undef,Dc+1) + face_gids[Dc+1] = cell_gids + return DistributedGridTopology(topos,face_gids) +end + local_views(a::DistributedGridTopology) = a.topos function Geometry.OrientationStyle( @@ -63,6 +76,39 @@ Geometry.num_cell_dims(::Type{<:DistributedGridTopology{Dc,Dp}}) where {Dc,Dp} = Geometry.num_point_dims(::DistributedGridTopology{Dc,Dp}) where {Dc,Dp} = Dp Geometry.num_point_dims(::Type{<:DistributedGridTopology{Dc,Dp}}) where {Dc,Dp} = Dp +function get_cell_gids(topo::DistributedGridTopology{Dc}) where Dc + topo.face_gids[Dc+1] +end + +function get_face_gids(topo::DistributedGridTopology,dim::Integer) + _setup_face_gids!(topo,dim) + return topo.face_gids[dim+1] +end + +function _setup_face_gids!(topo::DistributedGridTopology{Dc},dim) where {Dc} + Gridap.Helpers.@check 0 <= dim <= Dc + if !isassigned(topo.face_gids,dim+1) + cell_gids = topo.face_gids[Dc+1] + nlfaces = map(local_views(topo)) do topo + num_faces(topo,dim) + end + cell_lfaces = map(local_views(topo)) do topo + get_faces(topo, Dc, dim) + end + topo.face_gids[dim+1] = generate_gids(cell_gids,cell_lfaces,nlfaces) + end +end + +function Geometry.get_isboundary_face(topo::DistributedGridTopology, d::Integer) + face_gids = get_face_gids(topo, d) + is_local_boundary = map(local_views(topo)) do topo + get_isboundary_face(topo,d) + end + t = assemble!(&,PVector(is_local_boundary, partition(face_gids))) + is_global_boundary = partition(fetch(consistent!(fetch(t)))) + return is_global_boundary +end + """ """ struct DistributedFaceLabeling{A<:AbstractArray{<:FaceLabeling}} @@ -72,12 +118,40 @@ end local_views(a::DistributedFaceLabeling) = a.labels function Geometry.add_tag_from_tags!(labels::DistributedFaceLabeling, name, tags) - map(labels.labels) do labels + map(local_views(labels)) do labels add_tag_from_tags!(labels, name, tags) end end -# Dsitributed Discrete models +function Geometry.get_face_mask(labels::DistributedFaceLabeling, tags, d::Integer) + map(local_views(labels)) do labels + get_face_mask(labels, tags, d) + end +end + +function Geometry.FaceLabeling(topo::DistributedGridTopology) + D = num_cell_dims(topo) + labels = map(local_views(topo)) do topo + d_to_ndfaces = [ num_faces(topo,d) for d in 0:D ] + labels = FaceLabeling(d_to_ndfaces) + for d in 0:D + get_face_entity(labels,d) .= 1 # Interior as default + end + add_tag!(labels,"interior",[1]) + add_tag!(labels,"boundary",[2]) + return labels + end + for d in 0:D-1 + dface_to_is_boundary = get_isboundary_face(topo,d) # Global boundary + map(labels,dface_to_is_boundary) do labels, dface_to_is_boundary + dface_to_entity = get_face_entity(labels,d) + dface_to_entity .+= dface_to_is_boundary + end + end + return labels +end + +# Distributed Discrete models # We do not inherit from DiscreteModel on purpose. # This object cannot implement the DiscreteModel interface in a strict sense @@ -131,7 +205,7 @@ function Geometry.get_grid(model::DistributedDiscreteModel) end function Geometry.get_grid_topology(model::DistributedDiscreteModel) - DistributedGridTopology(map(get_grid_topology,local_views(model))) + DistributedGridTopology(map(get_grid_topology,local_views(model)),model.face_gids) end function Geometry.get_face_labeling(model::DistributedDiscreteModel) @@ -320,7 +394,7 @@ function _cartesian_model_with_periodic_bcs(pdesc::DistributedCartesianDescripto remove_boundary = map((p,n)->((p && (n!=1)) ? true : false),desc.isperiodic,parts) CartesianDiscreteModel(_desc,cmin,cmax,remove_boundary) end - + return models, global_partition end @@ -465,7 +539,7 @@ struct DistributedTriangulation{Dc,Dp,A,B,C} <: GridapType function DistributedTriangulation( trians::AbstractArray{<:Triangulation{Dc,Dp}}, model::DistributedDiscreteModel; - metadata = nothing + metadata = nothing ) where {Dc,Dp} A = typeof(trians) B = typeof(model) @@ -510,55 +584,74 @@ function Geometry.Triangulation(model::DistributedDiscreteModel;kwargs...) Triangulation(no_ghost,ReferenceFE{D},model;kwargs...) end -function Geometry.Triangulation(::Type{ReferenceFE{D}}, model::DistributedDiscreteModel;kwargs...) where D +function Geometry.Triangulation(::Type{ReferenceFE{D}},model::DistributedDiscreteModel;kwargs...) where D Triangulation(no_ghost, ReferenceFE{D}, model; kwargs...) end -function Geometry.BoundaryTriangulation(model::DistributedDiscreteModel;kwargs...) - BoundaryTriangulation(no_ghost,model;kwargs...) +function Geometry.Triangulation(portion, model::DistributedDiscreteModel;kwargs...) + D = num_cell_dims(model) + Triangulation(portion,ReferenceFE{D},model;kwargs...) end -function Geometry.BoundaryTriangulation(trian::DistributedTriangulation;kwargs...) - BoundaryTriangulation(no_ghost,trian;kwargs...) +function Geometry.Triangulation( + portion,::Type{ReferenceFE{D}},model::DistributedDiscreteModel;kwargs...) where D + gids = get_face_gids(model,D) + trians = map(local_views(model)) do model + Triangulation(ReferenceFE{D},model;kwargs...) + end + parent = DistributedTriangulation(trians,model) + return filter_cells_when_needed(portion,gids,parent) end -function Geometry.SkeletonTriangulation(model::DistributedDiscreteModel;kwargs...) - SkeletonTriangulation(no_ghost,model;kwargs...) +function Geometry.BoundaryTriangulation(model::DistributedDiscreteModel,args...;kwargs...) + BoundaryTriangulation(no_ghost,model,args...;kwargs...) end -function Geometry.SkeletonTriangulation(trian::DistributedTriangulation;kwargs...) - SkeletonTriangulation(no_ghost,trian;kwargs...) +function Geometry.BoundaryTriangulation(trian::DistributedTriangulation;kwargs...) + BoundaryTriangulation(no_ghost,trian;kwargs...) end -function Geometry.Triangulation(portion, model::DistributedDiscreteModel;kwargs...) - D = num_cell_dims(model) - Triangulation(portion,ReferenceFE{D},model;kwargs...) +function Geometry.BoundaryTriangulation( + portion,model::DistributedDiscreteModel;kwargs...) + labels = get_face_labeling(model) + Geometry.BoundaryTriangulation(portion,model,labels;kwargs...) end -function Geometry.Triangulation( - portion,::Type{ReferenceFE{Dt}},model::DistributedDiscreteModel{Dm};kwargs...) where {Dt,Dm} - # Generate global ordering for the faces of dimension Dt (if needed) - gids = get_face_gids(model,Dt) - trians = map(local_views(model)) do model - Triangulation(ReferenceFE{Dt},model;kwargs...) +function Geometry.BoundaryTriangulation( + portion,model::DistributedDiscreteModel,labels::DistributedFaceLabeling;tags=nothing) + Dc = num_cell_dims(model) + if isnothing(tags) + topo = get_grid_topology(model) + face_to_mask = get_isboundary_face(topo,Dc-1) # This is globally consistent + else + face_to_mask = get_face_mask(labels,tags,Dc-1) end - parent = DistributedTriangulation(trians,model) - return filter_cells_when_needed(portion,gids,parent) + Geometry.BoundaryTriangulation(portion,model,face_to_mask) end function Geometry.BoundaryTriangulation( - portion,model::DistributedDiscreteModel{Dc};kwargs...) where Dc - gids = get_face_gids(model,Dc) - trians = map(local_views(model)) do model - BoundaryTriangulation(model;kwargs...) + portion,model::DistributedDiscreteModel,face_to_mask::AbstractArray) + Dc = num_cell_dims(model) + gids = get_face_gids(model,Dc) + trians = map(local_views(model),face_to_mask) do model, face_to_mask + BoundaryTriangulation(model,face_to_mask) end parent = DistributedTriangulation(trians,model) return filter_cells_when_needed(portion,gids,parent) end +function Geometry.SkeletonTriangulation(model::DistributedDiscreteModel;kwargs...) + SkeletonTriangulation(no_ghost,model;kwargs...) +end + +function Geometry.SkeletonTriangulation(trian::DistributedTriangulation;kwargs...) + SkeletonTriangulation(no_ghost,trian;kwargs...) +end + function Geometry.SkeletonTriangulation( - portion,model::DistributedDiscreteModel{Dc};kwargs...) where Dc - gids = get_face_gids(model,Dc) + portion,model::DistributedDiscreteModel;kwargs...) + Dc = num_cell_dims(model) + gids = get_face_gids(model,Dc) trians = map(local_views(model)) do model SkeletonTriangulation(model;kwargs...) end @@ -566,7 +659,7 @@ function Geometry.SkeletonTriangulation( return filter_cells_when_needed(portion,gids,parent) end -# NOTE: The following constructors require adding back the ghost cells: +# NOTE: The following constructors require adding back the ghost cells: # Potentially, the input `trian` has had some/all of its ghost cells removed. If we do not # add them back, some skeleton facets might look like boundary facets to the local constructors... function Geometry.BoundaryTriangulation( @@ -621,7 +714,7 @@ function Geometry.InterfaceTriangulation(a::DistributedTriangulation,b::Distribu DistributedTriangulation(trians,a.model) end -# Filtering cells dispatch +# Filtering cells @inline function filter_cells_when_needed( portion::Union{WithGhost,FullyAssembledRows},cell_gids,trian) @@ -721,7 +814,7 @@ function add_ghost_cells(dmodel::DistributedDiscreteModel,dtrian::DistributedTri add_ghost_cells(dtrian.metadata,dmodel,dtrian) end -# We already have the parents saved up +# We already have the parents saved up function add_ghost_cells( metadata::RemoveGhostsMetadata, dmodel::DistributedDiscreteModel{Dm}, dtrian::DistributedTriangulation{Dt} ) where {Dm,Dt} @@ -844,4 +937,4 @@ function generate_cell_gids(dmodel::DistributedDiscreteModel{Dm}, tgids = PRange(indices) end return tgids -end +end \ No newline at end of file diff --git a/test.jl b/test.jl deleted file mode 100644 index b896f3f..0000000 --- a/test.jl +++ /dev/null @@ -1,122 +0,0 @@ -using Pkg; Pkg.activate("./GridapDistributed"); - -using Test - -using Gridap -using Gridap.Algebra -using Gridap.Arrays -using Gridap.Geometry -using Gridap.CellData -using Gridap.ReferenceFEs -using Gridap.FESpaces -using Gridap.MultiField - -using GridapDistributed -using PartitionedArrays - -function partial_trian(ranks,model) - cell_ids = get_cell_gids(model) - trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids - cell_mask = zeros(Bool, num_cells(model)) - if rank ∈ (1,2) - cell_mask[own_to_local(ids)] .= true - else - t = own_to_local(ids) - cell_mask[t[1:floor(Int,length(t)/2)]] .= true - end - Triangulation(model,cell_mask) - end - GridapDistributed.DistributedTriangulation(trians,model) -end - -function half_empty_trian(ranks,model) # edge case - cell_ids = get_cell_gids(model) - trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids - cell_mask = zeros(Bool, num_cells(model)) - if rank ∈ (3,4) - cell_mask[own_to_local(ids)] .= true - end - Triangulation(model,cell_mask) - end - GridapDistributed.DistributedTriangulation(trians,model) -end - -function trian_with_empty_procs(ranks,model) - cell_ids = get_cell_gids(model) - trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids - cell_mask = zeros(Bool, num_cells(model)) - if rank ∈ (1,2) - t = own_to_local(ids) - cell_mask[t[1:floor(Int,length(t)/2)]] .= true - end - Triangulation(model,cell_mask) - end - GridapDistributed.DistributedTriangulation(trians,model) -end - -np = (2,2) -ranks = with_debug() do distribute - distribute(LinearIndices((prod(np),))) -end - -model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(10,10)) -Ω = Triangulation(model) -Γ = BoundaryTriangulation(model) -dΩ = Measure(Ω,2) -dΓ = Measure(Γ,2) -V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1)) -V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1)) -V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1)) -X = MultiFieldFESpace([V1,V2,V3]) -uh = zero(X); - -f(xh) = ∫(xh[1]+xh[2]⋅xh[2]+xh[1]*xh[3])dΓ -df(v,xh) = ∫(v[1]+2*v[2]⋅xh[2]+v[1]*xh[3]+xh[1]*v[3])dΓ -du = gradient(f,uh) -du_vec = assemble_vector(du,X) -df_vec = assemble_vector(v->df(v,uh),X) - -@test df_vec ≈ du_vec - -f2(xh,yh) = ∫(xh[1]⋅yh[1]+xh[2]⋅yh[2]+xh[1]⋅xh[2]⋅yh[2]+xh[1]*xh[3]*yh[3])dΓ -dv = get_fe_basis(X); -j = jacobian(uh->f2(uh,dv),uh) -J = assemble_matrix(j,X,X) - -f2_jac(xh,dxh,yh) = ∫(dxh[1]⋅yh[1]+dxh[2]⋅yh[2]+dxh[1]⋅xh[2]⋅yh[2]+xh[1]⋅dxh[2]⋅yh[2]+dxh[1]*xh[3]*yh[3]+xh[1]*dxh[3]*yh[3])dΓ -op = FEOperator(f2,f2_jac,X,X) -J_fwd = jacobian(op,uh) - -@test reduce(&,map(≈,partition(J),partition(J_fwd))) - -# Skel -Γ = partial_trian(ranks,model) -# Γ = half_empty_trian(ranks,model) -V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) -V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) -V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1),conformity=:L2) -X = MultiFieldFESpace([V1,V2,V3]) -uh = zero(X); -Λ = SkeletonTriangulation(model) -dΛ = Measure(Λ,2) - -f(xh) = ∫(mean(xh[1])+mean(xh[2])⋅mean(xh[2])+mean(xh[1])*mean(xh[3]))dΛ -df(v,xh) = ∫(mean(v[1])+2*mean(v[2])⋅mean(xh[2])+mean(v[1])*mean(xh[3])+mean(xh[1])*mean(v[3]))dΛ -du = gradient(f,uh) -du_vec = assemble_vector(du,X) -df_vec = assemble_vector(v->df(v,uh),X) - -@test df_vec ≈ du_vec - -# Skel jac -f2(xh,yh) = ∫(mean(xh[1])⋅mean(yh[1])+mean(xh[2])⋅mean(yh[2])+mean(xh[1])⋅mean(xh[2])⋅mean(yh[2])+mean(xh[1])*mean(xh[3])*mean(yh[3]))dΛ -dv = get_fe_basis(X); -j = jacobian(uh->f2(uh,dv),uh); -J = assemble_matrix(j,X,X) - -f2_jac(xh,dxh,yh) = ∫(mean(dxh[1])⋅mean(yh[1])+mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])⋅mean(xh[2])⋅mean(yh[2]) + - mean(xh[1])⋅mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])*mean(xh[3])*mean(yh[3])+mean(xh[1])*mean(dxh[3])*mean(yh[3]))dΛ -op = FEOperator(f2,f2_jac,X,X) -J_fwd = jacobian(op,uh) - -@test reduce(&,map(≈,partition(J),partition(J_fwd))) \ No newline at end of file diff --git a/test/AutodiffTests.jl b/test/AutodiffTests.jl index 20e6795..f86a79d 100644 --- a/test/AutodiffTests.jl +++ b/test/AutodiffTests.jl @@ -44,7 +44,7 @@ function main_sf(distribute,parts) @test b ≈ b_AD # Skeleton AD - # I would like to compare the results, but we cannot be using FD in parallel... + # I would like to compare the results, but we cannot be using FD in parallel... Λ = SkeletonTriangulation(model) dΛ = Measure(Λ,2*k) g_Λ(v) = ∫(mean(v))*dΛ @@ -73,7 +73,7 @@ function main_mf(distribute,parts) Ω = Triangulation(model) dΩ = Measure(Ω,2*(k+1)) - + ν = 1.0 f = VectorValue(0.0,0.0) @@ -104,9 +104,113 @@ function main_mf(distribute,parts) @test b ≈ b_AD end +## MultiField AD with different triangulations for each field +function generate_trian(ranks,model,case) + cell_ids = get_cell_gids(model) + trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids + cell_mask = zeros(Bool, num_cells(model)) + if case == :partial_trian + if rank ∈ (1,2) + cell_mask[own_to_local(ids)] .= true + else + t = own_to_local(ids) + cell_mask[t[1:floor(Int,length(t)/2)]] .= true + end + elseif case == :half_empty_trian + if rank ∈ (3,4) + cell_mask[own_to_local(ids)] .= true + end + elseif case == :trian_with_empty_procs + if rank ∈ (1,2) + t = own_to_local(ids) + cell_mask[t[1:floor(Int,length(t)/2)]] .= true + end + else + error("Unknown case") + end + Triangulation(model,cell_mask) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + +function mf_different_fespace_trians(distribute,parts) + ranks = distribute(LinearIndices((prod(parts),))) + model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(10,10)) + V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1)) + V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1)) + for case in (:boundary,:partial_trian, :half_empty_trian, :trian_with_empty_procs) + if case == :boundary + Γ = BoundaryTriangulation(model) + else + Γ = generate_trian(ranks,model,case) + end + dΓ = Measure(Γ,2) + V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1)) + X = MultiFieldFESpace([V1,V2,V3]) + uh = zero(X); + + f(xh) = ∫(xh[1]+xh[2]⋅xh[2]+xh[1]*xh[3])dΓ + df(v,xh) = ∫(v[1]+2*v[2]⋅xh[2]+v[1]*xh[3]+xh[1]*v[3])dΓ + du = gradient(f,uh) + du_vec = assemble_vector(du,X) + df_vec = assemble_vector(v->df(v,uh),X) + + @test df_vec ≈ du_vec + + f2(xh,yh) = ∫(xh[1]⋅yh[1]+xh[2]⋅yh[2]+xh[1]⋅xh[2]⋅yh[2]+xh[1]*xh[3]*yh[3])dΓ + dv = get_fe_basis(X); + j = jacobian(uh->f2(uh,dv),uh) + J = assemble_matrix(j,X,X) + + f2_jac(xh,dxh,yh) = ∫(dxh[1]⋅yh[1]+dxh[2]⋅yh[2]+dxh[1]⋅xh[2]⋅yh[2]+xh[1]⋅dxh[2]⋅yh[2]+dxh[1]*xh[3]*yh[3]+xh[1]*dxh[3]*yh[3])dΓ + op = FEOperator(f2,f2_jac,X,X) + J_fwd = jacobian(op,uh) + + @test reduce(&,map(≈,partition(J),partition(J_fwd))) + end +end + +function skeleton_mf_different_fespace_trians(distribute,parts) + ranks = distribute(LinearIndices((prod(parts),))) + model = CartesianDiscreteModel(ranks,parts,(0,1,0,1),(10,10)) + for case in (:partial_trian, :half_empty_trian, :trian_with_empty_procs) + Γ = generate_trian(ranks,model,case) + V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) + V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) + V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1),conformity=:L2) + X = MultiFieldFESpace([V1,V2,V3]) + uh = zero(X); + Λ = SkeletonTriangulation(model) + dΛ = Measure(Λ,2) + + f(xh) = ∫(mean(xh[1])+mean(xh[2])⋅mean(xh[2])+mean(xh[1])*mean(xh[3]))dΛ + df(v,xh) = ∫(mean(v[1])+2*mean(v[2])⋅mean(xh[2])+mean(v[1])*mean(xh[3])+mean(xh[1])*mean(v[3]))dΛ + du = gradient(f,uh) + du_vec = assemble_vector(du,X) + df_vec = assemble_vector(v->df(v,uh),X) + + @test df_vec ≈ du_vec + + # Skel jac + f2(xh,yh) = ∫(mean(xh[1])⋅mean(yh[1])+mean(xh[2])⋅mean(yh[2])+mean(xh[1])⋅mean(xh[2])⋅mean(yh[2])+mean(xh[1])*mean(xh[3])*mean(yh[3]))dΛ + dv = get_fe_basis(X); + j = jacobian(uh->f2(uh,dv),uh); + J = assemble_matrix(j,X,X) + + f2_jac(xh,dxh,yh) = ∫(mean(dxh[1])⋅mean(yh[1])+mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])⋅mean(xh[2])⋅mean(yh[2]) + + mean(xh[1])⋅mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])*mean(xh[3])*mean(yh[3])+mean(xh[1])*mean(dxh[3])*mean(yh[3]))dΛ + op = FEOperator(f2,f2_jac,X,X) + J_fwd = jacobian(op,uh) + + @test reduce(&,map(≈,partition(J),partition(J_fwd))) + end +end + function main(distribute,parts) main_sf(distribute,parts) main_mf(distribute,parts) + mf_different_fespace_trians(distribute,parts) + skeleton_mf_different_fespace_trians(distribute,parts) end end \ No newline at end of file diff --git a/test_mwe.jl b/test_mwe.jl index 1825f4d..3a89942 100644 --- a/test_mwe.jl +++ b/test_mwe.jl @@ -18,7 +18,7 @@ function half_empty_trian(ranks,model) # edge case cell_ids = get_cell_gids(model) trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids cell_mask = zeros(Bool, num_cells(model)) - if rank ∈ (1,2) + if rank ∈ (1,2) # (1,2) cell_mask[own_to_local(ids)] .= true end Triangulation(model,cell_mask) @@ -31,7 +31,7 @@ ranks = with_debug() do distribute distribute(LinearIndices((prod(np),))) end -model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(5,5)) +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(10,10)) Γ = half_empty_trian(ranks,model) V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) From 8a55f9ebbbc05193115e6558e334cbafdf5e6021 Mon Sep 17 00:00:00 2001 From: zjwegert <60646897+zjwegert@users.noreply.github.com> Date: Fri, 19 Sep 2025 22:48:19 +1000 Subject: [PATCH 7/7] cleanup --- test_mwe.jl | 55 ----------------------------------------------------- 1 file changed, 55 deletions(-) delete mode 100644 test_mwe.jl diff --git a/test_mwe.jl b/test_mwe.jl deleted file mode 100644 index 3a89942..0000000 --- a/test_mwe.jl +++ /dev/null @@ -1,55 +0,0 @@ -using Pkg; Pkg.activate("./GridapDistributed"); - -using Test - -using Gridap -using Gridap.Algebra -using Gridap.Arrays -using Gridap.Geometry -using Gridap.CellData -using Gridap.ReferenceFEs -using Gridap.FESpaces -using Gridap.MultiField - -using GridapDistributed -using PartitionedArrays - -function half_empty_trian(ranks,model) # edge case - cell_ids = get_cell_gids(model) - trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids - cell_mask = zeros(Bool, num_cells(model)) - if rank ∈ (1,2) # (1,2) - cell_mask[own_to_local(ids)] .= true - end - Triangulation(model,cell_mask) - end - GridapDistributed.DistributedTriangulation(trians,model) -end - -np = (2,2) -ranks = with_debug() do distribute - distribute(LinearIndices((prod(np),))) -end - -model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(10,10)) -Γ = half_empty_trian(ranks,model) -V1 = FESpace(Γ,ReferenceFE(lagrangian,Float64,1),conformity=:L2) -V2 = FESpace(model,ReferenceFE(lagrangian,VectorValue{2,Float64},1),conformity=:L2) -V3 = FESpace(model,ReferenceFE(lagrangian,Float64,1),conformity=:L2) -X = MultiFieldFESpace([V1,V2,V3]) -uh = zero(X); -Λ = SkeletonTriangulation(model) -dΛ = Measure(Λ,2) - -# Skel jac -f2(xh,yh) = ∫(mean(xh[1])⋅mean(yh[1])+mean(xh[2])⋅mean(yh[2])+mean(xh[1])⋅mean(xh[2])⋅mean(yh[2])+mean(xh[1])*mean(xh[3])*mean(yh[3]))dΛ -dv = get_fe_basis(X); -j = jacobian(uh->f2(uh,dv),uh); -J = assemble_matrix(j,X,X) - -f2_jac(xh,dxh,yh) = ∫(mean(dxh[1])⋅mean(yh[1])+mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])⋅mean(xh[2])⋅mean(yh[2]) + - mean(xh[1])⋅mean(dxh[2])⋅mean(yh[2])+mean(dxh[1])*mean(xh[3])*mean(yh[3])+mean(xh[1])*mean(dxh[3])*mean(yh[3]))dΛ -op = FEOperator(f2,f2_jac,X,X) -J_fwd = jacobian(op,uh) - -@test reduce(&,map(≈,partition(J),partition(J_fwd))) \ No newline at end of file