diff --git a/NEWS.md b/NEWS.md index 08bf90b..d52db39 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). +## [0.4.2] 2024-07-4 + +### Added + +- Added uniform anisotropic refinement of distributed cartesian meshes. Since PR[#148](https://github.com/gridap/GridapDistributed.jl/pull/148). + ## [0.4.1] 2024-06-25 ### Fixed diff --git a/Project.toml b/Project.toml index 4353b29..9255dfd 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridapDistributed" uuid = "f9701e48-63b3-45aa-9a63-9bc6c271f355" authors = ["S. Badia ", "A. F. Martin ", "F. Verdugo "] -version = "0.4.1" +version = "0.4.2" [deps] BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e" diff --git a/src/Adaptivity.jl b/src/Adaptivity.jl index d648260..6e1fe74 100644 --- a/src/Adaptivity.jl +++ b/src/Adaptivity.jl @@ -3,13 +3,32 @@ const DistributedAdaptedDiscreteModel{Dc,Dp} = GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:AdaptedDiscreteModel{Dc,Dp}}} -function DistributedAdaptedDiscreteModel(model :: DistributedDiscreteModel, - parent :: DistributedDiscreteModel, - glue :: AbstractArray{<:AdaptivityGlue}) +function DistributedAdaptedDiscreteModel( + model :: DistributedDiscreteModel, + parent :: DistributedDiscreteModel, + glue :: AbstractArray{<:AdaptivityGlue}; +) models = map(local_views(model),local_views(parent),glue) do model, parent, glue AdaptedDiscreteModel(model,parent,glue) end - return GenericDistributedDiscreteModel(models,get_cell_gids(model)) + gids = get_cell_gids(model) + metadata = model.metadata + return GenericDistributedDiscreteModel(models,gids;metadata) +end + +function Adaptivity.get_model(model::DistributedAdaptedDiscreteModel) + GenericDistributedDiscreteModel( + map(get_model,local_views(model)), + get_cell_gids(model); + metadata=model.metadata + ) +end + +function Adaptivity.get_parent(model::DistributedAdaptedDiscreteModel) + msg = " Error: Cannot get global parent model. \n + We do not keep the global ids of the parent model within the children.\n + You can extract the local parents with map(get_parent,local_views(model))" + @notimplemented msg end function Adaptivity.get_adaptivity_glue(model::DistributedAdaptedDiscreteModel) @@ -45,6 +64,15 @@ end get_parts(g::RedistributeGlue) = get_parts(g.parts_rcv) +function Base.reverse(g::RedistributeGlue) + RedistributeGlue( + g.old_parts,g.new_parts, + g.parts_snd,g.parts_rcv, + g.lids_snd,g.lids_rcv, + g.new2old,g.old2new + ) +end + function get_old_and_new_parts(g::RedistributeGlue,reverse::Val{false}) return g.old_parts, g.new_parts end @@ -82,6 +110,17 @@ function redistribute(::DistributedDiscreteModel,args...;kwargs...) @abstractmethod end +function redistribute(model::DistributedAdaptedDiscreteModel,args...;kwargs...) + # Local cmodels are AdaptedDiscreteModels. To correctly dispatch, we need to + # extract the underlying models, then redistribute. + _model = GenericDistributedDiscreteModel( + map(get_model,local_views(model)), + get_cell_gids(model); + metadata=model.metadata + ) + return redistribute(_model,args...;kwargs...) +end + # Redistribution of cell-wise dofs, free values and FEFunctions function _allocate_cell_wise_dofs(cell_to_ldofs) @@ -177,12 +216,13 @@ function _num_dofs_x_cell(cell_dofs_array,lids) end end -function get_redistribute_cell_dofs_cache(cell_dof_values_old, - cell_dof_ids_new, - model_new, - glue::RedistributeGlue; - reverse=false) - +function get_redistribute_cell_dofs_cache( + cell_dof_values_old, + cell_dof_ids_new, + model_new, + glue::RedistributeGlue; + reverse=false +) lids_rcv, lids_snd, parts_rcv, parts_snd, new2old = get_glue_components(glue,Val(reverse)) cell_dof_values_old = change_parts(cell_dof_values_old,get_parts(glue);default=[]) @@ -199,25 +239,28 @@ function get_redistribute_cell_dofs_cache(cell_dof_values_old, return caches end -function redistribute_cell_dofs(cell_dof_values_old, - cell_dof_ids_new, - model_new, - glue::RedistributeGlue; - reverse=false) +function redistribute_cell_dofs( + cell_dof_values_old, + cell_dof_ids_new, + model_new, + glue::RedistributeGlue; + reverse=false +) caches = get_redistribute_cell_dofs_cache(cell_dof_values_old,cell_dof_ids_new,model_new,glue;reverse=reverse) return redistribute_cell_dofs!(caches,cell_dof_values_old,cell_dof_ids_new,model_new,glue;reverse=reverse) end -function redistribute_cell_dofs!(caches, - cell_dof_values_old, - cell_dof_ids_new, - model_new, - glue::RedistributeGlue; - reverse=false) +function redistribute_cell_dofs!( + caches, + cell_dof_values_old, + cell_dof_ids_new, + model_new, + glue::RedistributeGlue; + reverse=false +) snd_data, rcv_data, cell_dof_values_new = caches lids_rcv, lids_snd, parts_rcv, parts_snd, new2old = get_glue_components(glue,Val(reverse)) - old_parts, new_parts = get_old_and_new_parts(glue,Val(reverse)) cell_dof_values_old = change_parts(cell_dof_values_old,get_parts(glue);default=[]) cell_dof_ids_new = change_parts(cell_dof_ids_new,get_parts(glue);default=[[]]) @@ -231,15 +274,15 @@ function redistribute_cell_dofs!(caches, # We have to build the owned part of "cell_dof_values_new" out of # 1. cell_dof_values_old (for those cells s.t. new2old[:]!=0) # 2. cell_dof_values_new_rcv (for those cells s.t. new2old[:]=0) - _update_cell_dof_values_with_local_info!(cell_dof_values_new, - cell_dof_values_old, - new2old) - + _update_cell_dof_values_with_local_info!( + cell_dof_values_new, cell_dof_values_old, new2old + ) _unpack_rcv_data!(cell_dof_values_new,rcv_data,lids_rcv) # Now that every part knows it's new owned dofs, exchange ghosts - cell_dof_values_new = change_parts(cell_dof_values_new,new_parts) - if i_am_in(new_parts) + if !isnothing(model_new) + new_parts = get_parts(model_new) + cell_dof_values_new = change_parts(cell_dof_values_new,new_parts) cache = fetch_vector_ghost_values_cache(cell_dof_values_new,partition(get_cell_gids(model_new))) fetch_vector_ghost_values!(cell_dof_values_new,cache) |> wait end @@ -260,11 +303,13 @@ function _get_cell_dof_ids_inner_space(s::TrialFESpace) _get_cell_dof_ids_inner_space(s.space) end -function redistribute_fe_function(uh_old::Union{DistributedSingleFieldFEFunction,Nothing}, - Uh_new::Union{DistributedSingleFieldFESpace,Nothing}, - model_new, - glue::RedistributeGlue; - reverse=false) +function redistribute_fe_function( + uh_old::Union{DistributedSingleFieldFEFunction,Nothing}, + Uh_new::Union{DistributedSingleFieldFESpace,Nothing}, + model_new, + glue::RedistributeGlue; + reverse=false +) old_parts, new_parts = get_old_and_new_parts(glue,Val(reverse)) cell_dof_values_old = i_am_in(old_parts) ? map(get_cell_dof_values,local_views(uh_old)) : nothing @@ -290,31 +335,37 @@ for T in [:DistributedSingleFieldFESpace,:DistributedMultiFieldFESpace] end end -function redistribute_free_values(fv_new::Union{PVector,Nothing}, - Uh_new::Union{DistributedSingleFieldFESpace,DistributedMultiFieldFESpace,Nothing}, - fv_old::Union{PVector,Nothing}, - dv_old::Union{AbstractArray,Nothing}, - Uh_old::Union{DistributedSingleFieldFESpace,DistributedMultiFieldFESpace,Nothing}, - model_new, - glue::RedistributeGlue; - reverse=false) +function redistribute_free_values( + fv_new::Union{PVector,Nothing}, + Uh_new::Union{DistributedSingleFieldFESpace,DistributedMultiFieldFESpace,Nothing}, + fv_old::Union{PVector,Nothing}, + dv_old::Union{AbstractArray,Nothing}, + Uh_old::Union{DistributedSingleFieldFESpace,DistributedMultiFieldFESpace,Nothing}, + model_new, + glue::RedistributeGlue; + reverse=false +) caches = get_redistribute_free_values_cache(fv_new,Uh_new,fv_old,dv_old,Uh_old,model_new,glue;reverse=reverse) return redistribute_free_values!(caches,fv_new,Uh_new,fv_old,dv_old,Uh_old,model_new,glue;reverse=reverse) end -function get_redistribute_free_values_cache(fv_new,Uh_new, - fv_old,dv_old,Uh_old, - model_new,glue::RedistributeGlue; - reverse=false) +function get_redistribute_free_values_cache( + fv_new,Uh_new, + fv_old,dv_old,Uh_old, + model_new,glue::RedistributeGlue; + reverse=false +) T = _get_fe_type(Uh_new,Uh_old) get_redistribute_free_values_cache(T,fv_new,Uh_new,fv_old,dv_old,Uh_old,model_new,glue;reverse=reverse) end -function get_redistribute_free_values_cache(::Type{DistributedSingleFieldFESpace}, - fv_new,Uh_new, - fv_old,dv_old,Uh_old, - model_new,glue::RedistributeGlue; - reverse=false) +function get_redistribute_free_values_cache( + ::Type{DistributedSingleFieldFESpace}, + fv_new,Uh_new, + fv_old,dv_old,Uh_old, + model_new,glue::RedistributeGlue; + reverse=false +) old_parts, new_parts = get_old_and_new_parts(glue,Val(reverse)) cell_dof_values_old = i_am_in(old_parts) ? map(scatter_free_and_dirichlet_values,local_views(Uh_old),local_views(fv_old),dv_old) : nothing cell_dof_ids_new = i_am_in(new_parts) ? map(_get_cell_dof_ids_inner_space, local_views(Uh_new)) : nothing @@ -322,11 +373,13 @@ function get_redistribute_free_values_cache(::Type{DistributedSingleFieldFESpace return caches end -function get_redistribute_free_values_cache(::Type{DistributedMultiFieldFESpace}, - fv_new,Uh_new, - fv_old,dv_old,Uh_old, - model_new,glue::RedistributeGlue; - reverse=false) +function get_redistribute_free_values_cache( + ::Type{DistributedMultiFieldFESpace}, + fv_new,Uh_new, + fv_old,dv_old,Uh_old, + model_new,glue::RedistributeGlue; + reverse=false +) old_parts, new_parts = get_old_and_new_parts(glue,Val(reverse)) if i_am_in(old_parts) @@ -356,27 +409,30 @@ function get_redistribute_free_values_cache(::Type{DistributedMultiFieldFESpace} return caches end -function redistribute_free_values!(caches, - fv_new,Uh_new, - fv_old,dv_old,Uh_old, - model_new, - glue::RedistributeGlue; - reverse=false) +function redistribute_free_values!( + caches, + fv_new,Uh_new, + fv_old,dv_old,Uh_old, + model_new, + glue::RedistributeGlue; + reverse=false +) T = _get_fe_type(Uh_new,Uh_old) redistribute_free_values!(T,caches,fv_new,Uh_new,fv_old,dv_old,Uh_old,model_new,glue;reverse=reverse) end -function redistribute_free_values!(::Type{DistributedSingleFieldFESpace}, - caches, - fv_new::Union{PVector,Nothing}, - Uh_new::Union{DistributedSingleFieldFESpace,Nothing}, - fv_old::Union{PVector,Nothing}, - dv_old::Union{AbstractArray,Nothing}, - Uh_old::Union{DistributedSingleFieldFESpace,Nothing}, - model_new, - glue::RedistributeGlue; - reverse=false) - +function redistribute_free_values!( + ::Type{DistributedSingleFieldFESpace}, + caches, + fv_new::Union{PVector,Nothing}, + Uh_new::Union{DistributedSingleFieldFESpace,Nothing}, + fv_old::Union{PVector,Nothing}, + dv_old::Union{AbstractArray,Nothing}, + Uh_old::Union{DistributedSingleFieldFESpace,Nothing}, + model_new, + glue::RedistributeGlue; + reverse=false +) old_parts, new_parts = get_old_and_new_parts(glue,Val(reverse)) cell_dof_values_old = i_am_in(old_parts) ? map(scatter_free_and_dirichlet_values,local_views(Uh_old),local_views(fv_old),dv_old) : nothing cell_dof_ids_new = i_am_in(new_parts) ? map(_get_cell_dof_ids_inner_space, local_views(Uh_new)) : nothing @@ -389,17 +445,18 @@ function redistribute_free_values!(::Type{DistributedSingleFieldFESpace}, return fv_new end -function redistribute_free_values!(::Type{DistributedMultiFieldFESpace}, - caches, - fv_new::Union{PVector,Nothing}, - Uh_new::Union{DistributedMultiFieldFESpace,Nothing}, - fv_old::Union{PVector,Nothing}, - dv_old::Union{AbstractArray,Nothing}, - Uh_old::Union{DistributedMultiFieldFESpace,Nothing}, - model_new, - glue::RedistributeGlue; - reverse=false) - +function redistribute_free_values!( + ::Type{DistributedMultiFieldFESpace}, + caches, + fv_new::Union{PVector,Nothing}, + Uh_new::Union{DistributedMultiFieldFESpace,Nothing}, + fv_old::Union{PVector,Nothing}, + dv_old::Union{AbstractArray,Nothing}, + Uh_old::Union{DistributedMultiFieldFESpace,Nothing}, + model_new, + glue::RedistributeGlue; + reverse=false +) old_parts, new_parts = get_old_and_new_parts(glue,Val(reverse)) if i_am_in(old_parts) @@ -426,3 +483,257 @@ function redistribute_free_values!(::Type{DistributedMultiFieldFESpace}, redistribute_free_values!(DistributedSingleFieldFESpace,caches,fv_new_i,Uh_new_i,fv_old_i,dv_old_i,Uh_old_i,model_new,glue;reverse=reverse) end end + + +# Cartesian Model uniform refinement + +function Adaptivity.refine( + cmodel::DistributedDiscreteModel{Dc}, + refs::Integer = 2 +) where Dc + Adaptivity.refine(cmodel,Tuple(fill(refs,Dc))) +end + +function Adaptivity.refine( + cmodel::DistributedAdaptedDiscreteModel{Dc}, + refs::NTuple{Dc,<:Integer} +) where Dc + + # Local cmodels are AdaptedDiscreteModels. To correctly dispatch, we need to + # extract the underlying models, then refine. + _cmodel = get_model(cmodel) + _fmodel = refine(_cmodel,refs) + + # Now the issue is that the local parents are not pointing to local_views(cmodel). + # We have to fix that... + fmodel = GenericDistributedDiscreteModel( + map(get_model,local_views(_fmodel)), + get_cell_gids(_fmodel); + metadata=_fmodel.metadata + ) + glues = get_adaptivity_glue(_fmodel) + return DistributedAdaptedDiscreteModel(fmodel,cmodel,glues) +end + +function Adaptivity.refine( + cmodel::DistributedCartesianDiscreteModel{Dc}, + refs::NTuple{Dc,<:Integer}, +) where Dc + + ranks = linear_indices(local_views(cmodel)) + desc, parts = cmodel.metadata.descriptor, cmodel.metadata.mesh_partition + + # Create the new model + nC = desc.partition + domain = Adaptivity._get_cartesian_domain(desc) + nF = nC .* refs + fmodel = CartesianDiscreteModel( + ranks,parts,domain,nF;map=desc.map,isperiodic=desc.isperiodic + ) + + # The idea for the glue is the following: + # For each coarse local model (owned + ghost), we can use the serial code to create + # the glue. However, this glue is NOT fully correct. + # Why? Because all the children belonging to coarse ghost cells are in the glue. This + # is not correct, since we only want to keep the children which are ghosts in the new model. + # To this end, we have to remove the extra fine layers of ghosts from the glue. This we + # can do thanks to how predictable the Cartesian model is. + glues = map(ranks,local_views(cmodel)) do rank,cmodel + # Glue for the local models, of size nC_local .* ref + desc_local = get_cartesian_descriptor(cmodel) + nC_local = desc_local.partition + nF_local = nC_local .* refs + f2c_map, child_map = Adaptivity._create_cartesian_f2c_maps(nC_local,refs) + + # Remove extra fine layers of ghosts + p = Tuple(CartesianIndices(parts)[rank]) + periodic_ghosts = map((isp,np) -> (np==1) ? false : isp, desc.isperiodic, parts) + local_range = map(p,parts,periodic_ghosts,nF_local,refs) do p, np, pg, nFl, refs + has_ghost_start = (np > 1) && ((p != 1) || pg) + has_ghost_stop = (np > 1) && ((p != np) || pg) + # If has coarse ghost layer, remove all fine layers but one at each end + start = 1 + has_ghost_start*(refs-1) + stop = nFl - has_ghost_stop*(refs-1) + return start:stop + end + + _indices = LinearIndices(nF_local)[local_range...] + indices = reshape(_indices,length(_indices)) + f2c_cell_map = f2c_map[indices] + fcell_to_child_id = child_map[indices] + + # Create the glue + faces_map = [(d==Dc) ? f2c_cell_map : Int[] for d in 0:Dc] + poly = (Dc == 2) ? QUAD : HEX + reffe = LagrangianRefFE(Float64,poly,1) + rrules = RefinementRule(reffe,refs) + return AdaptivityGlue(faces_map,fcell_to_child_id,rrules) + end + + return DistributedAdaptedDiscreteModel(fmodel,cmodel,glues) +end + +# Cartesian Model redistribution + +@inline function redistribute( + old_model::Union{DistributedCartesianDiscreteModel,Nothing}, + pdesc::DistributedCartesianDescriptor; + old_ranks = nothing +) + redistribute_cartesian(old_model,pdesc;old_ranks) +end + +""" + redistribute_cartesian(old_model,new_ranks,new_parts) + redistribute_cartesian(old_model,pdesc::DistributedCartesianDescriptor) + + Redistributes a DistributedCartesianDiscreteModel to a new set of ranks and parts. + Only redistributes into a superset of the old_model ranks (i.e. towards more processors). +""" +function redistribute_cartesian( + old_model::Union{DistributedCartesianDiscreteModel,Nothing}, + new_ranks, + new_parts; + old_ranks = nothing +) + _pdesc = isnothing(old_model) ? nothing : old_model.metadata + pdesc = emit_cartesian_descriptor(_pdesc,new_ranks,new_parts) + redistribute_cartesian(old_model,pdesc;old_ranks) +end + +function redistribute_cartesian( + old_model::Union{DistributedCartesianDiscreteModel,Nothing}, + pdesc::DistributedCartesianDescriptor{Dc}; + old_ranks = nothing +) where Dc + new_ranks = pdesc.ranks + new_parts = pdesc.mesh_partition + + desc = pdesc.descriptor + ncells = desc.partition + domain = Adaptivity._get_cartesian_domain(desc) + new_model = CartesianDiscreteModel(new_ranks,new_parts,domain,ncells) + + rglue = get_cartesian_redistribute_glue(new_model,old_model;old_ranks) + return new_model, rglue +end + +function get_cartesian_owners(gids,nparts,ncells) + # This is currently O(sqrt(np)), but I believe we could make it + # O(ln(np)) if we ensured the search is sorted. Even faster if we sort + # the gids first, which progressively reduces the number of ranges to search. + ranges = map(nparts,ncells) do np, nc + map(p -> PartitionedArrays.local_range(p,np,nc,false,false), 1:np) + end + cart_ids = CartesianIndices(ncells) + owner_ids = LinearIndices(nparts) + owners = map(gids) do gid + gid_cart = Tuple(cart_ids[gid]) + owner_cart = map((r,g) -> findfirst(ri -> g ∈ ri,r),ranges,gid_cart) + return owner_ids[owner_cart...] + end + return owners +end + +function get_cartesian_redistribute_glue( + new_model::DistributedCartesianDiscreteModel{Dc}, + old_model::Union{DistributedCartesianDiscreteModel{Dc},Nothing}; + old_ranks = nothing +) where Dc + pdesc = new_model.metadata + desc = pdesc.descriptor + ncells = desc.partition + + # Components in the new partition + new_ranks = pdesc.ranks + new_parts = new_model.metadata.mesh_partition + new_ids = partition(get_cell_gids(new_model)) + new_models = local_views(new_model) + + # Components in the old partition (if present) + _old_parts = map(new_ranks) do r + (r == 1) ? Int[old_model.metadata.mesh_partition...] : Int[] + end + old_parts = Tuple(PartitionedArrays.getany(emit(_old_parts))) + _old_ids = isnothing(old_model) ? nothing : partition(get_cell_gids(old_model)) + old_ids = change_parts(_old_ids,new_ranks) + _old_models = isnothing(old_model) ? nothing : local_views(old_model) + old_models = change_parts(_old_models,new_ranks) + + # Produce the glue components + old2new,new2old,parts_rcv,parts_snd,lids_rcv,lids_snd = map( + new_ranks,new_models,old_models,new_ids,old_ids) do r, new_model, old_model, new_ids, old_ids + + if !isnothing(old_ids) + # If I'm in the old subprocessor, + # - I send all owned old cells that I don't own in the new model. + # - I receive all owned new cells that I don't own in the old model. + old2new = replace(find_local_to_local_map(old_ids,new_ids), -1 => 0) + new2old = replace(find_local_to_local_map(new_ids,old_ids), -1 => 0) + + new_l2o = local_to_own(new_ids) + old_l2o = local_to_own(old_ids) + mask_rcv = map(1:local_length(new_ids)) do new_lid + old_lid = new2old[new_lid] + A = !iszero(new_l2o[new_lid]) # I own this cell in the new model + B = (iszero(old_lid) || iszero(old_l2o[old_lid])) # I don't own this cell in the old model + return A && B + end + ids_rcv = findall(mask_rcv) + + mask_snd = map(1:local_length(old_ids)) do old_lid + new_lid = old2new[old_lid] + A = !iszero(old_l2o[old_lid]) # I own this cell in the old model + B = (iszero(new_lid) || iszero(new_l2o[new_lid])) # I don't own this cell in the new model + return A && B + end + ids_snd = findall(mask_snd) + else + # If I'm not in the old subprocessor, + # - I don't send anything. + # - I receive all my owned cells in the new model. + old2new = Int[] + new2old = fill(0,num_cells(new_model)) + ids_rcv = collect(own_to_local(new_ids)) + ids_snd = Int[] + end + + # When snd/rcv ids have been selected, we need to find their owners and prepare + # the snd/rcv buffers. + + # First, everyone can potentially receive stuff: + to_global!(ids_rcv,new_ids) + ids_rcv_to_part = get_cartesian_owners(ids_rcv,old_parts,ncells) + to_local!(ids_rcv,new_ids) + parts_rcv = unique(ids_rcv_to_part) + lids_rcv = map(parts_rcv) do nbor + ids_rcv[findall(x -> x == nbor, ids_rcv_to_part)] + end + lids_rcv = convert(Vector{Vector{Int}}, lids_rcv) + + # Then, only the old subprocessor can potentially send stuff: + if !isnothing(old_ids) + to_global!(ids_snd,old_ids) + ids_snd_to_part = get_cartesian_owners(ids_snd,new_parts,ncells) + to_local!(ids_snd,old_ids) + parts_snd = unique(ids_snd_to_part) + lids_snd = map(parts_snd) do nbor + ids_snd[findall(x -> x == nbor, ids_snd_to_part)] + end + lids_snd = convert(Vector{Vector{Int}}, lids_snd) + else + parts_snd = Int[] + lids_snd = [Int[]] + end + + return old2new, new2old, parts_rcv, parts_snd, JaggedArray(lids_rcv), JaggedArray(lids_snd) + end |> tuple_of_arrays + + # WARNING: This will fail if compared (===) with get_parts(old_model) + # Do we really require this in the glue? Could we remove the old ranks? + if isnothing(old_ranks) + old_ranks = generate_subparts(new_ranks,prod(old_parts)) + end + + return RedistributeGlue(new_ranks,old_ranks,parts_rcv,parts_snd,lids_rcv,lids_snd,old2new,new2old) +end diff --git a/src/Geometry.jl b/src/Geometry.jl index 3a3b356..3263e42 100644 --- a/src/Geometry.jl +++ b/src/Geometry.jl @@ -140,16 +140,21 @@ end """ """ -struct GenericDistributedDiscreteModel{Dc,Dp,A,B} <: DistributedDiscreteModel{Dc,Dp} +struct GenericDistributedDiscreteModel{Dc,Dp,A,B,C} <: DistributedDiscreteModel{Dc,Dp} models::A face_gids::B + metadata::C function GenericDistributedDiscreteModel( - models::AbstractArray{<:DiscreteModel{Dc,Dp}}, gids::PRange) where {Dc,Dp} - A = typeof(models) + models::AbstractArray{<:DiscreteModel{Dc,Dp}}, + gids::PRange; + metadata = nothing + ) where {Dc,Dp} face_gids=Vector{PRange}(undef,Dc+1) - face_gids[Dc+1]=gids + face_gids[Dc+1] = gids + A = typeof(models) B = typeof(face_gids) - new{Dc,Dp,A,B}(models,face_gids) + C = typeof(metadata) + new{Dc,Dp,A,B,C}(models,face_gids,metadata) end end @@ -186,6 +191,45 @@ function _setup_face_gids!(dmodel::GenericDistributedDiscreteModel{Dc},dim) wher end # CartesianDiscreteModel +struct DistributedCartesianDescriptor{A,B,C} + ranks::A + mesh_partition::B + descriptor::C + function DistributedCartesianDescriptor( + ranks::AbstractArray{<:Integer}, + mesh_partition :: NTuple{Dc,<:Integer}, + descriptor :: CartesianDescriptor{Dc} + ) where Dc + A = typeof(ranks) + B = typeof(mesh_partition) + C = typeof(descriptor) + new{A,B,C}(ranks,mesh_partition,descriptor) + end +end + +function emit_cartesian_descriptor( + pdesc::Union{<:DistributedCartesianDescriptor{Dc},Nothing}, + new_ranks::AbstractArray{<:Integer}, + new_mesh_partition +) where Dc + f(a) = Tuple(PartitionedArrays.getany(emit(a))) + a, b, c, d = map(new_ranks) do rank + if rank == 1 + desc = pdesc.descriptor + @assert desc.map === identity + Float64[desc.origin.data...], Float64[desc.sizes...], Int[desc.partition...], Bool[desc.isperiodic...] + else + Float64[], Float64[], Int[], Bool[] + end + end |> tuple_of_arrays + origin, sizes, partition, isperiodic = VectorValue(f(a)...), f(b), f(c), f(d) + new_desc = CartesianDescriptor(origin,sizes,partition;isperiodic) + return DistributedCartesianDescriptor(new_ranks,new_mesh_partition,new_desc) +end + +const DistributedCartesianDiscreteModel{Dc,Dp,A,B,C} = + GenericDistributedDiscreteModel{Dc,Dp,<:AbstractArray{<:CartesianDiscreteModel},B,<:DistributedCartesianDescriptor} + function Geometry.CartesianDiscreteModel( ranks::AbstractArray{<:Integer}, # Distributed array with the rank IDs parts::NTuple{N,<:Integer}, # Number of ranks (parts) in each direction @@ -211,7 +255,8 @@ function Geometry.CartesianDiscreteModel( CartesianDiscreteModel(desc,cmin,cmax) end gids = PRange(upartition) - return GenericDistributedDiscreteModel(models,gids) + metadata = DistributedCartesianDescriptor(ranks,parts,desc) + return GenericDistributedDiscreteModel(models,gids;metadata) end end @@ -258,7 +303,8 @@ function _cartesian_model_with_periodic_bcs(ranks,parts,desc) CartesianDiscreteModel(_desc,cmin,cmax,remove_boundary) end gids = PRange(global_partition) - return GenericDistributedDiscreteModel(models,gids) + metadata = DistributedCartesianDescriptor(ranks,parts,desc) + return GenericDistributedDiscreteModel(models,gids;metadata) end ## Helpers to partition a serial model diff --git a/test/CartesianAdaptivityTests.jl b/test/CartesianAdaptivityTests.jl new file mode 100644 index 0000000..3c2194c --- /dev/null +++ b/test/CartesianAdaptivityTests.jl @@ -0,0 +1,174 @@ +module CartesianAdaptivityTests +using Test + +using Gridap +using Gridap.Geometry +using Gridap.Adaptivity +using Gridap.FESpaces + +using MPI +using GridapDistributed +using PartitionedArrays + +using GridapDistributed: i_am_in, generate_subparts +using GridapDistributed: find_local_to_local_map +using GridapDistributed: DistributedAdaptedDiscreteModel, redistribute, redistribute_cartesian +using GridapDistributed: RedistributeGlue, redistribute_cell_dofs, redistribute_fe_function, redistribute_free_values + +function are_equal(a1::MPIArray,a2::MPIArray) + same = map(a1,a2) do a1,a2 + a1 ≈ a2 + end + return reduce(&,same,init=true) +end + +function are_equal(a1::DebugArray,a2::DebugArray) + same = map(a1,a2) do a1,a2 + a1 ≈ a2 + end + return reduce(&,same,init=true) +end + +function are_equal(a1::PVector,a2::PVector) + consistent!(a1) |> wait + consistent!(a2) |> wait + are_equal(own_values(a1),own_values(a2)) +end + +function test_redistribution(coarse_ranks, fine_ranks, model, redist_model, redist_glue) + sol(x) = cos(2π*x[1]) + cos(2π*x[2]) # We need this for periodic tests to work + reffe = ReferenceFE(lagrangian,Float64,1) + + if i_am_in(coarse_ranks) + space = FESpace(model,reffe) + u = interpolate(sol,space) + cell_dofs = map(get_cell_dof_values,local_views(u)) + free_values = get_free_dof_values(u) + dir_values = zero_dirichlet_values(space) + else + space = nothing; u = nothing; cell_dofs = nothing; free_values = nothing; dir_values = nothing; + end + + redist_space = FESpace(redist_model,reffe) + redist_u = interpolate(sol,redist_space) + redist_cell_dofs = map(get_cell_dof_values,local_views(redist_u)) + redist_free_values = get_free_dof_values(redist_u) + redist_dir_values = zero_dirichlet_values(redist_space) + + # Redistribute cell values, both ways + tmp_cell_dofs = deepcopy(redist_cell_dofs) + redistribute_cell_dofs(cell_dofs,tmp_cell_dofs,redist_model,redist_glue) + @test are_equal(redist_cell_dofs,tmp_cell_dofs) + + tmp_cell_dofs = i_am_in(coarse_ranks) ? deepcopy(cell_dofs) : nothing + redistribute_cell_dofs(redist_cell_dofs,tmp_cell_dofs,model,redist_glue;reverse=true) + if i_am_in(coarse_ranks) + @test are_equal(cell_dofs,tmp_cell_dofs) + end + + # Redistribute free values, both ways + tmp_free_values = deepcopy(redist_free_values) + redistribute_free_values(tmp_free_values,redist_space,free_values,dir_values,space,redist_model,redist_glue) + @test are_equal(redist_free_values,tmp_free_values) + + tmp_free_values = i_am_in(coarse_ranks) ? deepcopy(free_values) : nothing + redistribute_free_values(tmp_free_values,space,redist_free_values,redist_dir_values,redist_space,model,redist_glue;reverse=true) + if i_am_in(coarse_ranks) + @test are_equal(free_values,tmp_free_values) + end + + return true +end + +function test_adaptivity(ranks,cmodel,fmodel,glue) + if i_am_in(ranks) + sol(x) = cos(2π*x[1]) + cos(2π*x[2]) # We need this for periodic tests to work + order = 3 + qorder = 2*order + reffe = ReferenceFE(lagrangian,Float64,order) + amodel = fmodel + + Ωf = Triangulation(amodel) + dΩf = Measure(Ωf,qorder) + Vf = FESpace(amodel,reffe) + Uf = TrialFESpace(Vf) + uh_fine = interpolate(sol,Vf) + + Ωc = Triangulation(cmodel) + dΩc = Measure(Ωc,qorder) + Vc = FESpace(cmodel,reffe) + Uc = TrialFESpace(Vc) + uh_coarse = interpolate(sol,Vc) + + dΩcf = Measure(Ωc,Ωf,qorder) + + # Coarse to Fine projection + af(u,v) = ∫(u⋅v)*dΩf + lf(v) = ∫(uh_coarse*v)*dΩf + op = AffineFEOperator(af,lf,Uf,Vf) + uh_coarse_to_fine = solve(op) + + eh = uh_fine - uh_coarse_to_fine + @test sum(∫(eh⋅eh)*dΩf) < 1e-6 + + # Fine to Coarse projection + ac(u,v) = ∫(u⋅v)*dΩc + lc(v) = ∫(uh_fine*v)*dΩcf + op = AffineFEOperator(ac,lc,Uc,Vc) + uh_fine_to_coarse = solve(op) + + eh = uh_coarse - uh_fine_to_coarse + @test sum(∫(eh⋅eh)*dΩc) < 1e-6 + end + return true +end + +############################################################################################ + +function main(distribute,ncells,isperiodic) + fine_parts = (2,2) + fine_ranks = distribute(LinearIndices((4,))) + + coarse_parts = (2,1) + coarse_ranks = generate_subparts(fine_ranks,2) + + map_main(fine_ranks) do r + println("CartesianAdaptivityTests: ncells = ", ncells, ", isperiodic = ", isperiodic) + end + + # Create models and glues + if i_am_in(coarse_ranks) + parent = CartesianDiscreteModel(coarse_ranks,coarse_parts,(0,1,0,1),ncells;isperiodic) + child = refine(parent,(2,2)) + coarse_adaptivity_glue = get_adaptivity_glue(child) + else + parent = nothing; child = nothing; coarse_adaptivity_glue = nothing + end + + # Redistribute when you know it's cartesian + redist_parent, redist_glue_parent = redistribute_cartesian(parent,fine_ranks,fine_parts) + + redist_child_1 = refine(redist_parent,(2,2)) + fine_adaptivity_glue = get_adaptivity_glue(redist_child_1) + + # Redistribute by dispatching on the DistributedCartesianDescriptor + pdesc = redist_child_1.metadata + redist_child_2, redist_glue_child = redistribute(child,pdesc) + + # Tests + test_redistribution(coarse_ranks,fine_ranks,parent,redist_parent,redist_glue_parent) + test_redistribution(coarse_ranks,fine_ranks,child,redist_child_2,redist_glue_child) + + test_adaptivity(coarse_ranks,parent,child,coarse_adaptivity_glue) + test_adaptivity(fine_ranks,redist_parent,redist_child_1,fine_adaptivity_glue) + return +end + +function main(distribute) + main(distribute,(8,8),(false,false)) + main(distribute,(8,8),(true,false)) + main(distribute,(8,8),(false,true)) + main(distribute,(8,8),(true,true)) +end + +end # module AdaptivityTests \ No newline at end of file diff --git a/test/TestApp/src/TestApp.jl b/test/TestApp/src/TestApp.jl index fb8c52b..c9bcb2f 100644 --- a/test/TestApp/src/TestApp.jl +++ b/test/TestApp/src/TestApp.jl @@ -12,6 +12,7 @@ module TestApp include("../../HeatEquationTests.jl") include("../../StokesOpenBoundaryTests.jl") include("../../AdaptivityTests.jl") + include("../../CartesianAdaptivityTests.jl") include("../../AdaptivityMultiFieldTests.jl") include("../../BlockSparseMatrixAssemblersTests.jl") end \ No newline at end of file diff --git a/test/mpi/runtests_np4_body.jl b/test/mpi/runtests_np4_body.jl index 573ad03..c86a88a 100644 --- a/test/mpi/runtests_np4_body.jl +++ b/test/mpi/runtests_np4_body.jl @@ -42,6 +42,7 @@ function all_tests(distribute,parts) if prod(parts) == 4 TestApp.AdaptivityTests.main(distribute) + TestApp.CartesianAdaptivityTests.main(distribute) TestApp.AdaptivityMultiFieldTests.main(distribute) PArrays.toc!(t,"Adaptivity") end diff --git a/test/sequential/AdaptivityTests.jl b/test/sequential/AdaptivityTests.jl new file mode 100644 index 0000000..b858f73 --- /dev/null +++ b/test/sequential/AdaptivityTests.jl @@ -0,0 +1,10 @@ +module AdaptivityTestsSeq + +using PartitionedArrays + +include("../CartesianAdaptivityTests.jl") +with_debug() do distribute + CartesianAdaptivityTests.main(distribute) +end + +end # module \ No newline at end of file diff --git a/test/sequential/runtests.jl b/test/sequential/runtests.jl index dc1a6c0..d9285ff 100644 --- a/test/sequential/runtests.jl +++ b/test/sequential/runtests.jl @@ -40,4 +40,8 @@ end include("BlockSparseMatrixAssemblersTests.jl") end +@time @testset "AdaptivityTests" begin + include("AdaptivityTests.jl") +end + end # module