Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Optimizing patch-based smoothers #10

Merged
merged 4 commits into from
Mar 11, 2023
Merged

Optimizing patch-based smoothers #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
f5f1808
Optimized injection for patch-based fespaces
JordiManyer Mar 10, 2023
2b18131
Removed unnecessary caches in the solver
JordiManyer Mar 10, 2023
442f008
Further optimized patch-based injection
JordiManyer Mar 10, 2023
86a0af7
Optimized weights for distributed
JordiManyer Mar 10, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 10 additions & 23 deletions src/PatchBasedSmoothers/mpi/PatchFESpaces.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,32 +14,21 @@ function PatchFESpace(model::GridapDistributed.DistributedDiscreteModel,
Vh::GridapDistributed.DistributedSingleFieldFESpace)
root_gids = get_face_gids(model,get_patch_root_dim(patch_decomposition))

function f(model,patch_decomposition,Vh,partition)
spaces = map_parts(local_views(model),
local_views(patch_decomposition),
local_views(Vh),
root_gids.partition) do model, patch_decomposition, Vh, partition
patches_mask = fill(false,length(partition.lid_to_gid))
patches_mask[partition.hid_to_lid] .= true # Mask ghost patch roots
PatchFESpace(model,
reffe,
conformity,
patch_decomposition,
Vh;
patches_mask=patches_mask)
PatchFESpace(model,reffe,conformity,patch_decomposition,Vh;patches_mask=patches_mask)
end

spaces = map_parts(f,
local_views(model),
local_views(patch_decomposition),
local_views(Vh),
root_gids.partition)

parts = get_parts(model)
nodofs = map_parts(spaces) do space
num_free_dofs(space)
end
ngdofs = sum(nodofs)

first_gdof, _ = xscan(+,reduce,nodofs,init=1)
local_ndofs = map_parts(num_free_dofs,spaces)
global_ndofs = sum(local_ndofs)
first_gdof, _ = xscan(+,reduce,local_ndofs,init=1)
# This PRange has no ghost dofs
gids = PRange(parts,ngdofs,nodofs,first_gdof)
gids = PRange(parts,global_ndofs,local_ndofs,first_gdof)
return GridapDistributed.DistributedSingleFieldFESpace(spaces,gids,get_vector_type(Vh))
end

Expand Down Expand Up @@ -98,9 +87,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,Ph.Vh)
w .= _w
w_sums .= _w_sums
compute_weight_operators!(Ph,Ph.Vh,w,w_sums)
end

# partial sums -> global sums
Expand Down
10 changes: 3 additions & 7 deletions src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,13 +63,11 @@ end
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)
r = PVector(0.0,Vh.gids)
x = PVector(0.0,Vh.gids)
return rp_mat, dxp_mat, rp, dxp, r, x
return rp, dxp, r, x
end

function _allocate_col_vector(A::AbstractMatrix)
Expand Down Expand Up @@ -111,14 +109,12 @@ function Gridap.Algebra.solve!(x_mat::PVector,ns::PatchBasedSmootherNumericalSet

Ph = ns.solver.Ph
w, w_sums = weights
rp_mat, dxp_mat, rp, dxp, r, x = caches
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)
solve!(dxp,Ap_ns,rp)
inject!(x,Ph,dxp,w,w_sums)
copy!(x_mat,x)

Expand Down
155 changes: 95 additions & 60 deletions src/PatchBasedSmoothers/seq/PatchFESpaces.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@ struct PatchFESpace <: Gridap.FESpaces.SingleFieldFESpace
patch_cell_dofs_ids :: Gridap.Arrays.Table
Vh :: Gridap.FESpaces.SingleFieldFESpace
patch_decomposition :: PatchDecomposition
dof_to_pdof :: Gridap.Arrays.Table
end

# INPUT
Expand Down Expand Up @@ -68,7 +69,10 @@ function PatchFESpace(model::DiscreteModel,
cell_conformity,
patches_mask)

return PatchFESpace(num_dofs,patch_cell_dofs_ids,Vh,patch_decomposition)
dof_to_pdof = allocate_dof_to_pdof(Vh,patch_decomposition,patch_cell_dofs_ids)
generate_dof_to_pdof!(dof_to_pdof,Vh,patch_decomposition,patch_cell_dofs_ids)

return PatchFESpace(num_dofs,patch_cell_dofs_ids,Vh,patch_decomposition,dof_to_pdof)
end

Gridap.FESpaces.get_dof_value_type(a::PatchFESpace) = Gridap.FESpaces.get_dof_value_type(a.Vh)
Expand Down Expand Up @@ -106,7 +110,7 @@ function allocate_patch_cell_dofs_ids(num_cells_overlapped_mesh,cell_patches,cel
end
end

Gridap.Helpers.@check num_cells_overlapped_mesh+1 == gcell_overlapped_mesh
@check num_cells_overlapped_mesh+1 == gcell_overlapped_mesh
data = Vector{Int}(undef,ptrs[end]-1)
return Gridap.Arrays.Table(data,ptrs)
end
Expand Down Expand Up @@ -159,9 +163,6 @@ function generate_patch_cell_dofs_ids!(patch_cell_dofs_ids,
free_dofs_offset=1,
mask=false)

patch_global_space_cell_dofs_ids=
lazy_map(Broadcasting(Reindex(global_space_cell_dofs_ids)),patch_cells)

o = patch_cells_overlapped_mesh.ptrs[patch]
if mask
for lpatch_cell = 1:length(patch_cells)
Expand All @@ -186,15 +187,17 @@ function generate_patch_cell_dofs_ids!(patch_cell_dofs_ids,
cells_d_faces = Gridap.Geometry.get_faces(topology,Dc,d)
cell_d_face = cells_d_faces[patch_cell]

# 1) DoFs belonging to faces (Df < Dc)
for (lf,f) in enumerate(cell_d_face)
# If current face is on the patch boundary
# A) If current face is on the patch boundary
if (patch_cells_faces_on_boundary[d+1][cell_overlapped_mesh][lf])
# assign negative indices to DoFs owned by face
for ldof in cell_conformity.ctype_lface_own_ldofs[ctype][face_offset+lf]
gdof = global_space_cell_dofs_ids[patch_cell][ldof]
current_patch_cell_dofs_ids[ldof] = -1
end
else
# B) If current face is not in patch boundary,
# rely on the existing glued info (available at global_space_cell_dof_ids)
# (we will need a Dict{Int,Int} to hold the correspondence among global
# space and patch cell dofs IDs)
Expand All @@ -217,7 +220,7 @@ function generate_patch_cell_dofs_ids!(patch_cell_dofs_ids,
face_offset += cell_conformity.d_ctype_num_dfaces[d+1][ctype]
end

# Interior DoFs
# 2) Interior DoFs
for ldof in cell_conformity.ctype_lface_own_ldofs[ctype][face_offset+1]
current_patch_cell_dofs_ids[ldof] = free_dofs_offset
free_dofs_offset += 1
Expand All @@ -227,6 +230,68 @@ function generate_patch_cell_dofs_ids!(patch_cell_dofs_ids,
return free_dofs_offset
end

function allocate_dof_to_pdof(Vh,PD,patch_cell_dofs_ids)
touched = Dict{Int,Bool}()
cell_mesh_overlapped = 1
cache_patch_cells = array_cache(PD.patch_cells)
cell_dof_ids = get_cell_dof_ids(Vh)
cache_cell_dof_ids = array_cache(cell_dof_ids)

ptrs = fill(0,num_free_dofs(Vh)+1)
for patch = 1:length(PD.patch_cells)
current_patch_cells = getindex!(cache_patch_cells,PD.patch_cells,patch)
for cell in current_patch_cells
current_cell_dof_ids = getindex!(cache_cell_dof_ids,cell_dof_ids,cell)
s = patch_cell_dofs_ids.ptrs[cell_mesh_overlapped]
e = patch_cell_dofs_ids.ptrs[cell_mesh_overlapped+1]-1
current_patch_cell_dof_ids = view(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 ∈ keys(touched))
touched[dof] = true
ptrs[dof+1] += 1
end
end
cell_mesh_overlapped += 1
end
empty!(touched)
end
PartitionedArrays.length_to_ptrs!(ptrs)

data = fill(0,ptrs[end]-1)
return Gridap.Arrays.Table(data,ptrs)
end

function generate_dof_to_pdof!(dof_to_pdof,Vh,PD,patch_cell_dofs_ids)
touched = Dict{Int,Bool}()
cell_mesh_overlapped = 1
cache_patch_cells = array_cache(PD.patch_cells)
cell_dof_ids = get_cell_dof_ids(Vh)
cache_cell_dof_ids = array_cache(cell_dof_ids)

ptrs = dof_to_pdof.ptrs
data = dof_to_pdof.data
local_ptrs = fill(Int32(0),num_free_dofs(Vh))
for patch = 1:length(PD.patch_cells)
current_patch_cells = getindex!(cache_patch_cells,PD.patch_cells,patch)
for cell in current_patch_cells
current_cell_dof_ids = getindex!(cache_cell_dof_ids,cell_dof_ids,cell)
s = patch_cell_dofs_ids.ptrs[cell_mesh_overlapped]
e = patch_cell_dofs_ids.ptrs[cell_mesh_overlapped+1]-1
current_patch_cell_dof_ids = view(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 ∈ keys(touched))
touched[dof] = true
idx = ptrs[dof] + local_ptrs[dof]
@check idx < ptrs[dof+1]
data[idx] = pdof
local_ptrs[dof] += 1
end
end
cell_mesh_overlapped += 1
end
empty!(touched)
end
end

# x \in PatchFESpace
# y \in SingleFESpace
Expand Down Expand Up @@ -261,71 +326,41 @@ function inject!(x,Ph::PatchFESpace,y)
end

function inject!(x,Ph::PatchFESpace,y,w)
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,
patch)
for cell in current_patch_cells
current_cell_dof_ids = getindex!(cache_cell_dof_ids,cell_dof_ids,cell)
s = Ph.patch_cell_dofs_ids.ptrs[cell_mesh_overlapped]
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 ∈ keys(touched))
touched[dof] = true
x[dof] += y[pdof] * w[pdof]
end
end
cell_mesh_overlapped += 1
dof_to_pdof = Ph.dof_to_pdof

ptrs = dof_to_pdof.ptrs
data = dof_to_pdof.data
for dof in 1:length(dof_to_pdof)
x[dof] = 0.0
for k in ptrs[dof]:ptrs[dof+1]-1
pdof = data[k]
x[dof] += y[pdof] * w[pdof]
end
empty!(touched)
end
end

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)
dof_to_pdof = Ph.dof_to_pdof

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,
patch)
for cell in current_patch_cells
current_cell_dof_ids = getindex!(cache_cell_dof_ids,cell_dof_ids,cell)
s = Ph.patch_cell_dofs_ids.ptrs[cell_mesh_overlapped]
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 ∈ keys(touched))
touched[dof] = true
x[dof] += y[pdof] * w[pdof] / w_sums[dof]
end
end
cell_mesh_overlapped += 1
ptrs = dof_to_pdof.ptrs
data = dof_to_pdof.data
for dof in 1:length(dof_to_pdof)
x[dof] = 0.0
for k in ptrs[dof]:ptrs[dof+1]-1
pdof = data[k]
x[dof] += y[pdof] * w[pdof] / w_sums[dof]
end
empty!(touched)
end
end

function compute_weight_operators(Ph::PatchFESpace,Vh)
w = Fill(1.0,num_free_dofs(Ph))
w_sums = compute_partial_sums(Ph,Vh,w)
w_sums = zeros(num_free_dofs(Vh))
inject!(w_sums,Ph,w,Fill(1.0,num_free_dofs(Vh)))
return w, w_sums
end

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
function compute_weight_operators!(Ph::PatchFESpace,Vh,w,w_sums)
fill!(w,1.0)
inject!(w_sums,Ph,w,Fill(1.0,num_free_dofs(Ph)))
end