From 09db08df1b0694d31994c82a3051300655a24023 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 14 Feb 2023 11:27:07 +1100 Subject: [PATCH 01/14] Minor changes --- src/PatchBasedSmoothers/mpi/PatchDecompositions.jl | 14 +++++++------- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 2 +- src/PatchBasedSmoothers/seq/PatchFESpaces.jl | 14 +++++++------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl index 609b0028..c208ec4e 100644 --- a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl @@ -7,27 +7,27 @@ end function PatchDecomposition(model::GridapDistributed.DistributedDiscreteModel{Dc,Dp}; Dr=0, patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude()) where {Dc,Dp} - patch_decompositions=map_parts(model.models) do lmodel + patch_decompositions = map_parts(model.models) do lmodel PatchDecomposition(lmodel; Dr=Dr, patch_boundary_style=patch_boundary_style) end - A=typeof(patch_decompositions) - B=typeof(model) + A = typeof(patch_decompositions) + B = typeof(model) DistributedPatchDecomposition{Dc,Dp,A,B}(patch_decompositions,model) end function Gridap.Geometry.Triangulation(a::DistributedPatchDecomposition) - trians=map_parts(a.patch_decompositions) do a + trians = map_parts(a.patch_decompositions) do a Triangulation(a) end GridapDistributed.DistributedTriangulation(trians,a.model) end function get_patch_root_dim(a::DistributedPatchDecomposition) - patch_root_dim=0 + patch_root_dim = 0 map_parts(a.patch_decompositions) do patch_decomposition - patch_root_dim=patch_decomposition.Dr + patch_root_dim = patch_decomposition.Dr end - patch_root_dim + return patch_root_dim end diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index 76924b8b..4778e82c 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -12,7 +12,7 @@ function PatchFESpace(model::GridapDistributed.DistributedDiscreteModel, conformity::Gridap.FESpaces.Conformity, patch_decomposition::DistributedPatchDecomposition, Vh::GridapDistributed.DistributedSingleFieldFESpace) - root_gids=get_face_gids(model,get_patch_root_dim(patch_decomposition)) + root_gids = get_face_gids(model,get_patch_root_dim(patch_decomposition)) function f(model,patch_decomposition,Vh,partition) patches_mask = fill(false,length(partition.lid_to_gid)) diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl index efaa1802..22235ee5 100644 --- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl @@ -71,13 +71,13 @@ function PatchFESpace(model::DiscreteModel, return PatchFESpace(num_dofs,patch_cell_dofs_ids,Vh,patch_decomposition) end -Gridap.FESpaces.get_dof_value_type(a::PatchFESpace)=Gridap.FESpaces.get_dof_value_type(a.Vh) -Gridap.FESpaces.get_free_dof_ids(a::PatchFESpace)=Base.OneTo(a.num_dofs) -Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace)=a.patch_cell_dofs_ids -Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace,::Triangulation)=a.patch_cell_dofs_ids -Gridap.FESpaces.get_fe_basis(a::PatchFESpace)=get_fe_basis(a.Vh) -Gridap.FESpaces.ConstraintStyle(a::PatchFESpace)=Gridap.FESpaces.UnConstrained() -Gridap.FESpaces.get_vector_type(a::PatchFESpace)=get_vector_type(a.Vh) +Gridap.FESpaces.get_dof_value_type(a::PatchFESpace) = Gridap.FESpaces.get_dof_value_type(a.Vh) +Gridap.FESpaces.get_free_dof_ids(a::PatchFESpace) = Base.OneTo(a.num_dofs) +Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace) = a.patch_cell_dofs_ids +Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace,::Triangulation) = a.patch_cell_dofs_ids +Gridap.FESpaces.get_fe_basis(a::PatchFESpace) = get_fe_basis(a.Vh) +Gridap.FESpaces.ConstraintStyle(::PatchFESpace) = Gridap.FESpaces.UnConstrained() +Gridap.FESpaces.get_vector_type(a::PatchFESpace) = get_vector_type(a.Vh) function Gridap.FESpaces.scatter_free_and_dirichlet_values(f::PatchFESpace,free_values,dirichlet_values) cell_vals = Gridap.Fields.PosNegReindex(free_values,dirichlet_values) From 90e4db8eca85dc0f77338e93eb15008c873ac3df Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 15 Feb 2023 14:58:36 +1100 Subject: [PATCH 02/14] Weights and contributions are now redistributed --- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 43 ++++++++++-------- src/PatchBasedSmoothers/seq/PatchFESpaces.jl | 41 ++++++++++------- test/seq/DistributedPatchFESpacesTests.jl | 48 ++++++++++++++++++++ 3 files changed, 98 insertions(+), 34 deletions(-) create mode 100644 test/seq/DistributedPatchFESpacesTests.jl diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index 4778e82c..be1f53ab 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -1,6 +1,6 @@ # Rationale behind distributed PatchFESpace: # 1. Patches have an owner. Only owners compute subspace correction. -# If am not owner of a patch, all dofs in my patch become -1. +# If am not owner of a patch, all dofs in my patch become -1. [DONE] # 2. Subspace correction on an owned patch may affect DoFs which # are non-owned. These corrections should be sent to the owner # process. I.e., NO -> O (reversed) communication. [PENDING] @@ -48,10 +48,7 @@ end function prolongate!(x::PVector, Ph::GridapDistributed.DistributedSingleFieldFESpace, y::PVector) - parts=get_part_ids(x.owned_values) - Gridap.Helpers.@notimplementedif num_parts(parts)!=1 - - map_parts(x.owned_values,Ph.spaces,y.owned_values) do x,Ph,y + map_parts(x.values,Ph.spaces,y.values) do x,Ph,y prolongate!(x,Ph,y) end end @@ -59,22 +56,32 @@ end function inject!(x::PVector, Ph::GridapDistributed.DistributedSingleFieldFESpace, y::PVector, - w::PVector) - parts = get_part_ids(x.owned_values) - Gridap.Helpers.@notimplementedif num_parts(parts)!=1 - - map_parts(x.owned_values,Ph.spaces,y.owned_values,w.owned_values) do x,Ph,y,w - inject!(x,Ph,y,w) + w::PVector, + w_sums::PVector) + + map_parts(x.values,Ph.spaces,y.values,w.values,w_sums.values) do x,Ph,y,w,w_sums + inject!(x,Ph,y,w,w_sums) end + + # Exchange local contributions + assemble!(x) + exchange!(x) # TO CONSIDER: Is this necessary? Do we need ghosts for later? + return x end -function compute_weight_operators(Ph::GridapDistributed.DistributedSingleFieldFESpace) - parts = get_part_ids(Ph.spaces) - Gridap.Helpers.@notimplementedif num_parts(parts) != 1 - +function compute_weight_operators(Ph::GridapDistributed.DistributedSingleFieldFESpace,Vh) + # Local weights and partial sums w = PVector(0.0,Ph.gids) - map_parts(w.owned_values,Ph.spaces) do w, Ph - w .= compute_weight_operators(Ph) + w_sums = PVector(0.0,Vh.gids) + map_parts(w.values,w_sums.values,Ph.spaces) do w, w_sums, Ph + _w, _w_sums = compute_weight_operators(Ph) + w .= _w + w_sums .= _w_sums end - return w + + # partial sums -> global sums + assemble!(w_sums) # ghost -> owners + exchange!(w_sums) # repopulate ghosts with owner info + + return w, w_sums end diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl index 22235ee5..852e4c22 100644 --- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl @@ -78,6 +78,7 @@ Gridap.FESpaces.get_cell_dof_ids(a::PatchFESpace,::Triangulation) = a.patch_cell Gridap.FESpaces.get_fe_basis(a::PatchFESpace) = get_fe_basis(a.Vh) Gridap.FESpaces.ConstraintStyle(::PatchFESpace) = Gridap.FESpaces.UnConstrained() Gridap.FESpaces.get_vector_type(a::PatchFESpace) = get_vector_type(a.Vh) +Gridap.FESpaces.get_fe_dof_basis(a::PatchFESpace) = get_fe_dof_basis(a.Vh) function Gridap.FESpaces.scatter_free_and_dirichlet_values(f::PatchFESpace,free_values,dirichlet_values) cell_vals = Gridap.Fields.PosNegReindex(free_values,dirichlet_values) @@ -140,7 +141,7 @@ end # TO-THINK/STRESS: # 1. MultiFieldFESpace case? -# 2. FESpaces which are directly defined on physical space? We think this cased is covered by +# 2. FESpaces which are directly defined on physical space? We think this case is covered by # the fact that we are using a CellConformity instance to rely on ownership info. # free_dofs_offset : the ID from which we start to assign free DoF IDs upwards # Note: we do not actually need to generate a global numbering for Dirichlet DoFs. We can @@ -255,8 +256,8 @@ end # x \in SingleFESpace # y \in PatchFESpace function inject!(x,Ph::PatchFESpace,y) - w = compute_weight_operators(Ph) - inject!(x,Ph::PatchFESpace,y,w) + w, w_sums = compute_weight_operators(Ph) + inject!(x,Ph::PatchFESpace,y,w,w_sums) end function inject!(x,Ph::PatchFESpace,y,w) @@ -288,14 +289,14 @@ function inject!(x,Ph::PatchFESpace,y,w) end end -function compute_weight_operators(Ph::PatchFESpace) - cell_dof_ids = get_cell_dof_ids(Ph.Vh) - cache_cell_dof_ids = array_cache(cell_dof_ids) - cache_patch_cells = array_cache(Ph.patch_decomposition.patch_cells) - - w = zeros(num_free_dofs(Ph.Vh)) +function inject!(x,Ph::PatchFESpace,y,w,w_sums) touched = Dict{Int,Bool}() cell_mesh_overlapped = 1 + cache_patch_cells = array_cache(Ph.patch_decomposition.patch_cells) + cell_dof_ids = get_cell_dof_ids(Ph.Vh) + cache_cell_dof_ids = array_cache(cell_dof_ids) + + fill!(x,0.0) for patch = 1:length(Ph.patch_decomposition.patch_cells) current_patch_cells = getindex!(cache_patch_cells, Ph.patch_decomposition.patch_cells, @@ -306,17 +307,25 @@ function compute_weight_operators(Ph::PatchFESpace) e = Ph.patch_cell_dofs_ids.ptrs[cell_mesh_overlapped+1]-1 current_patch_cell_dof_ids = view(Ph.patch_cell_dofs_ids.data,s:e) for (dof,pdof) in zip(current_cell_dof_ids,current_patch_cell_dof_ids) - if pdof > 0 && !(dof in keys(touched)) + if pdof > 0 && !(dof ∈ keys(touched)) touched[dof] = true - w[dof] += 1.0 + x[dof] += y[pdof] * w[pdof] / w_sums[dof] end end - cell_mesh_overlapped+=1 + cell_mesh_overlapped += 1 end empty!(touched) end - w .= 1.0 ./ w - w_Ph = similar(w,num_free_dofs(Ph)) - prolongate!(w_Ph,Ph,w) - return w_Ph +end + +function compute_weight_operators(Ph::PatchFESpace) + w = Fill(1.0,num_free_dofs(Ph)) + w_sums = compute_partial_sums(Ph,w) + return w, w_sums +end + +function compute_partial_sums(Ph::PatchFESpace,x) + x_sums = zeros(num_free_dofs(Ph.Vh)) + inject!(x_sums,Ph,x,Fill(1.0,num_free_dofs(Ph))) + return x_sums end diff --git a/test/seq/DistributedPatchFESpacesTests.jl b/test/seq/DistributedPatchFESpacesTests.jl new file mode 100644 index 00000000..ee20d3f9 --- /dev/null +++ b/test/seq/DistributedPatchFESpacesTests.jl @@ -0,0 +1,48 @@ +module DistributedPatchFESpacesTests + +using Test +using Gridap +using GridapDistributed +using PartitionedArrays +using FillArrays + +include("../../src/PatchBasedSmoothers/PatchBasedSmoothers.jl") +import .PatchBasedSmoothers as PBS + +backend = SequentialBackend() +ranks = (1,2) +parts = get_part_ids(backend,ranks) + +domain = (0.0,1.0,0.0,1.0) +partition = (2,4) +model = CartesianDiscreteModel(parts,domain,partition) + +order = 1 +reffe = ReferenceFE(lagrangian,Float64,order) +Vh = TestFESpace(model,reffe) +PD = PBS.PatchDecomposition(model,patch_boundary_style=PBS.PatchBoundaryInclude()) +Ph = PBS.PatchFESpace(model,reffe,Gridap.ReferenceFEs.H1Conformity(),PD,Vh) + +w, w_sums = PBS.compute_weight_operators(Ph,Vh); + +xP = PVector(1.0,Ph.gids) +yP = PVector(0.0,Ph.gids) +x = PVector(1.0,Vh.gids) +y = PVector(0.0,Vh.gids) + +PBS.prolongate!(yP,Ph,x) +PBS.inject!(y,Ph,yP,w,w_sums) + + +assembler = SparseMatrixAssembler(Ph,Ph) +Ωₚ = Triangulation(PD) +dΩₚ = Measure(Ωₚ,2*order+1) +a(u,v) = ∫(∇(v)⋅∇(u))*dΩₚ +l(v) = ∫(1*v)*dΩₚ + +Ah = assemble_matrix(a,assembler,Ph,Ph) +fh = assemble_vector(l,assembler,Ph) + + + +end \ No newline at end of file From 13003bf79590160308ce9c18ef00b3e631fcfa2a Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 15 Feb 2023 17:49:43 +1100 Subject: [PATCH 03/14] Adapted the PatchBasedLinearSolvers --- .../mpi/PatchDecompositions.jl | 2 +- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 2 +- .../seq/PatchBasedLinearSolvers.jl | 34 +++++++++++-------- src/PatchBasedSmoothers/seq/PatchFESpaces.jl | 4 +-- test/seq/DistributedPatchFESpacesTests.jl | 26 +++++++++++--- 5 files changed, 44 insertions(+), 24 deletions(-) diff --git a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl index c208ec4e..5075c697 100644 --- a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl @@ -25,7 +25,7 @@ function Gridap.Geometry.Triangulation(a::DistributedPatchDecomposition) end function get_patch_root_dim(a::DistributedPatchDecomposition) - patch_root_dim = 0 + patch_root_dim = -1 map_parts(a.patch_decompositions) do patch_decomposition patch_root_dim = patch_decomposition.Dr end diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index be1f53ab..972a42f4 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -74,7 +74,7 @@ function compute_weight_operators(Ph::GridapDistributed.DistributedSingleFieldFE w = PVector(0.0,Ph.gids) w_sums = PVector(0.0,Vh.gids) map_parts(w.values,w_sums.values,Ph.spaces) do w, w_sums, Ph - _w, _w_sums = compute_weight_operators(Ph) + _w, _w_sums = compute_weight_operators(Ph,Ph.Vh) w .= _w w_sums .= _w_sums end diff --git a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl index a9d15761..1289b1b3 100644 --- a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl +++ b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl @@ -12,9 +12,10 @@ # (not 100% sure, to investigate) -struct PatchBasedLinearSolver{A} <: Gridap.Algebra.LinearSolver +struct PatchBasedLinearSolver{A,B} <: Gridap.Algebra.LinearSolver bilinear_form :: Function Ph :: A + Vh :: B M :: Gridap.Algebra.LinearSolver end @@ -23,29 +24,30 @@ struct PatchBasedSymbolicSetup <: Gridap.Algebra.SymbolicSetup end function Gridap.Algebra.symbolic_setup(ls::PatchBasedLinearSolver,mat::AbstractMatrix) - PatchBasedSymbolicSetup(ls) + return PatchBasedSymbolicSetup(ls) end -struct PatchBasedSmootherNumericalSetup{A,B,C,D,E} <: Gridap.Algebra.NumericalSetup +struct PatchBasedSmootherNumericalSetup{A,B,C,D,E,F} <: Gridap.Algebra.NumericalSetup solver :: PatchBasedLinearSolver Ap :: A nsAp :: B rp :: C dxp :: D w :: E + w_sums :: F end function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix) - Ph = ss.solver.Ph + Ph, Vh = ss.solver.Ph, ss.solver.Vh assembler = SparseMatrixAssembler(Ph,Ph) - Ap = assemble_matrix(ss.solver.bilinear_form,assembler,Ph,Ph) - solver = ss.solver.M - ssAp = symbolic_setup(solver,Ap) - nsAp = numerical_setup(ssAp,Ap) - rp = _allocate_row_vector(Ap) - dxp = _allocate_col_vector(Ap) - w = compute_weight_operators(Ph) - PatchBasedSmootherNumericalSetup(ss.solver,Ap,nsAp,rp,dxp,w) + Ap = assemble_matrix(ss.solver.bilinear_form,assembler,Ph,Ph) + solver = ss.solver.M + ssAp = symbolic_setup(solver,Ap) + nsAp = numerical_setup(ssAp,Ap) + rp = _allocate_row_vector(Ap) + dxp = _allocate_col_vector(Ap) + w, w_sums = compute_weight_operators(Ph,Vh) + return PatchBasedSmootherNumericalSetup(ss.solver,Ap,nsAp,rp,dxp,w,w_sums) end function _allocate_col_vector(A::AbstractMatrix) @@ -69,9 +71,11 @@ function Gridap.Algebra.numerical_setup!(ns::PatchBasedSmootherNumericalSetup, A end function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumericalSetup,r::AbstractVector) - Ap, nsAp, rp, dxp, w = ns.Ap, ns.nsAp, ns.rp, ns.dxp, ns.w + Ap, nsAp, rp, dxp, w, w_sums = ns.Ap, ns.nsAp, ns.rp, ns.dxp, ns.w, ns.w_sums + Ph = ns.solver.Ph - prolongate!(rp,ns.solver.Ph,r) + prolongate!(rp,Ph,r) solve!(dxp,nsAp,rp) - inject!(x,ns.solver.Ph,dxp,w) + inject!(x,Ph,dxp,w,w_sums) + return x end diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl index 852e4c22..36fc2395 100644 --- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl @@ -256,7 +256,7 @@ end # x \in SingleFESpace # y \in PatchFESpace function inject!(x,Ph::PatchFESpace,y) - w, w_sums = compute_weight_operators(Ph) + w, w_sums = compute_weight_operators(Ph,Ph.Vh) inject!(x,Ph::PatchFESpace,y,w,w_sums) end @@ -318,7 +318,7 @@ function inject!(x,Ph::PatchFESpace,y,w,w_sums) end end -function compute_weight_operators(Ph::PatchFESpace) +function compute_weight_operators(Ph::PatchFESpace,Vh) w = Fill(1.0,num_free_dofs(Ph)) w_sums = compute_partial_sums(Ph,w) return w, w_sums diff --git a/test/seq/DistributedPatchFESpacesTests.jl b/test/seq/DistributedPatchFESpacesTests.jl index ee20d3f9..5b1ee651 100644 --- a/test/seq/DistributedPatchFESpacesTests.jl +++ b/test/seq/DistributedPatchFESpacesTests.jl @@ -1,14 +1,22 @@ module DistributedPatchFESpacesTests +using LinearAlgebra using Test +using PartitionedArrays using Gridap +using Gridap.Helpers +using Gridap.Geometry using GridapDistributed -using PartitionedArrays using FillArrays include("../../src/PatchBasedSmoothers/PatchBasedSmoothers.jl") import .PatchBasedSmoothers as PBS +# This is needed for assembly +include("../../src/MultilevelTools/GridapFixes.jl") + +include("../../src/LinearSolvers/RichardsonSmoothers.jl") + backend = SequentialBackend() ranks = (1,2) parts = get_part_ids(backend,ranks) @@ -32,17 +40,25 @@ y = PVector(0.0,Vh.gids) PBS.prolongate!(yP,Ph,x) PBS.inject!(y,Ph,yP,w,w_sums) +@test x ≈ y +PBS.inject!(x,Ph,xP,w,w_sums) +PBS.prolongate!(yP,Ph,x) +@test xP ≈ yP -assembler = SparseMatrixAssembler(Ph,Ph) Ωₚ = Triangulation(PD) dΩₚ = Measure(Ωₚ,2*order+1) a(u,v) = ∫(∇(v)⋅∇(u))*dΩₚ l(v) = ∫(1*v)*dΩₚ -Ah = assemble_matrix(a,assembler,Ph,Ph) -fh = assemble_vector(l,assembler,Ph) - +assembler = SparseMatrixAssembler(Vh,Vh) +Ah = assemble_matrix(a,assembler,Vh,Vh) +fh = assemble_vector(l,assembler,Vh) +M = PBS.PatchBasedLinearSolver(a,Ph,Vh,LUSolver()) +s = RichardsonSmoother(M,10,1.0/3.0) +x = PBS._allocate_col_vector(Ah) +r = fh-Ah*x +solve!(x,s,Ah,r) end \ No newline at end of file From 2e4f9e6c72018f11709006a282a73a85144bda53 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 16 Feb 2023 17:04:43 +1100 Subject: [PATCH 04/14] Modifications to solver --- Project.toml | 11 ++-- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 2 + .../seq/PatchBasedLinearSolvers.jl | 66 +++++++++++++++---- test/seq/DistributedPatchFESpacesTests.jl | 40 ++++++++++- test/seq/PatchLinearSolverTests.jl | 8 +-- 5 files changed, 101 insertions(+), 26 deletions(-) diff --git a/Project.toml b/Project.toml index 50de03a2..62e5ddec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,8 +1,12 @@ +authors = ["Santiago Badia ", "Jordi Manyer ", "Alberto F. Martin ", "Javier Principe "] name = "GridapSolvers" uuid = "6d3209ee-5e3c-4db7-a716-942eb12ed534" -authors = ["Santiago Badia ", "Jordi Manyer ", "Alberto F. Martin ", "Javier Principe "] version = "0.1.0" +[compat] +PartitionedArrays = "0.2.15" +julia = "1.7" + [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" @@ -16,11 +20,8 @@ MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" -[compat] -PartitionedArrays = "0.2.15" -julia = "1.7" - [extras] +MPIPreferences = "3da0fdf6-3ccc-4f1b-acd9-58baa6c99267" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index 972a42f4..cdd0c40d 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -51,6 +51,7 @@ function prolongate!(x::PVector, map_parts(x.values,Ph.spaces,y.values) do x,Ph,y prolongate!(x,Ph,y) end + exchange!(x) end function inject!(x::PVector, @@ -59,6 +60,7 @@ function inject!(x::PVector, w::PVector, w_sums::PVector) + exchange!(y) map_parts(x.values,Ph.spaces,y.values,w.values,w_sums.values) do x,Ph,y,w,w_sums inject!(x,Ph,y,w,w_sums) end diff --git a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl index 1289b1b3..3758f434 100644 --- a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl +++ b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl @@ -27,27 +27,45 @@ function Gridap.Algebra.symbolic_setup(ls::PatchBasedLinearSolver,mat::AbstractM return PatchBasedSymbolicSetup(ls) end -struct PatchBasedSmootherNumericalSetup{A,B,C,D,E,F} <: Gridap.Algebra.NumericalSetup - solver :: PatchBasedLinearSolver - Ap :: A - nsAp :: B - rp :: C - dxp :: D - w :: E - w_sums :: F +struct PatchBasedSmootherNumericalSetup{A,B,C,D} <: Gridap.Algebra.NumericalSetup + solver :: PatchBasedLinearSolver + Ap :: A + Ap_ns :: B + weights :: C + caches :: D end function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix) Ph, Vh = ss.solver.Ph, ss.solver.Vh + weights = compute_weight_operators(Ph,Vh) + + # Assemble patch system assembler = SparseMatrixAssembler(Ph,Ph) Ap = assemble_matrix(ss.solver.bilinear_form,assembler,Ph,Ph) - solver = ss.solver.M - ssAp = symbolic_setup(solver,Ap) - nsAp = numerical_setup(ssAp,Ap) + + # Patch system solver + Ap_solver = ss.solver.M + Ap_ss = symbolic_setup(Ap_solver,Ap) + Ap_ns = numerical_setup(Ap_ss,Ap) + + # Caches + caches = _patch_based_solver_caches(Ph,Ap) + + return PatchBasedSmootherNumericalSetup(ss.solver,Ap,Ap_ns,weights,caches) +end + +function _patch_based_solver_caches(Ph::PatchFESpace,Ap::AbstractMatrix) rp = _allocate_row_vector(Ap) dxp = _allocate_col_vector(Ap) - w, w_sums = compute_weight_operators(Ph,Vh) - return PatchBasedSmootherNumericalSetup(ss.solver,Ap,nsAp,rp,dxp,w,w_sums) + return rp, dxp +end + +function _patch_based_solver_caches(Ph::GridapDistributed.DistributedSingleFieldFESpace,Ap::PSparseMatrix) + rp_mat = _allocate_row_vector(Ap) + dxp_mat = _allocate_col_vector(Ap) + rp = PVector(0.0,Ph.gids) + dxp = PVector(0.0,Ph.gids) + return rp_mat, dxp_mat, rp, dxp end function _allocate_col_vector(A::AbstractMatrix) @@ -71,11 +89,31 @@ function Gridap.Algebra.numerical_setup!(ns::PatchBasedSmootherNumericalSetup, A end function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumericalSetup,r::AbstractVector) - Ap, nsAp, rp, dxp, w, w_sums = ns.Ap, ns.nsAp, ns.rp, ns.dxp, ns.w, ns.w_sums + Ap_ns, weights, caches = ns.Ap_ns, ns.weights, ns.caches + Ph = ns.solver.Ph + w, w_sums = weights + rp, dxp = caches prolongate!(rp,Ph,r) solve!(dxp,nsAp,rp) inject!(x,Ph,dxp,w,w_sums) + + return x +end + +function Gridap.Algebra.solve!(x::PVector,ns::PatchBasedSmootherNumericalSetup,r::PVector) + Ap_ns, weights, caches = ns.Ap_ns, ns.weights, ns.caches + + Ph = ns.solver.Ph + w, w_sums = weights + rp_mat, dxp_mat, rp, dxp = caches + + prolongate!(rp,Ph,r) + copy!(rp_mat,rp) + solve!(dxp_mat,Ap_ns,rp_mat) + copy!(dxp,dxp_mat) + inject!(x,Ph,dxp,w,w_sums) + return x end diff --git a/test/seq/DistributedPatchFESpacesTests.jl b/test/seq/DistributedPatchFESpacesTests.jl index 5b1ee651..c9ee8232 100644 --- a/test/seq/DistributedPatchFESpacesTests.jl +++ b/test/seq/DistributedPatchFESpacesTests.jl @@ -1,5 +1,8 @@ module DistributedPatchFESpacesTests +ENV["JULIA_MPI_BINARY"] = "system" +ENV["JULIA_MPI_PATH"] = "/usr/lib/x86_64-linux-gnu" + using LinearAlgebra using Test using PartitionedArrays @@ -28,7 +31,7 @@ model = CartesianDiscreteModel(parts,domain,partition) order = 1 reffe = ReferenceFE(lagrangian,Float64,order) Vh = TestFESpace(model,reffe) -PD = PBS.PatchDecomposition(model,patch_boundary_style=PBS.PatchBoundaryInclude()) +PD = PBS.PatchDecomposition(model)#,patch_boundary_style=PBS.PatchBoundaryInclude()) Ph = PBS.PatchFESpace(model,reffe,Gridap.ReferenceFEs.H1Conformity(),PD,Vh) w, w_sums = PBS.compute_weight_operators(Ph,Vh); @@ -56,9 +59,40 @@ Ah = assemble_matrix(a,assembler,Vh,Vh) fh = assemble_vector(l,assembler,Vh) M = PBS.PatchBasedLinearSolver(a,Ph,Vh,LUSolver()) -s = RichardsonSmoother(M,10,1.0/3.0) +R = RichardsonSmoother(M,10,1.0/3.0) +Rss = symbolic_setup(R,Ah) +Rns = numerical_setup(Rss,Ah) + x = PBS._allocate_col_vector(Ah) r = fh-Ah*x -solve!(x,s,Ah,r) +exchange!(r) +solve!(x,Rns,r) + +Mss = symbolic_setup(M,Ah) +Mns = numerical_setup(Mss,Ah) +solve!(x,Mns,r) + +assembler_P = SparseMatrixAssembler(Ph,Ph) +Ahp = assemble_matrix(a,assembler_P,Ph,Ph) +fhp = assemble_vector(l,assembler_P,Ph) + +lu = LUSolver() +luss = symbolic_setup(lu,Ahp) +luns = numerical_setup(luss,Ahp) + +rp = PVector(0.0,Ph.gids) +PBS.prolongate!(rp,Ph,r) + +rp_mat = PVector(0.0,Ahp.cols) +copy!(rp_mat,rp) +xp_mat = PVector(0.0,Ahp.cols) + +solve!(xp_mat,luns,rp_mat) + +xp = PVector(0.0,Ph.gids) +copy!(xp,xp_mat) + +w, w_sums = PBS.compute_weight_operators(Ph,Vh); +PBS.inject!(x,Ph,xp,w,w_sums) end \ No newline at end of file diff --git a/test/seq/PatchLinearSolverTests.jl b/test/seq/PatchLinearSolverTests.jl index 9e38869f..015e5134 100644 --- a/test/seq/PatchLinearSolverTests.jl +++ b/test/seq/PatchLinearSolverTests.jl @@ -10,13 +10,13 @@ module PatchLinearSolverTests using GridapSolvers using GridapSolvers.PatchBasedSmoothers - const order=1 + order=1 - function returns_PD_Ph_xh_Vh(model) + function returns_PD_Ph_xh_Vh(model;style=GridapSolvers.PatchBasedSmoothers.PatchBoundaryExclude()) reffe = ReferenceFE(lagrangian,Float64,order) # reffe=ReferenceFE(lagrangian,VectorValue{2,Float64},order) @santiagobadia: For Vector Laplacian Vh = TestFESpace(model,reffe) - PD = PatchDecomposition(model) + PD = PatchDecomposition(model;patch_boundary_style=style) Ph = PatchFESpace(model,reffe,H1Conformity(),PD,Vh) assembler = SparseMatrixAssembler(Ph,Ph) Ωₚ = Triangulation(PD) @@ -51,7 +51,7 @@ module PatchLinearSolverTests dΩₚ = Measure(Ωₚ,2*order+1) a(u,v) = ∫(∇(v)⋅∇(u))*dΩₚ # α =1,0; a(u,v)=∫(v⋅u)dΩ+∫(α*∇(v)⊙∇(u))dΩ # @santiagobadia: For vector Laplacian - M = PatchBasedLinearSolver(a,Ph,LUSolver()) + M = PatchBasedLinearSolver(a,Ph,Vh,LUSolver()) s = RichardsonSmoother(M,10,1.0/3.0) x = GridapSolvers.PatchBasedSmoothers._allocate_col_vector(A) r = b-A*x From bd6d8cca1228d1f5b83be8568d707ed51f335d9e Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Wed, 22 Feb 2023 15:04:08 +1100 Subject: [PATCH 05/14] Fixed distributed Patch-Based smoothers --- .../seq/PatchBasedLinearSolvers.jl | 21 +++-- test/seq/DistributedPatchFESpacesTests.jl | 90 +++++++++++-------- 2 files changed, 69 insertions(+), 42 deletions(-) diff --git a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl index 3758f434..9178178d 100644 --- a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl +++ b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl @@ -49,23 +49,27 @@ function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractM Ap_ns = numerical_setup(Ap_ss,Ap) # Caches - caches = _patch_based_solver_caches(Ph,Ap) + caches = _patch_based_solver_caches(Ph,Vh,Ap) return PatchBasedSmootherNumericalSetup(ss.solver,Ap,Ap_ns,weights,caches) end -function _patch_based_solver_caches(Ph::PatchFESpace,Ap::AbstractMatrix) +function _patch_based_solver_caches(Ph::PatchFESpace,Vh::FESpace,Ap::AbstractMatrix) rp = _allocate_row_vector(Ap) dxp = _allocate_col_vector(Ap) return rp, dxp end -function _patch_based_solver_caches(Ph::GridapDistributed.DistributedSingleFieldFESpace,Ap::PSparseMatrix) +function _patch_based_solver_caches(Ph::GridapDistributed.DistributedSingleFieldFESpace, + Vh::GridapDistributed.DistributedSingleFieldFESpace, + Ap::PSparseMatrix) rp_mat = _allocate_row_vector(Ap) dxp_mat = _allocate_col_vector(Ap) rp = PVector(0.0,Ph.gids) dxp = PVector(0.0,Ph.gids) - return rp_mat, dxp_mat, rp, dxp + r = PVector(0.0,Vh.gids) + x = PVector(0.0,Vh.gids) + return rp_mat, dxp_mat, rp, dxp, r, x end function _allocate_col_vector(A::AbstractMatrix) @@ -102,18 +106,21 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumerical return x end -function Gridap.Algebra.solve!(x::PVector,ns::PatchBasedSmootherNumericalSetup,r::PVector) +function Gridap.Algebra.solve!(x_mat::PVector,ns::PatchBasedSmootherNumericalSetup,r_mat::PVector) Ap_ns, weights, caches = ns.Ap_ns, ns.weights, ns.caches Ph = ns.solver.Ph w, w_sums = weights - rp_mat, dxp_mat, rp, dxp = caches + rp_mat, dxp_mat, rp, dxp, r, x = caches + copy!(r,r_mat) + exchange!(r) prolongate!(rp,Ph,r) copy!(rp_mat,rp) solve!(dxp_mat,Ap_ns,rp_mat) copy!(dxp,dxp_mat) inject!(x,Ph,dxp,w,w_sums) + copy!(x_mat,x) - return x + return x_mat end diff --git a/test/seq/DistributedPatchFESpacesTests.jl b/test/seq/DistributedPatchFESpacesTests.jl index c9ee8232..dc0a7dcc 100644 --- a/test/seq/DistributedPatchFESpacesTests.jl +++ b/test/seq/DistributedPatchFESpacesTests.jl @@ -1,8 +1,5 @@ module DistributedPatchFESpacesTests -ENV["JULIA_MPI_BINARY"] = "system" -ENV["JULIA_MPI_PATH"] = "/usr/lib/x86_64-linux-gnu" - using LinearAlgebra using Test using PartitionedArrays @@ -12,13 +9,8 @@ using Gridap.Geometry using GridapDistributed using FillArrays -include("../../src/PatchBasedSmoothers/PatchBasedSmoothers.jl") -import .PatchBasedSmoothers as PBS - -# This is needed for assembly -include("../../src/MultilevelTools/GridapFixes.jl") - -include("../../src/LinearSolvers/RichardsonSmoothers.jl") +using GridapSolvers +import GridapSolvers.PatchBasedSmoothers as PBS backend = SequentialBackend() ranks = (1,2) @@ -34,6 +26,7 @@ Vh = TestFESpace(model,reffe) PD = PBS.PatchDecomposition(model)#,patch_boundary_style=PBS.PatchBoundaryInclude()) Ph = PBS.PatchFESpace(model,reffe,Gridap.ReferenceFEs.H1Conformity(),PD,Vh) +# ---- Testing Prolongation and Injection ---- # w, w_sums = PBS.compute_weight_operators(Ph,Vh); xP = PVector(1.0,Ph.gids) @@ -49,50 +42,77 @@ PBS.inject!(x,Ph,xP,w,w_sums) PBS.prolongate!(yP,Ph,x) @test xP ≈ yP -Ωₚ = Triangulation(PD) -dΩₚ = Measure(Ωₚ,2*order+1) -a(u,v) = ∫(∇(v)⋅∇(u))*dΩₚ -l(v) = ∫(1*v)*dΩₚ +# ---- Assemble systems ---- # +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order+1) +a(u,v) = ∫(v⋅u)*dΩ +l(v) = ∫(1*v)*dΩ assembler = SparseMatrixAssembler(Vh,Vh) Ah = assemble_matrix(a,assembler,Vh,Vh) fh = assemble_vector(l,assembler,Vh) -M = PBS.PatchBasedLinearSolver(a,Ph,Vh,LUSolver()) -R = RichardsonSmoother(M,10,1.0/3.0) -Rss = symbolic_setup(R,Ah) -Rns = numerical_setup(Rss,Ah) +sol_h = solve(LUSolver(),Ah,fh) -x = PBS._allocate_col_vector(Ah) -r = fh-Ah*x -exchange!(r) -solve!(x,Rns,r) +Ωₚ = Triangulation(PD) +dΩₚ = Measure(Ωₚ,2*order+1) +ap(u,v) = ∫(v⋅u)*dΩₚ +lp(v) = ∫(1*v)*dΩₚ + +assembler_P = SparseMatrixAssembler(Ph,Ph) +Ahp = assemble_matrix(ap,assembler_P,Ph,Ph) +fhp = assemble_vector(lp,assembler_P,Ph) +# ---- Define solvers ---- # +LU = LUSolver() +LUss = symbolic_setup(LU,Ahp) +LUns = numerical_setup(LUss,Ahp) + +M = PBS.PatchBasedLinearSolver(ap,Ph,Vh,LU) Mss = symbolic_setup(M,Ah) Mns = numerical_setup(Mss,Ah) -solve!(x,Mns,r) -assembler_P = SparseMatrixAssembler(Ph,Ph) -Ahp = assemble_matrix(a,assembler_P,Ph,Ph) -fhp = assemble_vector(l,assembler_P,Ph) +R = RichardsonSmoother(M,10,1.0/3.0) +Rss = symbolic_setup(R,Ah) +Rns = numerical_setup(Rss,Ah) -lu = LUSolver() -luss = symbolic_setup(lu,Ahp) -luns = numerical_setup(luss,Ahp) +# ---- Manual solve using LU ---- # -rp = PVector(0.0,Ph.gids) -PBS.prolongate!(rp,Ph,r) +x1_mat = PVector(0.5,Ah.cols) +r1_mat = fh-Ah*x1_mat +exchange!(r1_mat) +r1 = PVector(0.0,Vh.gids) +x1 = PVector(0.0,Vh.gids) +rp = PVector(0.0,Ph.gids) +xp = PVector(0.0,Ph.gids) rp_mat = PVector(0.0,Ahp.cols) -copy!(rp_mat,rp) xp_mat = PVector(0.0,Ahp.cols) -solve!(xp_mat,luns,rp_mat) +copy!(r1,r1_mat) +exchange!(r1) +PBS.prolongate!(rp,Ph,r1) -xp = PVector(0.0,Ph.gids) +copy!(rp_mat,rp) +solve!(xp_mat,LUns,rp_mat) copy!(xp,xp_mat) w, w_sums = PBS.compute_weight_operators(Ph,Vh); -PBS.inject!(x,Ph,xp,w,w_sums) +PBS.inject!(x1,Ph,xp,w,w_sums) +copy!(x1_mat,x1) + +# ---- Same using the PatchBasedSmoother ---- # + +x2_mat = PVector(0.5,Ah.cols) +r2_mat = fh-Ah*x2_mat +exchange!(r2_mat) +solve!(x2_mat,Mns,r2_mat) + +# ---- Smoother inside Richardson +x3_mat = PVector(0.5,Ah.cols) +r3_mat = fh-Ah*x3_mat +exchange!(r3_mat) +solve!(x3_mat,Rns,r3_mat) +exchange!(x3_mat) end \ No newline at end of file From 7b73cf44f57854ee13ca66f42f51842b43fd8737 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 23 Feb 2023 14:33:01 +1100 Subject: [PATCH 06/14] Bugfix: TransferOps not working for VectorValued FEFunctions --- .../DistributedGridTransferOperators.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/MultilevelTools/DistributedGridTransferOperators.jl b/src/MultilevelTools/DistributedGridTransferOperators.jl index 9521e614..d991bed6 100644 --- a/src/MultilevelTools/DistributedGridTransferOperators.jl +++ b/src/MultilevelTools/DistributedGridTransferOperators.jl @@ -79,7 +79,7 @@ function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode:: model_h = get_model_before_redist(mh,lev) Uh = get_fe_space_before_redist(sh,lev) Ωh = Triangulation(model_h) - fv_h = PVector(0.0,Uh.gids) + fv_h = zero_free_values(Uh) dv_h = (mode == :solution) ? get_dirichlet_dof_values(Uh) : zero_dirichlet_values(Uh) model_H = get_model(mh,lev+1) @@ -92,10 +92,15 @@ function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode:: aH(u,v) = ∫(v⋅u)*dΩH lH(v,uh) = ∫(v⋅uh)*dΩhH assem = SparseMatrixAssembler(UH,VH) - - u_dir = (mode == :solution) ? interpolate(0.0,UH) : interpolate_everywhere(0.0,UH) + + fv_H = zero_free_values(UH) + dv_H = zero_dirichlet_values(UH) + u0 = FEFunction(UH,fv_H,true) # Zero at free dofs + u00 = FEFunction(UH,fv_H,dv_H,true) # Zero everywhere + + u_dir = (mode == :solution) ? u0 : u00 u,v = get_trial_fe_basis(UH), get_fe_basis(VH) - data = collect_cell_matrix_and_vector(UH,VH,aH(u,v),lH(v,0.0),u_dir) + data = collect_cell_matrix_and_vector(UH,VH,aH(u,v),lH(v,u00),u_dir) AH,bH0 = assemble_matrix_and_vector(assem,data) xH = PVector(0.0,AH.cols) bH = copy(bH0) From c64b505c5ea91a94b5ae44d400140625a620e98f Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 23 Feb 2023 14:33:36 +1100 Subject: [PATCH 07/14] Added PatchBased methods for multilevel --- .../PatchBasedSmoothers.jl | 2 ++ .../mpi/PatchDecompositions.jl | 23 +++++++++++--- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 30 +++++++++++++++---- 3 files changed, 46 insertions(+), 9 deletions(-) diff --git a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl index dcaf0d6c..e6a03299 100644 --- a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl +++ b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl @@ -12,6 +12,8 @@ using Gridap.FESpaces using PartitionedArrays using GridapDistributed +using GridapSolvers.MultilevelTools + export PatchDecomposition export PatchFESpace export PatchBasedLinearSolver diff --git a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl index 5075c697..79a44488 100644 --- a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl @@ -4,24 +4,39 @@ struct DistributedPatchDecomposition{Dc,Dp,A,B} <: GridapType model::B end -function PatchDecomposition(model::GridapDistributed.DistributedDiscreteModel{Dc,Dp}; +GridapDistributed.local_views(a::DistributedPatchDecomposition) = a.patch_decompositions + +function PatchDecomposition(model::GridapDistributed.AbstractDistributedDiscreteModel{Dc,Dp}; Dr=0, patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude()) where {Dc,Dp} - patch_decompositions = map_parts(model.models) do lmodel + patch_decompositions = map_parts(local_views(model)) do lmodel PatchDecomposition(lmodel; Dr=Dr, patch_boundary_style=patch_boundary_style) end A = typeof(patch_decompositions) B = typeof(model) - DistributedPatchDecomposition{Dc,Dp,A,B}(patch_decompositions,model) + return DistributedPatchDecomposition{Dc,Dp,A,B}(patch_decompositions,model) +end + +function PatchDecomposition(mh::ModelHierarchy;kwargs...) + nlevs = num_levels(mh) + decompositions = Vector{DistributedPatchDecomposition}(undef,nlevs-1) + for lev in 1:nlevs-1 + parts = get_level_parts(mh,lev) + if i_am_in(parts) + model = get_model(mh,lev) + decompositions[lev] = PatchDecomposition(model;kwargs...) + end + end + return decompositions end function Gridap.Geometry.Triangulation(a::DistributedPatchDecomposition) trians = map_parts(a.patch_decompositions) do a Triangulation(a) end - GridapDistributed.DistributedTriangulation(trians,a.model) + return GridapDistributed.DistributedTriangulation(trians,a.model) end function get_patch_root_dim(a::DistributedPatchDecomposition) diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index cdd0c40d..2278eb52 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -7,7 +7,7 @@ # 3. Each processor needs to know how many patches "touch" its owned DoFs. # This requires NO->O communication as well. [PENDING] -function PatchFESpace(model::GridapDistributed.DistributedDiscreteModel, +function PatchFESpace(model::GridapDistributed.AbstractDistributedDiscreteModel, reffe::Tuple{<:Gridap.FESpaces.ReferenceFEName,Any,Any}, conformity::Gridap.FESpaces.Conformity, patch_decomposition::DistributedPatchDecomposition, @@ -26,12 +26,12 @@ function PatchFESpace(model::GridapDistributed.DistributedDiscreteModel, end spaces = map_parts(f, - model.models, - patch_decomposition.patch_decompositions, - Vh.spaces, + local_views(model), + local_views(patch_decomposition), + local_views(Vh), root_gids.partition) - parts = get_part_ids(model.models) + parts = get_parts(model) nodofs = map_parts(spaces) do space num_free_dofs(space) end @@ -43,6 +43,26 @@ function PatchFESpace(model::GridapDistributed.DistributedDiscreteModel, return GridapDistributed.DistributedSingleFieldFESpace(spaces,gids,get_vector_type(Vh)) end +function PatchFESpace(mh::ModelHierarchy, + reffe::Tuple{<:Gridap.FESpaces.ReferenceFEName,Any,Any}, + conformity::Gridap.FESpaces.Conformity, + patch_decompositions::AbstractArray{<:DistributedPatchDecomposition}, + sh::FESpaceHierarchy) + nlevs = num_levels(mh) + levels = Vector{MultilevelTools.FESpaceHierarchyLevel}(undef,nlevs) + for lev in 1:nlevs-1 + parts = get_level_parts(mh,lev) + if i_am_in(parts) + model = get_model(mh,lev) + space = MultilevelTools.get_fe_space(sh,lev) + decomp = patch_decompositions[lev] + patch_space = PatchFESpace(model,reffe,conformity,decomp,space) + levels[lev] = MultilevelTools.FESpaceHierarchyLevel(lev,nothing,patch_space) + end + end + return FESpaceHierarchy(mh,levels) +end + # x \in PatchFESpace # y \in SingleFESpace function prolongate!(x::PVector, From 17b2c132094b172ffec3806891cca4e5d5bbfad8 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Thu, 23 Feb 2023 14:34:30 +1100 Subject: [PATCH 08/14] Adapted HDiv test to new patch based solvers --- test/mpi/GMGLinearSolversHDivRTTests.jl | 34 ++++++++++++++++++++++--- 1 file changed, 30 insertions(+), 4 deletions(-) diff --git a/test/mpi/GMGLinearSolversHDivRTTests.jl b/test/mpi/GMGLinearSolversHDivRTTests.jl index 1a9bff6f..f9513dd7 100644 --- a/test/mpi/GMGLinearSolversHDivRTTests.jl +++ b/test/mpi/GMGLinearSolversHDivRTTests.jl @@ -13,11 +13,33 @@ using GridapP4est using GridapSolvers using GridapSolvers.LinearSolvers +using GridapSolvers.MultilevelTools +using GridapSolvers.PatchBasedSmoothers u(x) = VectorValue(x[1],x[2]) f(x) = VectorValue(2.0*x[2]*(1.0-x[1]*x[1]),2.0*x[1]*(1-x[2]*x[2])) +function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree) + mh = tests.mh + nlevs = num_levels(mh) + smoothers = Vector{RichardsonSmoother}(undef,nlevs-1) + for lev in 1:nlevs-1 + parts = get_level_parts(mh,lev) + if i_am_in(parts) + PD = patch_decompositions[lev] + Ph = get_fe_space(patch_spaces,lev) + Vh = get_fe_space(tests,lev) + Ω = Triangulation(PD) + dΩ = Measure(Ω,qdegree) + a(u,v) = biform(u,v,dΩ) + patch_smoother = PatchBasedLinearSolver(a,Ph,Vh,LUSolver()) + smoothers[lev] = RichardsonSmoother(patch_smoother,1,1.0/3.0) + end + end + return smoothers +end + function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, order, α) GridapP4est.with(parts) do domain = (0,1,0,1) @@ -32,12 +54,15 @@ function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, tests = TestFESpace(mh,reffe;dirichlet_tags="boundary") trials = TrialFESpace(tests,u) + patch_decompositions = PatchDecomposition(mh) + patch_spaces = PatchFESpace(mh,reffe,DivConformity(),patch_decompositions,tests) + biform(u,v,dΩ) = ∫(v⋅u)dΩ + ∫(α*divergence(v)⋅divergence(u))dΩ liform(v,dΩ) = ∫(v⋅f)dΩ smatrices, A, b = compute_hierarchy_matrices(trials,biform,liform,qdegree) # Preconditioner - smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0),num_levels-1) + smoothers = get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree) restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual) gmg = GMGLinearSolver(mh, @@ -57,9 +82,10 @@ function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, x = PVector(0.0,A.cols) x, history = IterativeSolvers.cg!(x,A,b; verbose=i_am_main(parts), - reltol=1.0e-12, + reltol=1.0e-8, Pl=ns, - log=true) + log=true, + maxiter=10) # Error norms and print solution model = get_model(mh,1) @@ -93,7 +119,7 @@ num_refs_coarse = 2 α = 1.0 num_parts_x_level = [4,2,1] ranks = num_parts_x_level[1] -#num_iters, num_free_dofs2 = with_backend(main,MPIBackend(),ranks,coarse_grid_partition,num_parts_x_level,num_refs_coarse,order,α) +num_iters, num_free_dofs2 = with_backend(main,MPIBackend(),ranks,coarse_grid_partition,num_parts_x_level,num_refs_coarse,order,α) """ From 31f58b96566aceee40993a3e8086894e32806068 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 28 Feb 2023 12:23:11 +1100 Subject: [PATCH 09/14] Fix Project.toml --- Project.toml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/Project.toml b/Project.toml index 02bd1f70..c4c8c073 100644 --- a/Project.toml +++ b/Project.toml @@ -3,10 +3,6 @@ name = "GridapSolvers" uuid = "6d3209ee-5e3c-4db7-a716-942eb12ed534" version = "0.1.0" -[compat] -PartitionedArrays = "0.2.15" -julia = "1.7" - [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" From 0fb4b121ae4a71a11682a98ae7d73e5fe1c405da Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 28 Feb 2023 12:26:16 +1100 Subject: [PATCH 10/14] Updated to GridapDistributed 0.2.7 --- src/PatchBasedSmoothers/mpi/PatchDecompositions.jl | 2 +- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl index 79a44488..20b45e7c 100644 --- a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl @@ -6,7 +6,7 @@ end GridapDistributed.local_views(a::DistributedPatchDecomposition) = a.patch_decompositions -function PatchDecomposition(model::GridapDistributed.AbstractDistributedDiscreteModel{Dc,Dp}; +function PatchDecomposition(model::GridapDistributed.DistributedDiscreteModel{Dc,Dp}; Dr=0, patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude()) where {Dc,Dp} patch_decompositions = map_parts(local_views(model)) do lmodel diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index 2278eb52..24232d87 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -7,7 +7,7 @@ # 3. Each processor needs to know how many patches "touch" its owned DoFs. # This requires NO->O communication as well. [PENDING] -function PatchFESpace(model::GridapDistributed.AbstractDistributedDiscreteModel, +function PatchFESpace(model::GridapDistributed.DistributedDiscreteModel, reffe::Tuple{<:Gridap.FESpaces.ReferenceFEName,Any,Any}, conformity::Gridap.FESpaces.Conformity, patch_decomposition::DistributedPatchDecomposition, From ffa839aa5834931adeedf4c023ba0d70912cc21a Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 3 Mar 2023 12:44:03 +1100 Subject: [PATCH 11/14] Bugfix: Edge-based patch-based smoothers When deciding whether a facet belongs to the patch boundary, the serial version was not correctly treating interface facets (which are neither marked as boundary nor have any neighboring cell not belonging to the patch). --- .../mpi/PatchDecompositions.jl | 38 ++- src/PatchBasedSmoothers/mpi/PatchFESpaces.jl | 4 +- .../seq/PatchDecompositions.jl | 46 ++-- src/PatchBasedSmoothers/seq/PatchFESpaces.jl | 6 +- test/mpi/GMGLinearSolversHDivRTTests.jl | 4 +- .../DistributedPatchFESpacesDebuggingTests.jl | 223 ++++++++++++++++++ test/seq/DistributedPatchFESpacesTests.jl | 20 +- 7 files changed, 307 insertions(+), 34 deletions(-) create mode 100644 test/seq/DistributedPatchFESpacesDebuggingTests.jl diff --git a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl index 20b45e7c..36c5dcba 100644 --- a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl @@ -9,10 +9,12 @@ GridapDistributed.local_views(a::DistributedPatchDecomposition) = a.patch_decomp function PatchDecomposition(model::GridapDistributed.DistributedDiscreteModel{Dc,Dp}; Dr=0, patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude()) where {Dc,Dp} + mark_interface_facets!(model) patch_decompositions = map_parts(local_views(model)) do lmodel PatchDecomposition(lmodel; Dr=Dr, - patch_boundary_style=patch_boundary_style) + patch_boundary_style=patch_boundary_style, + boundary_tag_names=["boundary","interface"]) end A = typeof(patch_decompositions) B = typeof(model) @@ -46,3 +48,37 @@ function get_patch_root_dim(a::DistributedPatchDecomposition) end return patch_root_dim end + +function mark_interface_facets!(model::GridapDistributed.DistributedDiscreteModel{Dc,Dp}) where {Dc,Dp} + face_labeling = get_face_labeling(model) + topo = get_grid_topology(model) + + map_parts(local_views(face_labeling),local_views(topo)) do face_labeling, topo + tag_to_name = face_labeling.tag_to_name + tag_to_entities = face_labeling.tag_to_entities + d_to_dface_to_entity = face_labeling.d_to_dface_to_entity + + # Create new tag & entity + interface_entity = maximum(map(maximum,tag_to_entities)) + 1 + push!(tag_to_entities,[interface_entity]) + push!(tag_to_name,"interface") + + # Interface faces should also be interior + interior_tag = findfirst(x->(x=="interior"),tag_to_name) + push!(tag_to_entities[interior_tag],interface_entity) + + # Select interface entities + boundary_tag = findfirst(x->(x=="boundary"),tag_to_name) + boundary_entities = tag_to_entities[boundary_tag] + + f2c_map = Geometry.get_faces(topo,Dc-1,Dc) + num_cells_around_facet = map(length,f2c_map) + mx = maximum(num_cells_around_facet) + for (f,nf) in enumerate(num_cells_around_facet) + is_boundary = (d_to_dface_to_entity[Dc][f] ∈ boundary_entities) + if !is_boundary && (nf != mx) + d_to_dface_to_entity[Dc][f] = interface_entity + end + end + end +end diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl index 24232d87..1f6614e9 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl @@ -74,13 +74,15 @@ function prolongate!(x::PVector, exchange!(x) end +# x \in SingleFESpace +# y \in PatchFESpace function inject!(x::PVector, Ph::GridapDistributed.DistributedSingleFieldFESpace, y::PVector, w::PVector, w_sums::PVector) - exchange!(y) + #exchange!(y) map_parts(x.values,Ph.spaces,y.values,w.values,w_sums.values) do x,Ph,y,w,w_sums inject!(x,Ph,y,w,w_sums) end diff --git a/src/PatchBasedSmoothers/seq/PatchDecompositions.jl b/src/PatchBasedSmoothers/seq/PatchDecompositions.jl index df093840..913f0c3d 100644 --- a/src/PatchBasedSmoothers/seq/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/seq/PatchDecompositions.jl @@ -20,7 +20,8 @@ Gridap.Geometry.num_cells(a::PatchDecomposition) = a.patch_cells_overlapped_mesh function PatchDecomposition( model::DiscreteModel{Dc,Dp}; Dr=0, - patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude()) where {Dc,Dp} + patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude(), + boundary_tag_names::AbstractArray{String}=["boundary"]) where {Dc,Dp} Gridap.Helpers.@check 0 <= Dr <= Dc-1 grid = get_grid(model) @@ -41,12 +42,13 @@ function PatchDecomposition( patch_cells, patch_cells_overlapped_mesh) - generate_patch_boundary_faces!(model, - patch_cells_faces_on_boundary, + generate_patch_boundary_faces!(patch_cells_faces_on_boundary, + model, patch_cells, patch_cells_overlapped_mesh, patch_facets, - patch_boundary_style) + patch_boundary_style, + boundary_tag_names) return PatchDecomposition{Dc,Dp}(model, Dr, patch_cells, @@ -112,47 +114,45 @@ function allocate_cell_overlapped_mesh_lface(::Type{T}, return Gridap.Arrays.Table(data,ptrs) end -function generate_patch_boundary_faces!(model, - patch_cells_faces_on_boundary, +function generate_patch_boundary_faces!(patch_cells_faces_on_boundary, + model::DiscreteModel, patch_cells, patch_cells_overlapped_mesh, patch_facets, - patch_boundary_style) - Dc = num_cell_dims(model) - topology = get_grid_topology(model) - labeling = get_face_labeling(model) - num_patches = length(patch_cells.ptrs)-1 + patch_boundary_style, + boundary_tag_names) + num_patches = length(patch_cells.ptrs)-1 cache_patch_cells = array_cache(patch_cells) cache_patch_facets = array_cache(patch_facets) for patch = 1:num_patches current_patch_cells = getindex!(cache_patch_cells,patch_cells,patch) current_patch_facets = getindex!(cache_patch_facets,patch_facets,patch) generate_patch_boundary_faces!(patch_cells_faces_on_boundary, - Dc, - topology, - labeling, + model, patch, current_patch_cells, patch_cells_overlapped_mesh, current_patch_facets, - patch_boundary_style) + patch_boundary_style, + boundary_tag_names) end end function generate_patch_boundary_faces!(patch_cells_faces_on_boundary, - Dc, - topology, - face_labeling, + model::DiscreteModel{Dc}, patch, patch_cells, patch_cells_overlapped_mesh, patch_facets, - patch_boundary_style) + patch_boundary_style, + boundary_tag_names) where Dc + face_labeling = get_face_labeling(model) + topology = get_grid_topology(model) - boundary_tag = findfirst(x->(x=="boundary"),face_labeling.tag_to_name) - Gridap.Helpers.@check !isa(boundary_tag, Nothing) - boundary_entities = face_labeling.tag_to_entities[boundary_tag] + boundary_tags = findall(x -> (x ∈ boundary_tag_names), face_labeling.tag_to_name) + Gridap.Helpers.@check !isempty(boundary_tags) + boundary_entities = vcat(face_labeling.tag_to_entities[boundary_tags]...) # Cells facets Df = Dc-1 @@ -184,7 +184,7 @@ function generate_patch_boundary_faces!(patch_cells_faces_on_boundary, facet_at_global_boundary = (facet_entity ∈ boundary_entities) A = (facet_at_global_boundary) && (facet ∉ patch_facets) B = (patch_boundary_style isa PatchBoundaryExclude) && cell_not_in_patch_found - facet_at_patch_boundary = A || B + facet_at_patch_boundary = (A || B) if (facet_at_patch_boundary) diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl index 36fc2395..9758e813 100644 --- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl @@ -320,12 +320,12 @@ end function compute_weight_operators(Ph::PatchFESpace,Vh) w = Fill(1.0,num_free_dofs(Ph)) - w_sums = compute_partial_sums(Ph,w) + w_sums = compute_partial_sums(Ph,Vh,w) return w, w_sums end -function compute_partial_sums(Ph::PatchFESpace,x) - x_sums = zeros(num_free_dofs(Ph.Vh)) +function compute_partial_sums(Ph::PatchFESpace,Vh,x) + x_sums = zeros(num_free_dofs(Vh)) inject!(x_sums,Ph,x,Fill(1.0,num_free_dofs(Ph))) return x_sums end diff --git a/test/mpi/GMGLinearSolversHDivRTTests.jl b/test/mpi/GMGLinearSolversHDivRTTests.jl index f9513dd7..a982d8da 100644 --- a/test/mpi/GMGLinearSolversHDivRTTests.jl +++ b/test/mpi/GMGLinearSolversHDivRTTests.jl @@ -101,7 +101,7 @@ function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, println("L2 error = ", e_l2) end - return history.iters, num_free_dofs(Vh) + return history.iters, num_free_dofs(Uh) end end @@ -112,7 +112,7 @@ if !MPI.Initialized() end # Parameters -order = 2 +order = 0 coarse_grid_partition = (2,2) num_refs_coarse = 2 diff --git a/test/seq/DistributedPatchFESpacesDebuggingTests.jl b/test/seq/DistributedPatchFESpacesDebuggingTests.jl new file mode 100644 index 00000000..ecf2202a --- /dev/null +++ b/test/seq/DistributedPatchFESpacesDebuggingTests.jl @@ -0,0 +1,223 @@ +module DistributedPatchFESpacesHDivTests + +using LinearAlgebra +using Test +using PartitionedArrays +using Gridap +using Gridap.Helpers +using Gridap.Arrays +using Gridap.Geometry +using Gridap.ReferenceFEs +using GridapDistributed +using FillArrays + +using GridapSolvers +import GridapSolvers.PatchBasedSmoothers as PBS + +#= +backend = SequentialBackend() +ranks = (1,2) +parts = get_part_ids(backend,ranks) +=# + +function run(parts) + domain = (0.0,1.0,0.0,1.0) + partition = (2,4) + model = CartesianDiscreteModel(parts,domain,partition) + + order = 0 + reffe = ReferenceFE(raviart_thomas,Float64,order) + Vh = TestFESpace(model,reffe) + PD = PBS.PatchDecomposition(model) + Ph = PBS.PatchFESpace(model,reffe,DivConformity(),PD,Vh) + + + # ---- Testing Prolongation and Injection ---- # + + w, w_sums = PBS.compute_weight_operators(Ph,Vh); + + xP = PVector(1.0,Ph.gids) + yP = PVector(0.0,Ph.gids) + x = PVector(1.0,Vh.gids) + y = PVector(0.0,Vh.gids) + + PBS.prolongate!(yP,Ph,x) + PBS.inject!(y,Ph,yP,w,w_sums) + @test x ≈ y + + PBS.inject!(x,Ph,xP,w,w_sums) + PBS.prolongate!(yP,Ph,x) + @test xP ≈ yP + + # ---- Assemble systems ---- # + + sol(x) = VectorValue(x[1],x[2]) + f(x) = VectorValue(2.0*x[2]*(1.0-x[1]*x[1]),2.0*x[1]*(1-x[2]*x[2])) + + α = 1.0 + biform(u,v,dΩ) = ∫(v⋅u)dΩ + ∫(α*divergence(v)⋅divergence(u))dΩ + liform(v,dΩ) = ∫(v⋅f)dΩ + + Ω = Triangulation(model) + dΩ = Measure(Ω,2*order+1) + a(u,v) = biform(u,v,dΩ) + l(v) = liform(v,dΩ) + + assembler = SparseMatrixAssembler(Vh,Vh) + Ah = assemble_matrix(a,assembler,Vh,Vh) + fh = assemble_vector(l,assembler,Vh) + + sol_h = solve(LUSolver(),Ah,fh) + + Ωₚ = Triangulation(PD) + dΩₚ = Measure(Ωₚ,2*order+1) + ap(u,v) = biform(u,v,dΩₚ) + lp(v) = liform(v,dΩₚ) + + assembler_P = SparseMatrixAssembler(Ph,Ph,FullyAssembledRows()) + Ahp = assemble_matrix(ap,assembler_P,Ph,Ph) + fhp = assemble_vector(lp,assembler_P,Ph) + + sol_hp = solve(LUSolver(),Ahp,fhp) + + # ---- Define solvers ---- # + + LU = LUSolver() + LUss = symbolic_setup(LU,Ahp) + LUns = numerical_setup(LUss,Ahp) + + M = PBS.PatchBasedLinearSolver(ap,Ph,Vh,LU) + Mss = symbolic_setup(M,Ah) + Mns = numerical_setup(Mss,Ah) + + R = RichardsonSmoother(M,10,1.0/3.0) + Rss = symbolic_setup(R,Ah) + Rns = numerical_setup(Rss,Ah) + + # ---- Manual solve using LU ---- # + + x1_mat = PVector(0.5,Ah.cols) + r1_mat = fh-Ah*x1_mat + exchange!(r1_mat) + + r1 = PVector(0.0,Vh.gids) + x1 = PVector(0.0,Vh.gids) + rp = PVector(0.0,Ph.gids) + xp = PVector(0.0,Ph.gids) + rp_mat = PVector(0.0,Ahp.cols) + xp_mat = PVector(0.0,Ahp.cols) + + copy!(r1,r1_mat) + exchange!(r1) + PBS.prolongate!(rp,Ph,r1) # OK + + copy!(rp_mat,rp) + exchange!(rp_mat) + solve!(xp_mat,LUns,rp_mat) + copy!(xp,xp_mat) # Some big numbers appear here.... + + w, w_sums = PBS.compute_weight_operators(Ph,Vh); + PBS.inject!(x1,Ph,xp,w,w_sums) # Problem here!! + copy!(x1_mat,x1) + + # ---- Same using the PatchBasedSmoother ---- # + + x2_mat = PVector(0.5,Ah.cols) + r2_mat = fh-Ah*x2_mat + exchange!(r2_mat) + solve!(x2_mat,Mns,r2_mat) + + # ---- Smoother inside Richardson + + x3_mat = PVector(0.5,Ah.cols) + r3_mat = fh-Ah*x3_mat + exchange!(r3_mat) + solve!(x3_mat,Rns,r3_mat) + exchange!(x3_mat) + + # Outputs + res = Dict{String,Any}() + res["sol_h"] = sol_h + res["sol_hp"] = sol_hp + + res["r1"] = r1 + res["x1"] = x1 + res["r1_mat"] = r1_mat + res["x1_mat"] = x1_mat + res["rp"] = rp + res["xp"] = xp + res["rp_mat"] = rp_mat + res["xp_mat"] = xp_mat + + res["w"] = w + res["w_sums"] = w_sums + + return model,PD,Ph,Vh,res +end + +backend = SequentialBackend() + +parts = get_part_ids(backend,(1,1)) +Ms,PDs,Phs,Vhs,res_single = run(parts); + +parts = get_part_ids(backend,(1,2)) +Mm,PDm,Phm,Vhm,res_multi = run(parts); + +println(repeat('#',80)) + +map_parts(local_views(Ms)) do model + cell_ids = get_cell_node_ids(model) + cell_coords = get_cell_coordinates(model) + display(reshape(cell_ids,length(cell_ids))) + display(reshape(cell_coords,length(cell_coords))) +end; +println(repeat('-',80)) + +cell_gids = get_cell_gids(Mm) +vertex_gids = get_face_gids(Mm,0) +edge_gids = get_face_gids(Mm,1) + +println(">>> Cell gids") +map_parts(cell_gids.partition) do p + println(transpose(p.lid_to_ohid)) +end; +println(repeat('-',80)) + +println(">>> Vertex gids") +map_parts(vertex_gids.partition) do p + println(transpose(p.lid_to_ohid)) +end; +println(repeat('-',80)) + +println(">>> Edge gids") +map_parts(edge_gids.partition) do p + println(transpose(p.lid_to_ohid)) +end; + +println(repeat('#',80)) + +map_parts(local_views(Phs)) do Ph + display(Ph.patch_cell_dofs_ids) +end; + +map_parts(local_views(Phm)) do Ph + display(Ph.patch_cell_dofs_ids) +end; + +println(repeat('#',80)) + +for key in keys(res_single) + val_s = res_single[key] + val_m = res_multi[key] + + println(">>> ", key) + map_parts(val_s.values) do v + println(transpose(v)) + end; + map_parts(val_m.owned_values) do v + println(transpose(v)) + end; + println(repeat('-',80)) +end + +end \ No newline at end of file diff --git a/test/seq/DistributedPatchFESpacesTests.jl b/test/seq/DistributedPatchFESpacesTests.jl index dc0a7dcc..18dc37c3 100644 --- a/test/seq/DistributedPatchFESpacesTests.jl +++ b/test/seq/DistributedPatchFESpacesTests.jl @@ -6,6 +6,7 @@ using PartitionedArrays using Gridap using Gridap.Helpers using Gridap.Geometry +using Gridap.ReferenceFEs using GridapDistributed using FillArrays @@ -20,13 +21,16 @@ domain = (0.0,1.0,0.0,1.0) partition = (2,4) model = CartesianDiscreteModel(parts,domain,partition) -order = 1 -reffe = ReferenceFE(lagrangian,Float64,order) +# order = 1 +# reffe = ReferenceFE(lagrangian,Float64,order) +order = 0 +reffe = ReferenceFE(raviart_thomas,Float64,order) Vh = TestFESpace(model,reffe) -PD = PBS.PatchDecomposition(model)#,patch_boundary_style=PBS.PatchBoundaryInclude()) -Ph = PBS.PatchFESpace(model,reffe,Gridap.ReferenceFEs.H1Conformity(),PD,Vh) +PD = PBS.PatchDecomposition(model) +Ph = PBS.PatchFESpace(model,reffe,DivConformity(),PD,Vh) # ---- Testing Prolongation and Injection ---- # + w, w_sums = PBS.compute_weight_operators(Ph,Vh); xP = PVector(1.0,Ph.gids) @@ -42,7 +46,9 @@ PBS.inject!(x,Ph,xP,w,w_sums) PBS.prolongate!(yP,Ph,x) @test xP ≈ yP + # ---- Assemble systems ---- # + Ω = Triangulation(model) dΩ = Measure(Ω,2*order+1) a(u,v) = ∫(v⋅u)*dΩ @@ -63,7 +69,9 @@ assembler_P = SparseMatrixAssembler(Ph,Ph) Ahp = assemble_matrix(ap,assembler_P,Ph,Ph) fhp = assemble_vector(lp,assembler_P,Ph) + # ---- Define solvers ---- # + LU = LUSolver() LUss = symbolic_setup(LU,Ahp) LUns = numerical_setup(LUss,Ahp) @@ -76,6 +84,7 @@ R = RichardsonSmoother(M,10,1.0/3.0) Rss = symbolic_setup(R,Ah) Rns = numerical_setup(Rss,Ah) + # ---- Manual solve using LU ---- # x1_mat = PVector(0.5,Ah.cols) @@ -101,6 +110,7 @@ w, w_sums = PBS.compute_weight_operators(Ph,Vh); PBS.inject!(x1,Ph,xp,w,w_sums) copy!(x1_mat,x1) + # ---- Same using the PatchBasedSmoother ---- # x2_mat = PVector(0.5,Ah.cols) @@ -108,7 +118,9 @@ r2_mat = fh-Ah*x2_mat exchange!(r2_mat) solve!(x2_mat,Mns,r2_mat) + # ---- Smoother inside Richardson + x3_mat = PVector(0.5,Ah.cols) r3_mat = fh-Ah*x3_mat exchange!(r3_mat) From fffaeaf64c71f9b5ec687a3dca994b160a0ad558 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 3 Mar 2023 14:57:36 +1100 Subject: [PATCH 12/14] Small bugfix --- src/MultilevelTools/GridapDistributedExtensions.jl | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/MultilevelTools/GridapDistributedExtensions.jl b/src/MultilevelTools/GridapDistributedExtensions.jl index b7e163a0..56fa3264 100644 --- a/src/MultilevelTools/GridapDistributedExtensions.jl +++ b/src/MultilevelTools/GridapDistributedExtensions.jl @@ -12,6 +12,10 @@ function Gridap.CellData.Measure(tt::GridapDistributed.DistributedTriangulation{ return GridapDistributed.DistributedMeasure(measures) end +function GridapDistributed.get_parts(x::GridapDistributed.DistributedFESpace) + PartitionedArrays.get_part_ids(local_views(x)) +end + # change_parts function change_parts(x::Union{AbstractPData,Nothing}, new_parts; default=nothing) From c0134382a99d3d89f3a72b2cd02020e372c0fc26 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 3 Mar 2023 16:57:19 +1100 Subject: [PATCH 13/14] Small bugfix in tests --- test/mpi/RefinementToolsTests.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/mpi/RefinementToolsTests.jl b/test/mpi/RefinementToolsTests.jl index 235d9e6e..06a02165 100644 --- a/test/mpi/RefinementToolsTests.jl +++ b/test/mpi/RefinementToolsTests.jl @@ -32,15 +32,17 @@ function run(parts,num_parts_x_level,coarse_grid_partition,num_refs_coarse) cparts = get_level_parts(mh,lev+1) if i_am_in(cparts) + model_h = get_model_before_redist(mh,lev) Vh = get_fe_space_before_redist(tests,lev) Uh = get_fe_space_before_redist(trials,lev) - Ωh = get_triangulation(Uh,get_model_before_redist(mh,lev)) + Ωh = get_triangulation(model_h) dΩh = Measure(Ωh,quad_order) uh = interpolate(sol,Uh) + model_H = get_model(mh,lev+1) VH = get_fe_space(tests,lev+1) UH = get_fe_space(trials,lev+1) - ΩH = get_triangulation(UH,get_model(mh,lev+1)) + ΩH = get_triangulation(model_H) dΩH = Measure(ΩH,quad_order) uH = interpolate(sol,UH) dΩhH = Measure(ΩH,Ωh,quad_order) From 5d7ce2101bc138cac0fd03f9da471dcead8b2c81 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Fri, 3 Mar 2023 22:29:39 +1030 Subject: [PATCH 14/14] Minor bugfix --- src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl | 2 +- test/seq/PatchLinearSolverTests.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl index 9178178d..867f0f5d 100644 --- a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl +++ b/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl @@ -100,7 +100,7 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumerical rp, dxp = caches prolongate!(rp,Ph,r) - solve!(dxp,nsAp,rp) + solve!(dxp,Ap_ns,rp) inject!(x,Ph,dxp,w,w_sums) return x diff --git a/test/seq/PatchLinearSolverTests.jl b/test/seq/PatchLinearSolverTests.jl index 015e5134..5df48824 100644 --- a/test/seq/PatchLinearSolverTests.jl +++ b/test/seq/PatchLinearSolverTests.jl @@ -82,7 +82,7 @@ module PatchLinearSolverTests A,b = compute_matrix_vector(model,Vh) x = test_smoother(PD,Ph,Vh,A,b) - parts = get_part_ids(sequential,(1,1)) + parts = get_part_ids(SequentialBackend(),(1,1)) dmodel = CartesianDiscreteModel(parts,domain,partition) dPD,dPh,dxh,dVh = returns_PD_Ph_xh_Vh(dmodel); dA,db = compute_matrix_vector(dmodel,dVh);