From ae2baec7125a9244e15c39143f36626a0273fb53 Mon Sep 17 00:00:00 2001 From: Francesc Verdugo Date: Tue, 3 Oct 2023 15:31:53 +0200 Subject: [PATCH] Starting to add poisson with 2 level bddc --- examples/poisson_bddc.jl | 162 ++++++++++++++++++++++++++++++++++++++ src/GalerkinToolkit.jl | 71 ----------------- test/poisson_bddc_test.jl | 31 ++++++++ 3 files changed, 193 insertions(+), 71 deletions(-) create mode 100644 examples/poisson_bddc.jl create mode 100644 test/poisson_bddc_test.jl diff --git a/examples/poisson_bddc.jl b/examples/poisson_bddc.jl new file mode 100644 index 00000000..0b6cf859 --- /dev/null +++ b/examples/poisson_bddc.jl @@ -0,0 +1,162 @@ +module PoissonBDDC + +import GalerkinToolkit as glk +import ForwardDiff +using StaticArrays +using LinearAlgebra +using SparseArrays +using WriteVTK +using PartitionedArrays +using Metis + +const CORNER = 0 +const EDGE = 1 +const FACE = 2 + +function main(params_in) + + # Process params + params_default = default_params() + params = add_default_params(params_in,params_default) + + # Setup + pmesh = params[:pmesh] + bddc_type = params[:bddc_type] + u = params[:u] + ranks = linear_indices(pmesh) + np = length(ranks) + + # Setup dirichlet_bcs + dirichlet_tags = params[:dirichlet_tags] + dirichlet_bcs = map(pmesh) do mesh + node_to_tag = zeros(glk.num_nodes(mesh)) + tag_to_name = dirichlet_tags + glk.classify_mesh_nodes!(node_to_tag,mesh,tag_to_name) + free_and_dirichlet_nodes = glk.partition_from_mask(i->i==0,node_to_tag) + node_to_x = glk.node_coordinates(mesh) + dirichlet_nodes = last(free_and_dirichlet_nodes) + x_dirichlet = view(node_to_x,dirichlet_nodes) + u_dirichlet = u.(x_dirichlet) + (;free_and_dirichlet_nodes,u_dirichlet,node_to_tag) + end + + coarse_dofs = setup_coarse_dofs(pmesh,dirichlet_bcs,bddc_type) + + display(coarse_dofs) + + # Visualize mesh + example_path = params[:example_path] + map(pmesh,ranks) do mesh,rank + pvtk_grid(example_path,glk.vtk_args(mesh)...;part=rank,nparts=np) do vtk + glk.vtk_physical_groups!(vtk,mesh) + vtk["piece",VTKCellData()] = fill(rank,sum(glk.num_faces(mesh))) + vtk["owner",VTKPointData()] = local_to_owner(glk.local_nodes(mesh)) + vtk["interface",VTKPointData()] = map(colors->Int(length(colors)!=1),glk.local_node_colors(mesh)) + end + end + +end + +function add_default_params(params_in,params_default) + UmD = setdiff(keys(params_in),keys(params_default)) + for k in UmD + @warn "Parameter with key :$k is unused" + end + UiD = intersect(keys(params_default),keys(params_in)) + DmU = setdiff(keys(params_default),keys(params_in)) + a = [ k=>params_in[k] for k in UiD] + b = [ k=>params_default[k] for k in DmU] + Dict(vcat(a,b)) +end + +function default_params() + mesh = glk.cartesian_mesh((0,1,0,1),(10,10),complexify=false) + np = 4 + ranks = DebugArray(LinearIndices((np,))) + pmesh = glk.partition_mesh(Metis.partition,ranks,mesh,via=:cells) + outdir = mkpath(joinpath(@__DIR__,"..","output")) + params = Dict{Symbol,Any}() + params[:pmesh] = glk.cartesian_mesh((0,1,0,1),(10,10)) + params[:u] = (x) -> sum(x) + params[:f] = (x) -> 0.0 + params[:g] = (x) -> 0.0 + params[:dirichlet_tags] = ["boundary"] + params[:neumann_tags] = String[] + params[:integration_degree] = 2 + params[:solve] = \ + params[:example_path] = joinpath(outdir,"poisson_parallel") + params[:bddc_type] = (CORNER,EDGE,FACE) + params +end + +function setup_coarse_dofs(pmesh,dirichlet_bcs,types) + ranks = linear_indices(pmesh) + lnode_to_colors = map(glk.local_node_colors,pmesh) + node_partition = map(glk.local_nodes,pmesh) + cldof_to_ldofs, cldof_to_owner, nown = map(ranks,lnode_to_colors,dirichlet_bcs) do rank,lnode_to_colors,dirichlet_bcs + free_and_dirichlet_lnodes = dirichlet_bcs.free_and_dirichlet_nodes + nldofs = length(first(free_and_dirichlet_lnodes)) + lnode_to_ldof = glk.permutation(free_and_dirichlet_lnodes) + nlnodes = length(lnode_to_colors) + interface_node_to_lnode = findall(1:nlnodes) do lnode + colors = lnode_to_colors[lnode] + ldof = lnode_to_ldof[lnode] + length(colors)!=1 && ldof <= nldofs + end + interface_node_to_colors = view(lnode_to_colors,interface_node_to_lnode) + # assumes sorted colors + clnode_to_colors = unique(interface_node_to_colors) + clnode_to_owner = map(maximum,clnode_to_colors) + nclnodes = length(clnode_to_colors) + interface_node_to_clnode = indexin(interface_node_to_colors,clnode_to_colors) + clnode_to_interface_nodes = glk.inverse_index_map(interface_node_to_clnode,nclnodes) + f = interface_node->lnode_to_ldof[interface_node_to_lnode[interface_node]] + clnode_to_interface_nodes.data .= f.(clnode_to_interface_nodes.data) + clnode_to_ldofs = clnode_to_interface_nodes + clnode_to_type = fill(EDGE,nclnodes) + clnode_to_type[ map(i->length(i) == 2,clnode_to_colors) ] .= FACE + clnode_to_type[ map(i->length(i) == 1,clnode_to_ldofs) ] .= CORNER + mask = map(i->i in types,clnode_to_type) + my_cldof_to_ldofs = JaggedArray(clnode_to_ldofs[mask]) + my_cldof_to_owner = clnode_to_owner[mask] + my_nown = count(owner->owner==rank,my_cldof_to_owner) + my_cldof_to_ldofs, my_cldof_to_owner, my_nown + end |> tuple_of_arrays + ncdofs = sum(nown) + cdof_partition = variable_partition(nown,ncdofs) + v = PVector{Vector{Int}}(undef,node_partition) + map(ranks,cdof_partition,cldof_to_owner,cldof_to_ldofs,local_values(v),dirichlet_bcs) do rank, codofs, cldof_to_owner, cldof_to_ldofs, lnode_to_v, dirichlet_bcs + ldof_to_lnode = first(dirichlet_bcs.free_and_dirichlet_nodes) + codof_to_cdof = own_to_global(codofs) + codof = 0 + for (cldof,owner) in enumerate(cldof_to_owner) + if owner != rank + continue + end + codof += 1 + cdof = codof_to_cdof[codof] + ldofs = cldof_to_ldofs[cldof] + lnodes = view(ldof_to_lnode,ldofs) + lnode_to_v[lnodes] .= cdof + end + end + consistent!(v) |> wait + cldof_to_cdof = map(ranks,cldof_to_ldofs,local_values(v),dirichlet_bcs) do rank, cldof_to_ldofs, lnode_to_v,dirichlet_bcs + ldof_to_lnode = first(dirichlet_bcs.free_and_dirichlet_nodes) + ncldofs = length(cldof_to_ldofs) + my_cldof_to_cdof = zeros(Int,ncldofs) + for cldof in 1:ncldofs + ldofs = cldof_to_ldofs[cldof] + cdof = lnode_to_v[ldof_to_lnode[first(ldofs)]] + my_cldof_to_cdof[cldof] = cdof + end + my_cldof_to_cdof + end + rank_to_cdofs = gather(cldof_to_cdof,destination=MAIN) + map(cldof_to_ldofs,rank_to_cdofs) do cldof_to_ldofs,rank_to_cdofs + (;cldof_to_ldofs, rank_to_cdofs, num_cdofs = ncdofs) + end +end + + +end # module diff --git a/src/GalerkinToolkit.jl b/src/GalerkinToolkit.jl index ea4dcad2..7ea961d0 100644 --- a/src/GalerkinToolkit.jl +++ b/src/GalerkinToolkit.jl @@ -3241,77 +3241,6 @@ function partition_mesh_via_cells(colorize,ranks,mesh) pmesh end -function coarse_objects_from_pmesh(pmesh) - ranks = linear_indices(pmesh) - lnode_to_colors = map(local_node_colors,pmesh) - node_partition = map(local_nodes,pmesh) - clnode_to_lnodes, clnode_to_owner, nown, clnode_to_type = map(ranks,lnode_to_colors) do rank,lnode_to_colors - interface_node_to_lnode = findall(colors->length(colors)!=1,lnode_to_colors) - interface_node_to_colors = view(lnode_to_colors,interface_node_to_lnode) - # assumes sorted colors - clnode_to_colors = unique(interface_node_to_colors) - nclnodes = length(clnode_to_colors) - interface_node_to_clnode = indexin(interface_node_to_colors,clnode_to_colors) - clnode_to_interface_nodes = inverse_index_map(interface_node_to_clnode,nclnodes) - f = interface_node->interface_node_to_lnode[interface_node] - clnode_to_interface_nodes.data .= f.(clnode_to_interface_nodes.data) - my_clnode_to_lnodes = clnode_to_interface_nodes - my_clnode_to_owner = map(maximum,clnode_to_colors) - CORNER = 0 - EDGE = 1 - FACE = 2 - my_clnode_to_type = fill(EDGE,nclnodes) - my_clnode_to_type[ map(i->length(i) == 2,clnode_to_colors) ] .= FACE - my_clnode_to_type[ map(i->length(i) == 1,my_clnode_to_lnodes) ] .= CORNER - my_nown = count(owner->owner==rank,my_clnode_to_owner) - my_clnode_to_lnodes, my_clnode_to_owner, my_nown, my_clnode_to_type - end |> tuple_of_arrays - ncnodes = sum(nown) - cnode_partition = variable_partition(nown,ncnodes) - v = PVector{Vector{Int}}(undef,node_partition) - map(ranks,cnode_partition,clnode_to_owner,clnode_to_lnodes,local_values(v)) do rank, conodes, clnode_to_owner, clnode_to_lnodes, lnode_to_v - conode_to_cnode = own_to_global(conodes) - conode = 0 - for (clnode,owner) in enumerate(clnode_to_owner) - if owner != rank - continue - end - conode += 1 - cnode = conode_to_cnode[conode] - lnodes = clnode_to_lnodes[clnode] - lnode_to_v[lnodes] .= cnode - end - end - consistent!(v) |> wait - clnode_to_cnode = map(ranks,clnode_to_lnodes,local_values(v)) do rank, clnode_to_lnodes, lnode_to_v - nclnodes = length(clnode_to_lnodes) - my_clnode_to_cnode = zeros(Int,nclnodes) - for clnode in 1:nclnodes - lnodes = clnode_to_lnodes[clnode] - cnode = lnode_to_v[first(lnodes)] - my_clnode_to_cnode[clnode] = cnode - end - my_clnode_to_cnode - end - rank_to_cnodes = gather(clnode_to_cnode,destination=MAIN) - rank_to_types = gather(clnode_to_type,destination=MAIN) - cnode_to_type = map(rank_to_types,rank_to_cnodes,ranks) do rank_clnode_to_type,rank_clnode_to_cnode,rank - if rank == MAIN - my_cnode_to_type = fill(Int(INVALID_ID),ncnodes) - for part in length(rank_clnode_to_type) - my_cnode_to_type[rank_clnode_to_cnode[part]] = rank_clnode_to_type[part] - end - my_cnode_to_type - else - my_cnode_to_type = Int[] - end - end - - - map(clnode_to_lnodes,rank_to_cnodes,clnode_to_type,cnode_to_type) do clnode_to_lnodes,rank_to_cnodes,clnode_to_type,cnode_to_type - (;clnode_to_lnodes, rank_to_cnodes, num_cnodes = ncnodes,clnode_to_type,cnode_to_type) - end -end function inverse_index_map(a_to_b,nb) b_to_as_ptrs = zeros(Int32,nb+1) diff --git a/test/poisson_bddc_test.jl b/test/poisson_bddc_test.jl new file mode 100644 index 00000000..fad0f2ee --- /dev/null +++ b/test/poisson_bddc_test.jl @@ -0,0 +1,31 @@ +module PoissonParallelTests + +using Test +import GalerkinToolkit as gt +using PartitionedArrays +using Metis + +include("../examples/poisson_bddc.jl") + +using .PoissonBDDC: CORNER,FACE + +with_debug() do distribute + + params = Dict{Symbol,Any}() + mesh = gt.cartesian_mesh((0,2,0,2),(20,20),complexify=false) + np = 4 + ranks = distribute(LinearIndices((np,))) + pmesh = gt.partition_mesh(Metis.partition,ranks,mesh,via=:cells) + params[:pmesh] = pmesh + params[:dirichlet_tags] = ["1-face-1","1-face-3","1-face-4"] + params[:neumann_tags] = ["1-face-2"] + params[:u] = (x) -> sum(x) + params[:f] = (x) -> 0.0 + params[:g] = (x) -> 1.0 + params[:bddc_type] = (FACE,) + @time results = PoissonBDDC.main(params) + +end + + +end # module