From b6321f167db3e840d2d7d56c7b8678d0777e0855 Mon Sep 17 00:00:00 2001 From: "Alberto F. Martin" Date: Tue, 5 Mar 2024 22:21:50 +1100 Subject: [PATCH] * Some bug-fixes for 2D(adaptive)+1D(uniform) refinement * Added missing methods for uniformly refining/coarsening in the regular quad/octree --- ...callyAdapted3DDistributedDiscreteModels.jl | 66 +++++++++++-------- src/PXestTypeMethods.jl | 56 +++++++++++++--- test/PoissonAnisotropicOctreeModelsTests.jl | 8 ++- test/runtests.jl | 3 + 4 files changed, 96 insertions(+), 37 deletions(-) diff --git a/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl b/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl index 35a9e21..363a0b5 100644 --- a/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl +++ b/src/AnisotropicallyAdapted3DDistributedDiscreteModels.jl @@ -37,7 +37,7 @@ function AnisotropicallyAdapted3DDistributedDiscreteModel( ptr_pXest_connectivity, ptr_pXest, pXest_type, - nothing, + PXestHorizontalRefinementRuleType(), true, nothing) end @@ -407,11 +407,14 @@ function generate_face_labeling(pXest_type::P6estType, polytope = HEX + + update_face_to_entity_with_ghost_data!(vertex_to_entity, cell_prange, num_faces(polytope,0), cell_to_faces(topology,Dc,0)) + update_face_to_entity_with_ghost_data!(edget_to_entity, cell_prange, num_faces(polytope,1), @@ -421,11 +424,13 @@ function generate_face_labeling(pXest_type::P6estType, cell_prange, num_faces(polytope,Dc-1), cell_to_faces(topology,Dc,Dc-1)) + update_face_to_entity_with_ghost_data!(cell_to_entity, cell_prange, num_faces(polytope,Dc), cell_to_faces(topology,Dc,Dc)) + faces_to_entity=[vertex_to_entity,edget_to_entity,facet_to_entity,cell_to_entity] @@ -476,6 +481,21 @@ face_labeling = face_labeling end + + # map(partition(cell_prange)) do indices + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] l2g=$(local_to_global(indices))") + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] l2o=$(local_to_own(indices))") + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] l2p=$(local_to_owner(indices))") + + + # cache=PartitionedArrays.assembly_cache(indices) + + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] ns=$(cache.neighbors_snd[])") + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] nr=$(cache.neighbors_rcv[])") + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] li_snd=$(cache.local_indices_snd[])") + # println("[$(MPI.Comm_rank(MPI.COMM_WORLD))] li_rcv=$(cache.local_indices_rcv[])") + # end + face_labeling end @@ -507,9 +527,7 @@ function _horizontally_refine_coarsen_balance!(model::OctreeDistributedDiscreteM map(refinement_and_coarsening_flags,num_cols) do flags, num_cols # The length of the local flags array has to match the number of locally owned columns in the model @assert num_cols==length(flags) - println("BBB: $(flags)") pXest_reset_data!(pXest_type, model.ptr_pXest, Cint(sizeof(Cint)), init_fn_callback_c, pointer(flags)) - println("CCC: $(flags)") end @@ -559,16 +577,23 @@ function horizontally_adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, pXest_lnodes_destroy(model.pXest_type,ptr_pXest_lnodes) pXest_refinement_rule_type = PXestHorizontalRefinementRuleType() - _extruded_flags = _extrude_refinement_and_coarsening_flags(_refinement_and_coarsening_flags, - model.ptr_pXest, - ptr_new_pXest) + extruded_ref_coarsen_flags=map(partition(get_cell_gids(model)),refinement_and_coarsening_flags) do indices, flags + similar(flags, length(local_to_global(indices))) + end + + _extrude_refinement_and_coarsening_flags!(extruded_ref_coarsen_flags, + refinement_and_coarsening_flags, + model.ptr_pXest, + ptr_new_pXest) + stride = pXest_stride_among_children(model.pXest_type,pXest_refinement_rule_type,model.ptr_pXest) + adaptivity_glue = _compute_fine_to_coarse_model_glue(model.pXest_type, pXest_refinement_rule_type, model.parts, model.dmodel, fmodel, - _extruded_flags, + extruded_ref_coarsen_flags, stride) adaptive_models = map(local_views(model), local_views(fmodel), @@ -590,10 +615,11 @@ function horizontally_adapt(model::OctreeDistributedDiscreteModel{Dc,Dp}, return ref_model, adaptivity_glue end - function _extrude_refinement_and_coarsening_flags( - refinement_and_coarsening_flags::MPIArray{<:Vector{<:Integer}}, - ptr_pXest_old, - ptr_pXest_new) +function _extrude_refinement_and_coarsening_flags!( + extruded_flags::MPIArray{<:Vector{<:Integer}}, + flags::MPIArray{<:Vector{<:Integer}}, + ptr_pXest_old, + ptr_pXest_new) pXest_old = ptr_pXest_old[] pXest_new = ptr_pXest_new[] @@ -602,19 +628,9 @@ end num_trees = Cint(pXest_old.columns[].connectivity[].num_trees) @assert num_trees == Cint(pXest_new.columns[].connectivity[].num_trees) - map(refinement_and_coarsening_flags) do flags + map(flags,extruded_flags) do flags,extruded_flags current_old_quad=1 - current_cell_old=1 - total=0 - for itree=0:num_trees-1 - tree = pXest_tree_array_index(pXest_type,pXest_old,itree)[] - num_quads = Cint(tree.quadrants.elem_count) - q = pXest_quadrant_array_index(pXest_type,tree,0) - f,l=P6EST_COLUMN_GET_RANGE(q[]) - num_layers=l-f - total+=num_quads*num_layers - end - extruded_flags=similar(flags, total) + current_cell_old=1 # Go over trees for itree=0:num_trees-1 @@ -654,7 +670,5 @@ end end end end - - extruded_flags end -end \ No newline at end of file +end \ No newline at end of file diff --git a/src/PXestTypeMethods.jl b/src/PXestTypeMethods.jl index 97ba91c..7dd3205 100644 --- a/src/PXestTypeMethods.jl +++ b/src/PXestTypeMethods.jl @@ -2409,12 +2409,52 @@ end function pXest_stride_among_children(pXest_type::P6estType, ::PXestHorizontalRefinementRuleType, ptr_pXest) - # Here we are assuming: - # (1) Each processor has at least one column. - # (2) The number of layers in each column is the same within and accross processors. - pXest = ptr_pXest[] - tree = pXest_tree_array_index(pXest_type, pXest, 0)[] - q = pXest_quadrant_array_index(pXest_type, tree, 0) - f,l=P6EST_COLUMN_GET_RANGE(q[]) - return l-f + # Here we are assuming: + # (1) Each processor has at least one column. + # (2) The number of layers in each column is the same within and accross processors. + num_trees = ptr_pXest[].columns[].connectivity[].num_trees + for itree = 0:num_trees-1 + tree = pXest_tree_array_index(pXest_type, ptr_pXest[], itree)[] + if tree.quadrants.elem_count>0 + q = pXest_quadrant_array_index(pXest_type, tree, 0) + f,l=P6EST_COLUMN_GET_RANGE(q[]) + return l-f + end + end + return 0 +end + +function pXest_uniformly_refine!(::P4estType, ptr_pXest) + # Refine callbacks + function refine_fn_2d(::Ptr{p4est_t},which_tree::p4est_topidx_t,quadrant::Ptr{p4est_quadrant_t}) + return Cint(1) + end + refine_fn_c = @cfunction($refine_fn_2d,Cint,(Ptr{p4est_t}, p4est_topidx_t, Ptr{p4est_quadrant_t})) + p4est_refine(ptr_pXest, Cint(0), refine_fn_c, C_NULL) +end + +function pXest_uniformly_refine!(::P8estType, ptr_pXest) + # Refine callbacks + function refine_fn_3d(::Ptr{p8est_t},which_tree::p4est_topidx_t,quadrant::Ptr{p8est_quadrant_t}) + return Cint(1) + end + # C-callable refine callback 3D + refine_fn_c = @cfunction($refine_fn_3d,Cint,(Ptr{p8est_t}, p4est_topidx_t, Ptr{p8est_quadrant_t})) + p8est_refine(ptr_pXest, Cint(0), refine_fn_c, C_NULL) +end + +function pXest_uniformly_coarsen!(::P4estType, ptr_pXest) + function coarsen_fn_2d(::Ptr{p4est_t},::p4est_topidx_t,::Ptr{Ptr{p4est_quadrant_t}}) + return Cint(1) + end + coarsen_fn_c=@cfunction($coarsen_fn_2d,Cint,(Ptr{p4est_t}, p4est_topidx_t, Ptr{Ptr{p4est_quadrant_t}})) + p4est_coarsen(ptr_pXest, Cint(0), coarsen_fn_c, C_NULL) +end + +function pXest_uniformly_coarsen!(::P8estType, ptr_pXest) + function coarsen_fn_3d(::Ptr{p8est_t},::p4est_topidx_t,::Ptr{Ptr{p8est_quadrant_t}}) + return Cint(1) + end + coarsen_fn_c=@cfunction($coarsen_fn_3d,Cint,(Ptr{p8est_t}, p4est_topidx_t, Ptr{Ptr{p8est_quadrant_t}})) + p8est_coarsen(ptr_pXest, Cint(0), coarsen_fn_c, C_NULL) end \ No newline at end of file diff --git a/test/PoissonAnisotropicOctreeModelsTests.jl b/test/PoissonAnisotropicOctreeModelsTests.jl index c480559..16d1c79 100644 --- a/test/PoissonAnisotropicOctreeModelsTests.jl +++ b/test/PoissonAnisotropicOctreeModelsTests.jl @@ -35,6 +35,8 @@ module PoissonAnisotropicOctreeModelsTests UH=TrialFESpace(VH,u) num_local_cols=GridapP4est.num_locally_owned_columns(dmodel) ref_coarse_flags=map(ranks,num_local_cols) do rank,num_local_cols + println("[rank:$(rank)] $(num_local_cols)") + flags=zeros(Cint,num_local_cols) flags.=nothing_flag @@ -237,7 +239,7 @@ module PoissonAnisotropicOctreeModelsTests function test_3d(ranks,order,T::Type;num_amr_steps=5) coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) - dmodel=AnisotropicallyAdapted3DDistributedDiscreteModel(ranks, coarse_model, 1, 1) + dmodel=AnisotropicallyAdapted3DDistributedDiscreteModel(ranks, coarse_model, 2, 2) #writevtk(dmodel.dmodel,"model") #test_refine_and_coarsen_at_once(ranks,dmodel,order,T) rdmodel=dmodel @@ -261,8 +263,8 @@ module PoissonAnisotropicOctreeModelsTests end end function run(distribute) - debug_logger = ConsoleLogger(stderr, Logging.Debug) - global_logger(debug_logger); # Enable the debug logger globally + # debug_logger = ConsoleLogger(stderr, Logging.Debug) + # global_logger(debug_logger); # Enable the debug logger globally ranks = distribute(LinearIndices((MPI.Comm_size(MPI.COMM_WORLD),))) for perm=1:4, order=1:4, scalar_or_vector in (:scalar,) test(ranks,perm,order,_field_type(Val{3}(),scalar_or_vector)) diff --git a/test/runtests.jl b/test/runtests.jl index 8a3f799..e8819a1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -47,6 +47,9 @@ function run_tests(testdir) elseif f in ["PoissonNonConformingOctreeModelsTests.jl"] np = [1,2,4] extra_args = "" + elseif f in ["PoissonAnisotropicOctreeModelsTests.jl"] + np = [1,4] + extra_args = "" else np = [nprocs] extra_args = ""