diff --git a/src/LinearSolvers/GMGLinearSolvers.jl b/src/LinearSolvers/GMGLinearSolvers.jl index 744cc47b..f38600f5 100644 --- a/src/LinearSolvers/GMGLinearSolvers.jl +++ b/src/LinearSolvers/GMGLinearSolvers.jl @@ -68,7 +68,8 @@ struct GMGNumericalSetup{A,B,C,D,E} <: Gridap.Algebra.NumericalSetup else post_smoothers_caches = pre_smoothers_caches end - coarsest_solver_cache = setup_coarsest_solver_cache(mh,coarsest_solver,smatrices) + #coarsest_solver_cache = setup_coarsest_solver_cache(mh,coarsest_solver,smatrices) + coarsest_solver_cache = coarse_solver_caches(mh,coarsest_solver,smatrices) A = typeof(finest_level_cache) B = typeof(pre_smoothers_caches) @@ -109,6 +110,17 @@ function setup_smoothers_caches(mh::ModelHierarchy,smoothers::AbstractVector{<:L return caches end +function coarse_solver_caches(mh,s,mats) + cache = nothing + nlevs = num_levels(mh) + parts = get_level_parts(mh,nlevs) + if i_am_in(parts) + mat = mats[nlevs] + cache = numerical_setup(symbolic_setup(s, mat), mat) + end + return cache +end + function setup_coarsest_solver_cache(mh::ModelHierarchy,coarsest_solver::LinearSolver,smatrices::Vector{<:AbstractMatrix}) cache = nothing nlevs = num_levels(mh) @@ -214,9 +226,10 @@ function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVec if i_am_in(parts) if (lev == num_levels(mh)) ## Coarsest level - coarsest_solver = ns.solver.coarsest_solver - coarsest_solver_cache = ns.coarsest_solver_cache - solve_coarsest_level!(parts,coarsest_solver,xh,rh,coarsest_solver_cache) + #coarsest_solver = ns.solver.coarsest_solver + #coarsest_solver_cache = ns.coarsest_solver_cache + #solve_coarsest_level!(parts,coarsest_solver,xh,rh,coarsest_solver_cache) + solve!(xh, ns.coarsest_solver_cache, rh) else ## General case Ah = ns.solver.smatrices[lev] diff --git a/src/LinearSolvers/RichardsonSmoothers.jl b/src/LinearSolvers/RichardsonSmoothers.jl index 86cfb176..b3832f54 100644 --- a/src/LinearSolvers/RichardsonSmoothers.jl +++ b/src/LinearSolvers/RichardsonSmoothers.jl @@ -48,12 +48,12 @@ function Gridap.Algebra.solve!(x::AbstractVector,ns::RichardsonSmootherNumerical iter = 1 while iter <= ns.smoother.num_smooth_steps - solve!(dx,Mns,r) - dx .= ns.smoother.damping_factor .* dx - x .= x .+ dx - mul!(Adx, ns.A, dx) - r .= r .- Adx - iter += 1 + solve!(dx,Mns,r) + dx .= ns.smoother.damping_factor .* dx + x .= x .+ dx + mul!(Adx, ns.A, dx) + r .= r .- Adx + iter += 1 end end diff --git a/src/MultilevelTools/Algebra.jl b/src/MultilevelTools/Algebra.jl deleted file mode 100644 index 013ee4f9..00000000 --- a/src/MultilevelTools/Algebra.jl +++ /dev/null @@ -1,29 +0,0 @@ - -# Row/Col vector allocations for serial -function allocate_row_vector(A::AbstractMatrix{T}) where T - return zeros(T,size(A,1)) -end - -function allocate_col_vector(A::AbstractMatrix{T}) where T - return zeros(T,size(A,2)) -end - -# Row/Col vector allocations for parallel -function allocate_row_vector(A::PSparseMatrix) - T = eltype(A) - return pfill(zero(T),partition(axes(A,1))) -end - -function allocate_col_vector(A::PSparseMatrix) - T = eltype(A) - return pfill(zero(T),partition(axes(A,2))) -end - -# Row/Col vector allocations for blocks -function allocate_row_vector(A::AbstractBlockMatrix) - return mortar(map(Aii->allocate_row_vector(Aii),blocks(A)[:,1])) -end - -function allocate_col_vector(A::AbstractBlockMatrix) - return mortar(map(Aii->allocate_col_vector(Aii),blocks(A)[1,:])) -end diff --git a/src/MultilevelTools/DistributedGridTransferOperators.jl b/src/MultilevelTools/DistributedGridTransferOperators.jl index 9e882c4e..84f492e8 100644 --- a/src/MultilevelTools/DistributedGridTransferOperators.jl +++ b/src/MultilevelTools/DistributedGridTransferOperators.jl @@ -25,7 +25,8 @@ function ProlongationOperator(lev::Int,sh::FESpaceHierarchy,qdegree::Int;kwargs. end function DistributedGridTransferOperator(lev::Int,sh::FESpaceHierarchy,qdegree::Int,op_type::Symbol; - mode::Symbol=:solution,restriction_method::Symbol=:projection) + mode::Symbol=:solution,restriction_method::Symbol=:projection, + solver=LUSolver()) mh = sh.mh @check lev < num_levels(mh) @check op_type ∈ [:restriction, :prolongation] @@ -35,13 +36,16 @@ function DistributedGridTransferOperator(lev::Int,sh::FESpaceHierarchy,qdegree:: # Refinement if (op_type == :prolongation) || (restriction_method ∈ [:interpolation,:dof_mask]) cache_refine = _get_interpolation_cache(lev,sh,qdegree,mode) - else + elseif mode == :solution cache_refine = _get_projection_cache(lev,sh,qdegree,mode) + else + cache_refine = _get_dual_projection_cache(lev,sh,qdegree,solver) + restriction_method = :dual_projection end # Redistribution redist = has_redistribution(mh,lev) - cache_redist = _get_redistribution_cache(lev,sh,mode,op_type,cache_refine) + cache_redist = _get_redistribution_cache(lev,sh,mode,op_type,restriction_method,cache_refine) cache = cache_refine, cache_redist return DistributedGridTransferOperator(op_type,redist,restriction_method,sh,cache) @@ -115,7 +119,38 @@ function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode:: return cache_refine end -function _get_redistribution_cache(lev::Int,sh::FESpaceHierarchy,mode::Symbol,op_type::Symbol,cache_refine) +function _get_dual_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,solver) + mh = sh.mh + cparts = get_level_parts(mh,lev+1) + + if i_am_in(cparts) + model_h = get_model_before_redist(mh,lev) + Uh = get_fe_space_before_redist(sh,lev) + Ωh = Triangulation(model_h) + dΩh = Measure(Ωh,qdegree) + uh = FEFunction(Uh,zero_free_values(Uh),zero_dirichlet_values(Uh)) + + model_H = get_model(mh,lev+1) + UH = get_fe_space(sh,lev+1) + ΩH = Triangulation(model_H) + dΩhH = Measure(ΩH,Ωh,qdegree) + + Mh = assemble_matrix((u,v)->∫(v⋅u)*dΩh,Uh,Uh) + Mh_ns = numerical_setup(symbolic_setup(solver,Mh),Mh) + + assem = SparseMatrixAssembler(UH,UH) + rh = allocate_col_vector(Mh) + cache_refine = model_h, Uh, UH, Mh_ns, rh, uh, assem, dΩhH + else + model_h = get_model_before_redist(mh,lev) + Uh = get_fe_space_before_redist(sh,lev) + cache_refine = model_h, Uh, nothing, nothing, nothing, nothing, nothing, nothing + end + + return cache_refine +end + +function _get_redistribution_cache(lev::Int,sh::FESpaceHierarchy,mode::Symbol,op_type::Symbol,restriction_method::Symbol,cache_refine) mh = sh.mh redist = has_redistribution(mh,lev) if !redist @@ -130,10 +165,14 @@ function _get_redistribution_cache(lev::Int,sh::FESpaceHierarchy,mode::Symbol,op glue = mh.levels[lev].red_glue if op_type == :prolongation - model_h, Uh, fv_h, dv_h = cache_refine + model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine cache_exchange = get_redistribute_free_values_cache(fv_h_red,Uh_red,fv_h,dv_h,Uh,model_h_red,glue;reverse=false) + elseif restriction_method == :projection + model_h, Uh, fv_h, dv_h, VH, AH, lH, xH, bH, bH0, assem = cache_refine + cache_exchange = get_redistribute_free_values_cache(fv_h,Uh,fv_h_red,dv_h_red,Uh_red,model_h,glue;reverse=true) else - model_h, Uh, fv_h, dv_h = cache_refine + model_h, Uh, UH, Mh, rh, uh, assem, dΩhH = cache_refine + fv_h = isa(uh,Nothing) ? nothing : get_free_dof_values(uh) cache_exchange = get_redistribute_free_values_cache(fv_h,Uh,fv_h_red,dv_h_red,Uh_red,model_h,glue;reverse=true) end @@ -307,4 +346,43 @@ function LinearAlgebra.mul!(y::Union{PVector,Nothing},A::DistributedGridTransfer end return y +end + +############################################################### + +function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dual_projection}},x::PVector) + cache_refine, cache_redist = A.cache + model_h, Uh, VH, Mh_ns, rh, uh, assem, dΩhH = cache_refine + fv_h = get_free_dof_values(uh) + + solve!(rh,Mh_ns,x) + copy!(fv_h,rh) + consistent!(fv_h) |> fetch + v = get_fe_basis(VH) + assemble_vector!(y,assem,collect_cell_vector(VH,∫(v⋅uh)*dΩhH)) + + return y +end + +function LinearAlgebra.mul!(y::Union{PVector,Nothing},A::DistributedGridTransferOperator{Val{:restriction},Val{true},Val{:dual_projection}},x::PVector) + cache_refine, cache_redist = A.cache + model_h, Uh, VH, Mh_ns, rh, uh, assem, dΩhH = cache_refine + fv_h_red, dv_h_red, Uh_red, model_h_red, glue, cache_exchange = cache_redist + fv_h = isa(uh,Nothing) ? nothing : get_free_dof_values(uh) + + # 1 - Redistribute from fine partition to coarse partition + copy!(fv_h_red,x) + consistent!(fv_h_red) |> fetch + redistribute_free_values!(cache_exchange,fv_h,Uh,fv_h_red,dv_h_red,Uh_red,model_h,glue;reverse=true) + + # 2 - Solve f2c projection coarse partition + if !isa(y,Nothing) + solve!(rh,Mh_ns,fv_h) + copy!(fv_h,rh) + consistent!(fv_h) |> fetch + v = get_fe_basis(VH) + assemble_vector!(y,assem,collect_cell_vector(VH,∫(v⋅uh)*dΩhH)) + end + + return y end \ No newline at end of file diff --git a/src/MultilevelTools/MultilevelTools.jl b/src/MultilevelTools/MultilevelTools.jl index f71e6e30..98f6bed8 100644 --- a/src/MultilevelTools/MultilevelTools.jl +++ b/src/MultilevelTools/MultilevelTools.jl @@ -21,7 +21,6 @@ using GridapDistributed: redistribute_fe_function using GridapDistributed: get_old_and_new_parts using GridapDistributed: generate_subparts, local_views -export allocate_col_vector, allocate_row_vector export change_parts, num_parts, i_am_in export generate_level_parts, generate_subparts @@ -38,7 +37,6 @@ export RestrictionOperator, ProlongationOperator export setup_transfer_operators export mul! -include("Algebra.jl") include("SubpartitioningTools.jl") include("GridapFixes.jl") include("RefinementTools.jl") diff --git a/test/LinearSolvers/BlockDiagonalSmoothersTests.jl b/test/LinearSolvers/BlockDiagonalSmoothersTests.jl index edd87475..bd73281e 100644 --- a/test/LinearSolvers/BlockDiagonalSmoothersTests.jl +++ b/test/LinearSolvers/BlockDiagonalSmoothersTests.jl @@ -118,7 +118,7 @@ function main_driver(D,model,solvers) BDSss = symbolic_setup(BDS,A) BDSns = numerical_setup(BDSss,A) - x = GridapSolvers.allocate_col_vector(A) + x = allocate_col_vector(A) x = cg!(x,A,b;verbose=true,Pl=BDSns,reltol=1.0e-12) @test is_same_vector(x,x_star,Xb,X) @@ -127,7 +127,7 @@ function main_driver(D,model,solvers) BDSss = symbolic_setup(BDS,A) BDSns = numerical_setup(BDSss,A) - x = GridapSolvers.allocate_col_vector(A) + x = allocate_col_vector(A) x = cg!(x,A,b;verbose=true,Pl=BDSns,reltol=1.0e-12) @test is_same_vector(x,x_star,Xb,X) end diff --git a/test/LinearSolvers/GMGTests.jl b/test/LinearSolvers/GMGTests.jl index fcdc0144..6666ea4f 100644 --- a/test/LinearSolvers/GMGTests.jl +++ b/test/LinearSolvers/GMGTests.jl @@ -31,7 +31,7 @@ function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdeg Ω = Triangulation(PD) dΩ = Measure(Ω,qdegree) a(u,v) = biform(u,v,dΩ) - local_solver = LUSolver() # IS_ConjugateGradientSolver(;reltol=1.e-6) + local_solver = IS_ConjugateGradientSolver(;reltol=1.e-6) patch_smoother = PatchBasedLinearSolver(a,Ph,Vh,local_solver) smoothers[lev] = RichardsonSmoother(patch_smoother,10,0.2) end @@ -39,7 +39,7 @@ function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdeg return smoothers end -function gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) +function gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u) tests, trials = spaces tic!(t;barrier=true) @@ -48,7 +48,11 @@ function gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restrict # Preconditioner coarse_solver = LUSolver() - restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual,restriction_method=restriction_method) + restrictions, prolongations = setup_transfer_operators(trials, + qdegree; + mode=:residual, + restriction_method=:projection, + solver=IS_ConjugateGradientSolver(;reltol=1.e-6)) gmg = GMGLinearSolver(mh, smatrices, prolongations, @@ -88,7 +92,7 @@ function gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restrict return e_l2 end -function gmg_poisson_driver(t,parts,mh,order,restriction_method) +function gmg_poisson_driver(t,parts,mh,order) tic!(t;barrier=true) u(x) = x[1] + x[2] f(x) = -Δ(u)(x) @@ -103,14 +107,14 @@ function gmg_poisson_driver(t,parts,mh,order,restriction_method) spaces = tests, trials toc!(t,"FESpaces") - return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u) end -function gmg_laplace_driver(t,parts,mh,order,restriction_method) +function gmg_laplace_driver(t,parts,mh,order) tic!(t;barrier=true) α = 1.0 u(x) = x[1] + x[2] - f(x) = -Δ(u)(x) + f(x) = u(x) - α * Δ(u)(x) biform(u,v,dΩ) = ∫(v*u)dΩ + ∫(α*∇(v)⋅∇(u))dΩ liform(v,dΩ) = ∫(v*f)dΩ qdegree = 2*order+1 @@ -122,18 +126,19 @@ function gmg_laplace_driver(t,parts,mh,order,restriction_method) spaces = tests, trials toc!(t,"FESpaces") - return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u) end -function gmg_vector_laplace_driver(t,parts,mh,order,restriction_method) +function gmg_vector_laplace_driver(t,parts,mh,order) tic!(t;barrier=true) + Dc = num_cell_dims(get_model(mh,1)) α = 1.0 - u(x) = VectorValue(x[1],x[2]) - f(x) = VectorValue(2.0*x[2]*(1.0-x[1]*x[1]),2.0*x[1]*(1-x[2]*x[2])) + u(x) = (Dc==2) ? VectorValue(x[1],x[2]) : VectorValue(x[1],x[2],x[3]) + f(x) = (Dc==2) ? VectorValue(x[1],x[2]) : VectorValue(x[1],x[2],x[3]) biform(u,v,dΩ) = ∫(v⋅u)dΩ + ∫(α*∇(v)⊙∇(u))dΩ liform(v,dΩ) = ∫(v⋅f)dΩ qdegree = 2*order+1 - reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) + reffe = ReferenceFE(lagrangian,VectorValue{Dc,Float64},order) smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0),num_levels(mh)-1) tests = TestFESpace(mh,reffe,dirichlet_tags="boundary") @@ -141,14 +146,15 @@ function gmg_vector_laplace_driver(t,parts,mh,order,restriction_method) spaces = tests, trials toc!(t,"FESpaces") - return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u) end -function gmg_hdiv_driver(t,parts,mh,order,restriction_method) +function gmg_hdiv_driver(t,parts,mh,order) tic!(t;barrier=true) - α = 1.0 - u(x) = VectorValue(x[1],x[2]) - f(x) = VectorValue(2.0*x[2]*(1.0-x[1]*x[1]),2.0*x[1]*(1-x[2]*x[2])) + Dc = num_cell_dims(get_model(mh,1)) + α = 1.0 + u(x) = (Dc==2) ? VectorValue(x[1],x[2]) : VectorValue(x[1],x[2],x[3]) + f(x) = (Dc==2) ? VectorValue(x[1],x[2]) : VectorValue(x[1],x[2],x[3]) biform(u,v,dΩ) = ∫(v⋅u)dΩ + ∫(α*divergence(v)⋅divergence(u))dΩ liform(v,dΩ) = ∫(v⋅f)dΩ qdegree = 2*(order+1) @@ -165,31 +171,26 @@ function gmg_hdiv_driver(t,parts,mh,order,restriction_method) smoothers = get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree) toc!(t,"Patch Decomposition") - return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u) end -function main_gmg_driver(parts,mh,order,restriction_method,pde) +function main_gmg_driver(parts,mh,order,pde) t = PTimer(parts,verbose=true) if pde == :poisson - gmg_poisson_driver(t,parts,mh,order,restriction_method) + gmg_poisson_driver(t,parts,mh,order) elseif pde == :laplace - gmg_laplace_driver(t,parts,mh,order,restriction_method) + gmg_laplace_driver(t,parts,mh,order) elseif pde == :vector_laplace - gmg_vector_laplace_driver(t,parts,mh,order,restriction_method) + gmg_vector_laplace_driver(t,parts,mh,order) elseif pde == :hdiv - gmg_hdiv_driver(t,parts,mh,order,restriction_method) + gmg_hdiv_driver(t,parts,mh,order) end end -function get_mesh_hierarchy(parts,Dc,np_per_level,num_refs_coarse) - if Dc == 2 - domain = (0,1,0,1) - nc = (2,2) - else - @assert Dc == 3 - domain = (0,1,0,1,0,1) - nc = (2,2,2) - end +function get_mesh_hierarchy(parts,nc,np_per_level) + Dc = length(nc) + domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) + num_refs_coarse = (Dc == 2) ? 1 : 0 num_levels = length(np_per_level) cparts = generate_subparts(parts,np_per_level[num_levels]) @@ -199,27 +200,18 @@ function get_mesh_hierarchy(parts,Dc,np_per_level,num_refs_coarse) return mh end -function main(distribute,np,Dc,np_per_level) - parts = distribute(LinearIndices((prod(np),))) - - num_refs_coarse = 2 - mh = get_mesh_hierarchy(parts,Dc,np_per_level,num_refs_coarse) +function main(distribute,np::Integer,nc::Tuple,np_per_level::Vector) + parts = distribute(LinearIndices((np,))) + mh = get_mesh_hierarchy(parts,nc,np_per_level) for pde in [:poisson,:laplace,:vector_laplace,:hdiv] - methods = (pde !== :hdiv) ? [:projection,:interpolation] : [:projection] - for restriction_method in methods - if i_am_main(parts) - println(repeat("=",80)) - println("Testing GMG with Dc=$Dc, PDE=$pde and restriction_method=$restriction_method") - end - order = (pde !== :hdiv) ? 1 : 0 - main_gmg_driver(parts,mh,order,restriction_method,pde) + if i_am_main(parts) + println(repeat("=",80)) + println("Testing GMG with Dc=$(length(nc)), PDE=$pde") end + order = (pde !== :hdiv) ? 1 : 0 + main_gmg_driver(parts,mh,order,pde) end end -with_mpi() do distribute - main(distribute,4,2,[4,2,1]) -end - end # module GMGTests \ No newline at end of file diff --git a/test/LinearSolvers/IterativeSolversWrappersTests.jl b/test/LinearSolvers/IterativeSolversWrappersTests.jl index 6c42782d..9fba4611 100644 --- a/test/LinearSolvers/IterativeSolversWrappersTests.jl +++ b/test/LinearSolvers/IterativeSolversWrappersTests.jl @@ -7,6 +7,7 @@ using LinearAlgebra using SparseArrays using PartitionedArrays +using GridapDistributed using GridapSolvers using GridapSolvers.LinearSolvers @@ -17,7 +18,7 @@ function test_solver(solver,op,Uh,dΩ) A, b = get_matrix(op), get_vector(op); ns = numerical_setup(symbolic_setup(solver,A),A) - x = LinearSolvers.allocate_col_vector(A) + x = allocate_col_vector(A) solve!(x,ns,b) u = interpolate(sol,Uh) diff --git a/test/LinearSolvers/KrylovSolversTests.jl b/test/LinearSolvers/KrylovSolversTests.jl index 9bf8fe51..31c0abb4 100644 --- a/test/LinearSolvers/KrylovSolversTests.jl +++ b/test/LinearSolvers/KrylovSolversTests.jl @@ -16,7 +16,7 @@ function test_solver(solver,op,Uh,dΩ) A, b = get_matrix(op), get_vector(op); ns = numerical_setup(symbolic_setup(solver,A),A) - x = LinearSolvers.allocate_col_vector(A) + x = allocate_col_vector(A) solve!(x,ns,b) u = interpolate(sol,Uh) diff --git a/test/LinearSolvers/SchurComplementSolversTests.jl b/test/LinearSolvers/SchurComplementSolversTests.jl index ca69f078..2e9b4e2c 100644 --- a/test/LinearSolvers/SchurComplementSolversTests.jl +++ b/test/LinearSolvers/SchurComplementSolversTests.jl @@ -110,7 +110,7 @@ function main(distribute,np) gmres = GMRESSolver(20;Pr=psc_solver,rtol=1.e-10,verbose=i_am_main(parts)) gmres_ns = numerical_setup(symbolic_setup(gmres,sysmat),sysmat) - x = LinearSolvers.allocate_col_vector(sysmat) + x = allocate_col_vector(sysmat) solve!(x,gmres_ns,sysvec) xh = FEFunction(X,x) diff --git a/test/LinearSolvers/SmoothersTests.jl b/test/LinearSolvers/SmoothersTests.jl index d5a8be02..15b94301 100644 --- a/test/LinearSolvers/SmoothersTests.jl +++ b/test/LinearSolvers/SmoothersTests.jl @@ -32,7 +32,7 @@ function smoothers_driver(parts,model,P) ss = symbolic_setup(P,A) ns = numerical_setup(ss,A) - x = LinearSolvers.allocate_col_vector(A) + x = allocate_col_vector(A) x, history = IterativeSolvers.cg!(x,A,b; verbose=i_am_main(parts), reltol=1.0e-8, diff --git a/test/LinearSolvers/mpi/GMGTests.jl b/test/LinearSolvers/mpi/GMGTests.jl index eefff1c7..ebc07f30 100644 --- a/test/LinearSolvers/mpi/GMGTests.jl +++ b/test/LinearSolvers/mpi/GMGTests.jl @@ -3,8 +3,8 @@ using MPI, PartitionedArrays include("../GMGTests.jl") with_mpi() do distribute - GMGTests.main(distribute,4,2,[4,2,1]) # 2D - # GMGTests.main(distribute,4,3,[4,2,1]) # 3D + GMGTests.main(distribute,4,(2,2),[4,2,1]) # 2D + GMGTests.main(distribute,4,(2,2,2),[4,2,1]) # 3D end end \ No newline at end of file diff --git a/test/_dev/GMG/GMGDebug.jl b/test/_dev/GMG/GMGDebug.jl index d6f60cb1..9d77691b 100644 --- a/test/_dev/GMG/GMGDebug.jl +++ b/test/_dev/GMG/GMGDebug.jl @@ -64,6 +64,24 @@ xH_star = get_free_dof_values(solve(opH)) Ah, bh = get_matrix(oph), get_vector(oph); AH, bH = get_matrix(opH), get_vector(opH); + +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Vh,Vh) + +function compute_MhH() + MhH = zeros(num_free_dofs(Vh),num_free_dofs(VH)) + xHi = fill(0.0,num_free_dofs(VH)) + for iH in 1:num_free_dofs(VH) + fill!(xHi,0.0); xHi[iH] = 1.0 + vHi = FEFunction(VH,xHi) + vH = assemble_vector((v)->∫(v⋅vHi)*dΩh,Vh) + MhH[:,iH] .= vH + end + return MhH +end + +MhH = compute_MhH() + + # Projection operators function Λ_project(xh) @@ -180,9 +198,9 @@ function smooth!(x,r) niter = 10 for i in 1:niter _r = bh - Λ_project(x) - prolongate!(rp,Ph,_r,w,w_sums) + PatchBasedSmoothers.prolongate!(rp,Ph,r,w,w_sums) dxp = Ap\rp - inject!(dx,Ph,dxp,w,w_sums) + PatchBasedSmoothers.inject!(dx,Ph,dxp,w,w_sums) x .+= β*dx r .-= β*A*dx @@ -196,6 +214,12 @@ niters = 10 wH = randn(size(AH,2)) wh = project_c2f(wH) +function project_f2c_bis(rh) + Qrh = Mhh\rh + uh = FEFunction(Vh,Qrh) + assemble_vector(v->∫(v⋅uh)*dΩHh,VH) +end + iter = 0 error = norm(bh - Ah*xh) while iter < niters && error > 1.0e-10 @@ -210,7 +234,14 @@ while iter < niters && error > 1.0e-10 println(" > norm(xh) = ",norm(xh)) println(" > norm(rh) = ",norm(rh)) - rH = project_f2c(rh) + #rH = project_f2c(rh) + #xH = project_f2c(xh) + #rH = bH - AH*xH + + #Qrh = Mhh\rh + #rH = transpose(MhH)*Qrh + rH = project_f2c_bis(rh) + qH = AH\rH qh = project_c2f(qH) @@ -264,3 +295,16 @@ end using IterativeSolvers x = zeros(size(Ah,2)) cg!(x,Ah,bh;Pl=smoother_ns.Mns,verbose=true) + + + +###################################################################################### + +_s = IS_ConjugateGradientSolver(maxiter=50,reltol=1.e-16,verbose=true) +_ns = numerical_setup(symbolic_setup(_s,Mhh),Mhh) + +x = zeros(size(Mhh,2)) +_rh = copy(rh) +solve!(x,_ns,_rh) + +norm(rh - Mhh*x) diff --git a/test/_dev/GMG/GMGLinearSolversHDivRTTests.jl b/test/_dev/GMG/GMGLinearSolversHDivRTTests.jl index f5a949f7..16648860 100644 --- a/test/_dev/GMG/GMGLinearSolversHDivRTTests.jl +++ b/test/_dev/GMG/GMGLinearSolversHDivRTTests.jl @@ -76,7 +76,10 @@ function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, # Preconditioner tic!(t;barrier=true) smoothers = get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree) - restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual)#,restriction_method=:interpolation) + restrictions, prolongations = setup_transfer_operators(trials, + qdegree; + mode=:residual, + solver=IS_ConjugateGradientSolver(reltol=1.e-6)) gmg = GMGLinearSolver(mh, smatrices, diff --git a/test/_dev/GMG/GMGPatchBasedTesting.jl b/test/_dev/GMG/GMGPatchBasedTesting.jl index 703f19c8..8e99623c 100644 --- a/test/_dev/GMG/GMGPatchBasedTesting.jl +++ b/test/_dev/GMG/GMGPatchBasedTesting.jl @@ -28,8 +28,8 @@ end function test_solver(s,D_j) ns = numerical_setup(symbolic_setup(s,D_j),D_j) - b = GridapSolvers.allocate_col_vector(D_j) - x = GridapSolvers.allocate_col_vector(D_j) + b = allocate_col_vector(D_j) + x = allocate_col_vector(D_j) fill!(b,1.0) solve!(x,ns,b) @@ -40,9 +40,9 @@ end function test_smoother(s,D_j) ns = numerical_setup(symbolic_setup(s,D_j),D_j) - b = GridapSolvers.allocate_col_vector(D_j) - x = GridapSolvers.allocate_col_vector(D_j) - r = GridapSolvers.allocate_row_vector(D_j) + b = allocate_col_vector(D_j) + x = allocate_col_vector(D_j) + r = allocate_row_vector(D_j) fill!(b,1.0) fill!(x,1.0) mul!(r,D_j,x) @@ -144,7 +144,7 @@ gmg = GMGLinearSolver(mh, solver = FGMRESSolver(100,gmg;rtol=1e-6,verbose=true) ns = numerical_setup(symbolic_setup(solver,A),A) -x = GridapSolvers.allocate_col_vector(A) +x = allocate_col_vector(A) solve!(x,ns,b) @@ -154,5 +154,5 @@ test_smoother(smoothers[1],A) Pl = LinearSolvers.IdentitySolver() solver2 = GMRESSolver(1000;Pl=Pl,rtol=1e-6,verbose=true) ns2 = numerical_setup(symbolic_setup(solver2,A),A) -x2 = GridapSolvers.allocate_col_vector(A) +x2 = allocate_col_vector(A) solve!(x2,ns2,b)