diff --git a/Project.toml b/Project.toml index edf8d4a..02841ff 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridapP4est" uuid = "c2c8e14b-f5fd-423d-9666-1dd9ad120af9" authors = ["Alberto F. Martin "] -version = "0.3.6" +version = "0.3.7" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" @@ -17,8 +17,8 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [compat] ArgParse = "1" FillArrays = "0.8.4, 0.9, 0.10, 0.11, 0.12, 1" -Gridap = "0.17.22, 0.18" -GridapDistributed = "0.3.1, 0.4" +Gridap = "0.18.2" +GridapDistributed = "0.4" MPI = "0.20" P4est_wrapper = "0.2.2" PartitionedArrays = "0.3.3" diff --git a/src/FESpaces.jl b/src/FESpaces.jl index 0a2df8f..93b36c2 100644 --- a/src/FESpaces.jl +++ b/src/FESpaces.jl @@ -1,3 +1,7 @@ + +const LagragianOrNedelec = + Union{Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any},Tuple{Gridap.ReferenceFEs.Nedelec,Any,Any}} + function _generate_ref_constraints(modelH,modelh,cell_reffe) VH=TestFESpace(modelH,cell_reffe) Vh=TestFESpace(modelh,cell_reffe) @@ -31,7 +35,7 @@ end function _build_constraint_coefficients_matrix_in_ref_space(::PXestUniformRefinementRuleType, Dc, - reffe::Tuple{<:Lagrangian,Any,Any}) + reffe::LagragianOrNedelec) cell_polytope = Dc == 2 ? QUAD : HEX basis, reffe_args, reffe_kwargs = reffe cell_reffe = ReferenceFE(cell_polytope, basis, reffe_args...; reffe_kwargs...) @@ -94,7 +98,7 @@ function _fill_face_subface_ldof_to_cell_ldof!(face_subface_ldof_to_cell_ldof, end function _generate_face_subface_ldof_to_cell_ldof(ref_rule::PXestUniformRefinementRuleType, - Df,Dc,reffe::Tuple{<:Lagrangian,Any,Any}) + Df,Dc,reffe::LagragianOrNedelec) cell_polytope = (Dc == 2) ? QUAD : HEX _coarse_faces_to_child_ids = coarse_faces_to_child_ids(ref_rule,Df,Dc) @@ -276,9 +280,6 @@ function coarse_faces_to_child_ids(::PXestVerticalRefinementRuleType,Df,Dc) end end - - - function coarse_faces_to_child_ids(::PXestHorizontalRefinementRuleType,Df,Dc) @assert Dc==3 if (Df==2) @@ -419,7 +420,10 @@ function _restrict_face_dofs_to_face_dim(cell_reffe,Df) first_face = Gridap.ReferenceFEs.get_offset(polytope,Df)+1 face_own_dofs=get_face_own_dofs(cell_reffe) first_face_faces = Gridap.ReferenceFEs.get_faces(polytope)[first_face] - face_dofs_to_face_dim = face_own_dofs[first_face_faces] + face_dofs_to_face_dim = Vector{Vector{Int}}() + for face in first_face_faces + push!(face_dofs_to_face_dim,copy(face_own_dofs[face])) + end touched=Dict{Int,Int}() dofs=Int64[] current=1 @@ -612,7 +616,22 @@ end function get_face_dofs_permutations( reffe::Gridap.ReferenceFEs.GenericRefFE{Gridap.ReferenceFEs.RaviartThomas, Dc},Df::Integer) where Dc first_face = get_offset(get_polytope(reffe),Df) - order = length(get_face_dofs(reffe)[first_face])-1 + order = length(get_face_dofs(reffe)[first_face+1])-1 + nfaces=num_faces(reffe,Df) + if (Df==Dc-1) + facet_polytope = Dc == 2 ? SEGMENT : QUAD + nodes, _ = Gridap.ReferenceFEs.compute_nodes(facet_polytope, [order for i = 1:Dc-1]) + Fill(Gridap.ReferenceFEs._compute_node_permutations(facet_polytope, nodes),nfaces) + elseif (Dc == 3 && Df==1) + nodes, _ = Gridap.ReferenceFEs.compute_nodes(SEGMENT, [order for i = 1:Dc-2]) + Fill(Gridap.ReferenceFEs._compute_node_permutations(SEGMENT, nodes),nfaces) + end +end + +function get_face_dofs_permutations( + reffe::Gridap.ReferenceFEs.GenericRefFE{Gridap.ReferenceFEs.Nedelec, Dc},Df::Integer) where Dc + first_face = get_offset(get_polytope(reffe),Df) + order = length(get_face_dofs(reffe)[first_face+1])-1 nfaces=num_faces(reffe,Df) if (Df==Dc-1) facet_polytope = Dc == 2 ? SEGMENT : QUAD @@ -897,7 +916,7 @@ end # Generates a new DistributedSingleFieldFESpace composed # by local FE spaces with linear multipoint constraints added function Gridap.FESpaces.FESpace(model::OctreeDistributedDiscreteModel{Dc}, - reffe::Tuple{Gridap.ReferenceFEs.Lagrangian,Any,Any}; + reffe::LagragianOrNedelec; kwargs...) where {Dc} spaces_wo_constraints = map(local_views(model)) do m FESpace(m, reffe; kwargs...) diff --git a/test/DarcyNonConformingOctreeModelsTests.jl b/test/DarcyNonConformingOctreeModelsTests.jl index 2d513cf..d1e909c 100644 --- a/test/DarcyNonConformingOctreeModelsTests.jl +++ b/test/DarcyNonConformingOctreeModelsTests.jl @@ -142,7 +142,7 @@ module DarcyNonConformingOctreeModelsTests dΩh = Measure(trian,degree) el2 = sqrt(sum( ∫( e⋅e )*dΩh )) - tol=1e-7 + tol=1e-6 @assert el2 < tol end @@ -247,25 +247,11 @@ module DarcyNonConformingOctreeModelsTests ep_l2 = l2(ep) ep_h1 = h1(ep) - tol = 1.0e-9 + tol = 1.0e-6 @test eu_l2 < tol @test ep_l2 < tol @test ep_h1 < tol - end - - function driver(ranks,coarse_model,order) - model=OctreeDistributedDiscreteModel(ranks,coarse_model,1) - ref_coarse_flags=map(ranks,partition(get_cell_gids(model.dmodel))) do rank,indices - flags=zeros(Cint,length(indices)) - flags.=nothing_flag - #flags[1:2^2].=coarsen_flag - flags[own_length(indices)]=refine_flag - flags - end - fmodel,glue=Gridap.Adaptivity.adapt(model,ref_coarse_flags) - xh,Xh = solve_darcy(fmodel,order) - check_error_darcy(fmodel,order,xh) - end + end function run(distribute) # debug_logger = ConsoleLogger(stderr, Logging.Debug) @@ -276,7 +262,7 @@ module DarcyNonConformingOctreeModelsTests order=1 test_2d(ranks,order) - order=1 + order=2 test_3d(ranks,order) end end \ No newline at end of file diff --git a/test/MaxwellNonConformingOctreeModelsTests.jl b/test/MaxwellNonConformingOctreeModelsTests.jl new file mode 100644 index 0000000..3abc6f8 --- /dev/null +++ b/test/MaxwellNonConformingOctreeModelsTests.jl @@ -0,0 +1,241 @@ +module MaxwellNonConformingOctreeModelsTests + using P4est_wrapper + using GridapP4est + using Gridap + using PartitionedArrays + using GridapDistributed + using MPI + using Gridap.FESpaces + using FillArrays + using Logging + using Test + + function test_transfer_ops_and_redistribute(ranks, + dmodel::GridapDistributed.DistributedDiscreteModel{Dc}, + order) where Dc + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices + flags=zeros(Cint,length(indices)) + flags.=nothing_flag + + flags[1]=refine_flag + flags[own_length(indices)]=refine_flag + + # To create some unbalance + if (rank%2==0 && own_length(indices)>1) + flags[div(own_length(indices),2)]=refine_flag + end + flags + end + fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); + + # Solve coarse + uH,UH=solve_maxwell(dmodel,order) + check_error_maxwell(dmodel,order,uH) + + # Solve fine + uh,Uh=solve_maxwell(fmodel,order) + check_error_maxwell(fmodel,order,uh) + + Ωh = Triangulation(fmodel) + degree = 2*(order+1) + dΩh = Measure(Ωh,degree) + + # prolongation via interpolation + uHh=interpolate(uH,Uh) + e = uh - uHh + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + tol=1e-6 + println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + # prolongation via L2-projection + # Coarse FEFunction -> Fine FEFunction, by projection + ahp(u,v) = ∫(v⋅u)*dΩh + lhp(v) = ∫(v⋅uH)*dΩh + oph = AffineFEOperator(ahp,lhp,Uh,Uh) + uHh = solve(oph) + e = uh - uHh + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + println("[L2 PROJECTION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + # restriction via interpolation + uhH=interpolate(uh,UH) + e = uH - uhH + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + println("[INTERPOLATION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + # restriction via L2-projection + ΩH = Triangulation(dmodel) + degree = 2*(order+1) + dΩH = Measure(ΩH,degree) + + dΩhH = Measure(ΩH,Ωh,2*order) + aHp(u,v) = ∫(v⋅u)*dΩH + lHp(v) = ∫(v⋅uh)*dΩhH + oph = AffineFEOperator(aHp,lHp,UH,UH) + uhH = solve(oph) + e = uH - uhH + el2 = sqrt(sum( ∫( e⋅e )*dΩH )) + + fmodel_red, red_glue=GridapDistributed.redistribute(fmodel); + uh_red,Uh_red=solve_maxwell(fmodel_red,order) + check_error_maxwell(fmodel_red,order,uh_red) + + trian = Triangulation(fmodel_red) + degree = 2*(order+1) + dΩhred = Measure(trian,degree) + + u_ex, f_ex=get_analytical_functions(Dc) + + uhred2 = GridapDistributed.redistribute_fe_function(uh,Uh_red,fmodel_red,red_glue) + + e = u_ex - uhred2 + el2 = sqrt(sum( ∫( e⋅e )*dΩhred )) + println("[REDISTRIBUTE SOLUTION] el2 < tol: $(el2) < $(tol)") + @assert el2 < tol + + fmodel_red + end + + function test_refine_and_coarsen_at_once(ranks, + dmodel::OctreeDistributedDiscreteModel{Dc}, + order) where Dc + degree = 2*order+1 + ref_coarse_flags=map(ranks,partition(get_cell_gids(dmodel.dmodel))) do rank,indices + flags=zeros(Cint,length(indices)) + flags.=nothing_flag + if (rank==1) + flags[1:min(2^Dc,own_length(indices))].=coarsen_flag + end + flags[own_length(indices)]=refine_flag + flags + end + fmodel,glue=Gridap.Adaptivity.adapt(dmodel,ref_coarse_flags); + + # Solve coarse + uH,UH=solve_maxwell(dmodel,order) + check_error_maxwell(dmodel,order,uH) + + # Solve fine + uh,Uh=solve_maxwell(fmodel,order) + check_error_maxwell(fmodel,order,uh) + + # # prolongation via interpolation + uHh=interpolate(uH,Uh) + e = uh - uHh + + trian = Triangulation(fmodel) + degree = 2*(order+1) + dΩh = Measure(trian,degree) + + el2 = sqrt(sum( ∫( e⋅e )*dΩh )) + tol=1e-6 + @assert el2 < tol + end + + function test_2d(ranks,order) + coarse_model=CartesianDiscreteModel((0,1,0,1),(1,1)) + dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,1) + test_refine_and_coarsen_at_once(ranks,dmodel,order) + rdmodel=dmodel + for i=1:5 + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) + end + end + + function test_3d(ranks,order) + coarse_model=CartesianDiscreteModel((0,1,0,1,0,1),(2,2,2)) + dmodel=OctreeDistributedDiscreteModel(ranks,coarse_model,0) + test_refine_and_coarsen_at_once(ranks,dmodel,order) + rdmodel=dmodel + for i=1:5 + rdmodel=test_transfer_ops_and_redistribute(ranks,rdmodel,order) + end + end + + u_ex_2D((x,y)) = 2*VectorValue(-y,x) + f_ex_2D(x) = u_ex_2D(x) + u_ex_3D((x,y,z)) = 2*VectorValue(-y,x,0.) - VectorValue(0.,-z,y) + f_ex_3D(x) = u_ex_3D(x) + + function get_analytical_functions(Dc) + if Dc==2 + return u_ex_2D, f_ex_2D + else + @assert Dc==3 + return u_ex_3D, f_ex_3D + end + end + + include("CoarseDiscreteModelsTools.jl") + + function solve_maxwell(model::GridapDistributed.DistributedDiscreteModel{Dc},order) where {Dc} + u_ex, f_ex=get_analytical_functions(Dc) + + V = FESpace(model, + ReferenceFE(nedelec,order), + conformity=:Hcurl, + dirichlet_tags="boundary") + + U = TrialFESpace(V,u_ex) + + trian = Triangulation(model) + degree = 2*(order+1) + dΩ = Measure(trian,degree) + + a(u,v) = ∫( (∇×u)⋅(∇×v) + u⋅v )dΩ + l(v) = ∫(f_ex⋅v)dΩ + + op = AffineFEOperator(a,l,U,V) + if (num_free_dofs(U)==0) + # UMFPACK cannot handle empty linear systems + uh = zero(U) + else + uh = solve(op) + end + + # uh_ex=interpolate(u_ex_3D,U) + # map(local_views(get_free_dof_values(uh_ex)), local_views(op.op.matrix), local_views(op.op.vector)) do U_ex, A, b + # r_ex = A*U_ex - b + # println(r_ex) + # end + uh,U + end + + function check_error_maxwell(model::GridapDistributed.DistributedDiscreteModel{Dc},order,uh) where {Dc} + trian = Triangulation(model) + degree = 2*(order+1) + dΩ = Measure(trian,degree) + + u_ex, f_ex = get_analytical_functions(Dc) + + eu = u_ex - uh + + l2(v) = sqrt(sum(∫(v⋅v)*dΩ)) + hcurl(v) = sqrt(sum(∫(v⋅v + (∇×v)⋅(∇×v))*dΩ)) + + eu_l2 = l2(eu) + eu_hcurl = hcurl(eu) + + tol = 1.0e-6 + @test eu_l2 < tol + @test eu_hcurl < tol + end + + function run(distribute) + # debug_logger = ConsoleLogger(stderr, Logging.Debug) + # global_logger(debug_logger); # Enable the debug logger globally + np = MPI.Comm_size(MPI.COMM_WORLD) + ranks = distribute(LinearIndices((np,))) + + for order=0:2 + test_2d(ranks,order) + end + + for order=0:2 + test_3d(ranks,order) + end + end +end diff --git a/test/mpi/MaxwellNonConformingOctreeModelsTests.jl b/test/mpi/MaxwellNonConformingOctreeModelsTests.jl new file mode 100644 index 0000000..69d3f81 --- /dev/null +++ b/test/mpi/MaxwellNonConformingOctreeModelsTests.jl @@ -0,0 +1,17 @@ + +using MPI +using PartitionedArrays + +include("../MaxwellNonConformingOctreeModelsTests.jl") +import .MaxwellNonConformingOctreeModelsTests as TestModule + +if !MPI.Initialized() + MPI.Init() +end + +with_mpi() do distribute + TestModule.run(distribute) +end + +MPI.Finalize() + diff --git a/test/runtests.jl b/test/runtests.jl index 28168c1..bb038db 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -44,6 +44,9 @@ function run_tests(testdir) elseif f in ["DarcyNonConformingOctreeModelsTests.jl"] np = [1,4] extra_args = "" + elseif f in ["MaxwellNonConformingOctreeModelsTests.jl"] + np = [1,4] + extra_args = "" elseif f in ["PoissonNonConformingOctreeModelsTests.jl"] np = [1,2,4] extra_args = ""