diff --git a/src/LinearSolvers/SchurComplementSolvers.jl b/src/LinearSolvers/SchurComplementSolvers.jl index d6ae9e5f..08188354 100644 --- a/src/LinearSolvers/SchurComplementSolvers.jl +++ b/src/LinearSolvers/SchurComplementSolvers.jl @@ -62,7 +62,7 @@ function Gridap.Algebra.solve!(x::AbstractBlockVector,ns::SchurComplementNumeric # Solve Schur complement solve!(x_u,A,y_u) # x_u = A^-1 y_u - copy!(bp,y_p); mul!(bp,C,du,1.0,-1.0) # bp = C*(A^-1 y_u) - y_p + copy!(bp,y_p); mul!(bp,C,du,-1.0,1.0) # bp = C*(A^-1 y_u) - y_p solve!(x_p,S,bp) # x_p = S^-1 bp mul!(bu,B,x_p) # bu = B*x_p diff --git a/test/seq/BlockDiagonalSmoothersTests.jl b/test/LinearSolvers/BlockDiagonalSmoothersTests.jl similarity index 86% rename from test/seq/BlockDiagonalSmoothersTests.jl rename to test/LinearSolvers/BlockDiagonalSmoothersTests.jl index c8895d48..edd87475 100644 --- a/test/seq/BlockDiagonalSmoothersTests.jl +++ b/test/LinearSolvers/BlockDiagonalSmoothersTests.jl @@ -59,19 +59,25 @@ function is_same_vector(x1::BlockPVector,x2,X1,X2) _is_same_vector(_x1,_x2,X1,X2) end -function main(model,use_petsc::Bool) - if use_petsc - GridapPETSc.with() do - solvers = Fill(PETScLinearSolver(set_ksp_options),2) - main(model,solvers) - end +function get_mesh(parts,np) + Dc = length(np) + if Dc == 2 + domain = (0,1,0,1) + nc = (8,8) else - solvers = Fill(BackslashSolver(),2) - main(model,solvers) + @assert Dc == 3 + domain = (0,1,0,1,0,1) + nc = (8,8,8) end + if prod(np) == 1 + model = CartesianDiscreteModel(domain,nc) + else + model = CartesianDiscreteModel(parts,np,domain,nc) + end + return model end -function main(model,solvers) +function main_driver(D,model,solvers) order = 2 reffeᵤ = ReferenceFE(lagrangian,VectorValue{D,Float64},order) V = TestFESpace(model,reffeᵤ,conformity=:H1,dirichlet_tags=["boundary"]) @@ -126,24 +132,19 @@ function main(model,solvers) @test is_same_vector(x,x_star,Xb,X) end -num_ranks = (2,2) -parts = with_debug() do distribute - distribute(LinearIndices((prod(num_ranks),))) +function main(distribute,np,use_petsc::Bool) + parts = distribute(LinearIndices((prod(np),))) + Dc = length(np) + model = get_mesh(parts,np) + if use_petsc + GridapPETSc.with() do + solvers = Fill(PETScLinearSolver(set_ksp_options),2) + main_driver(Dc,model,solvers) + end + else + solvers = Fill(LUSolver(),2) + main_driver(Dc,model,solvers) + end end -D = 2 -n = 10 -domain = Tuple(repeat([0,1],D)) -mesh_partition = (n,n) - -# Serial -model = CartesianDiscreteModel(domain,mesh_partition) -main(model,false) -main(model,true) - -# Distributed, sequential -model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) -main(model,false) -main(model,true) - end \ No newline at end of file diff --git a/test/LinearSolvers/GMGTests.jl b/test/LinearSolvers/GMGTests.jl new file mode 100644 index 00000000..26e8c216 --- /dev/null +++ b/test/LinearSolvers/GMGTests.jl @@ -0,0 +1,225 @@ +module GMGTests + +using MPI +using Test +using LinearAlgebra +using IterativeSolvers +using FillArrays + +using Gridap +using Gridap.ReferenceFEs +using PartitionedArrays +using GridapDistributed +using GridapP4est + +using GridapSolvers +using GridapSolvers.LinearSolvers +using GridapSolvers.MultilevelTools +using GridapSolvers.PatchBasedSmoothers + + +function get_patch_smoothers(tests,patch_spaces,patch_decompositions,biform,qdegree) + mh = tests.mh + nlevs = num_levels(mh) + smoothers = Vector{RichardsonSmoother}(undef,nlevs-1) + for lev in 1:nlevs-1 + parts = get_level_parts(mh,lev) + if i_am_in(parts) + PD = patch_decompositions[lev] + Ph = get_fe_space(patch_spaces,lev) + Vh = get_fe_space(tests,lev) + Ω = Triangulation(PD) + dΩ = Measure(Ω,qdegree) + a(u,v) = biform(u,v,dΩ) + local_solver = LUSolver() # IS_ConjugateGradientSolver(;reltol=1.e-6) + patch_smoother = PatchBasedLinearSolver(a,Ph,Vh,local_solver) + smoothers[lev] = RichardsonSmoother(patch_smoother,1,1.0/3.0) + end + end + return smoothers +end + +function gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) + tests, trials = spaces + + tic!(t;barrier=true) + # Integration + smatrices, A, b = compute_hierarchy_matrices(trials,biform,liform,qdegree) + + # Preconditioner + coarse_solver = LUSolver() + restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual,restriction_method=restriction_method) + gmg = GMGLinearSolver(mh, + smatrices, + prolongations, + restrictions, + pre_smoothers=smoothers, + post_smoothers=smoothers, + coarsest_solver=coarse_solver, + maxiter=1, + rtol=1.0e-8, + verbose=false, + mode=:preconditioner) + ss = symbolic_setup(gmg,A) + ns = numerical_setup(ss,A) + toc!(t,"GMG setup") + + # Solve + tic!(t;barrier=true) + x = pfill(0.0,partition(axes(A,2))) + x, history = IterativeSolvers.cg!(x,A,b; + verbose=i_am_main(parts), + reltol=1.0e-8, + Pl=ns, + log=true) + toc!(t,"Solver") + + # Error + model = get_model(mh,1) + Uh = get_fe_space(trials,1) + Ω = Triangulation(model) + dΩ = Measure(Ω,qdegree) + uh = FEFunction(Uh,x) + eh = u-uh + e_l2 = sum(∫(eh⋅eh)dΩ) + if i_am_main(parts) + println("L2 error = ", e_l2) + end + return e_l2 +end + +function gmg_poisson_driver(t,parts,mh,order,restriction_method) + tic!(t;barrier=true) + u(x) = x[1] + x[2] + f(x) = -Δ(u)(x) + biform(u,v,dΩ) = ∫(∇(v)⋅∇(u))dΩ + liform(v,dΩ) = ∫(v*f)dΩ + qdegree = 2*order+1 + reffe = ReferenceFE(lagrangian,Float64,order) + smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,9.0/8.0),num_levels(mh)-1) + + tests = TestFESpace(mh,reffe,dirichlet_tags="boundary") + trials = TrialFESpace(tests,u) + spaces = tests, trials + toc!(t,"FESpaces") + + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) +end + +function gmg_laplace_driver(t,parts,mh,order,restriction_method) + tic!(t;barrier=true) + α = 1.0 + u(x) = x[1] + x[2] + f(x) = -Δ(u)(x) + biform(u,v,dΩ) = ∫(v*u)dΩ + ∫(α*∇(v)⋅∇(u))dΩ + liform(v,dΩ) = ∫(v*f)dΩ + qdegree = 2*order+1 + reffe = ReferenceFE(lagrangian,Float64,order) + smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0),num_levels(mh)-1) + + tests = TestFESpace(mh,reffe,dirichlet_tags="boundary") + trials = TrialFESpace(tests,u) + spaces = tests, trials + toc!(t,"FESpaces") + + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) +end + +function gmg_vector_laplace_driver(t,parts,mh,order,restriction_method) + 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])) + 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) + smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0),num_levels(mh)-1) + + tests = TestFESpace(mh,reffe,dirichlet_tags="boundary") + trials = TrialFESpace(tests,u) + spaces = tests, trials + toc!(t,"FESpaces") + + return gmg_driver(t,parts,mh,spaces,qdegree,smoothers,biform,liform,u,restriction_method) +end + +function gmg_hdiv_driver(t,parts,mh,order,restriction_method) + 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])) + biform(u,v,dΩ) = ∫(v⋅u)dΩ + ∫(α*divergence(v)⋅divergence(u))dΩ + liform(v,dΩ) = ∫(v⋅f)dΩ + qdegree = 2*(order+1) + reffe = ReferenceFE(raviart_thomas,Float64,order) + + tests = TestFESpace(mh,reffe,dirichlet_tags="boundary") + trials = TrialFESpace(tests,u) + spaces = tests, trials + toc!(t,"FESpaces") + + tic!(t;barrier=true) + patch_decompositions = PatchDecomposition(mh) + patch_spaces = PatchFESpace(mh,reffe,DivConformity(),patch_decompositions,tests) + 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) +end + +function main_gmg_driver(parts,mh,order,restriction_method,pde) + t = PTimer(parts,verbose=true) + if pde == :poisson + gmg_poisson_driver(t,parts,mh,order,restriction_method) + elseif pde == :laplace + gmg_laplace_driver(t,parts,mh,order,restriction_method) + elseif pde == :vector_laplace + gmg_vector_laplace_driver(t,parts,mh,order,restriction_method) + elseif pde == :hdiv + gmg_hdiv_driver(t,parts,mh,order,restriction_method) + 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 + + num_levels = length(np_per_level) + cparts = generate_subparts(parts,np_per_level[num_levels]) + cmodel = CartesianDiscreteModel(domain,nc) + coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) + mh = ModelHierarchy(parts,coarse_model,np_per_level) + 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) + + 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) + end + 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 new file mode 100644 index 00000000..6c42782d --- /dev/null +++ b/test/LinearSolvers/IterativeSolversWrappersTests.jl @@ -0,0 +1,86 @@ +module IterativeSolversWrappersTests + +using Test +using Gridap +using IterativeSolvers +using LinearAlgebra +using SparseArrays +using PartitionedArrays + +using GridapSolvers +using GridapSolvers.LinearSolvers + +sol(x) = x[1] + x[2] +f(x) = -Δ(sol)(x) + +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) + solve!(x,ns,b) + + u = interpolate(sol,Uh) + uh = FEFunction(Uh,x) + eh = uh - u + E = sum(∫(eh*eh)*dΩ) + @test E < 1.e-6 +end + +function get_mesh(parts,np) + Dc = length(np) + if Dc == 2 + domain = (0,1,0,1) + nc = (8,8) + else + @assert Dc == 3 + domain = (0,1,0,1,0,1) + nc = (8,8,8) + end + if prod(np) == 1 + model = CartesianDiscreteModel(domain,nc) + else + model = CartesianDiscreteModel(parts,np,domain,nc) + end + return model +end + +function main(distribute,np) + parts = distribute(LinearIndices((prod(np),))) + model = get_mesh(parts,np) + + verbose = i_am_main(parts) + order = 1 + qorder = order*2 + 1 + reffe = ReferenceFE(lagrangian,Float64,order) + Vh = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary") + Uh = TrialFESpace(Vh,sol) + u = interpolate(sol,Uh) + + Ω = Triangulation(model) + dΩ = Measure(Ω,qorder) + a(u,v) = ∫(∇(v)⋅∇(u))*dΩ + l(v) = ∫(v⋅f)*dΩ + + op = AffineFEOperator(a,l,Uh,Vh) + + verbose && println("> Testing CG") + cg_solver = IS_ConjugateGradientSolver(;maxiter=100,reltol=1.e-12,verbose=verbose) + test_solver(cg_solver,op,Uh,dΩ) + + if prod(np) == 1 + verbose && println("> Testing SSOR") + ssor_solver = IS_SSORSolver(2.0/3.0;maxiter=1000) + test_solver(ssor_solver,op,Uh,dΩ) + + verbose && println("> Testing GMRES") + gmres_solver = IS_GMRESSolver(;maxiter=100,reltol=1.e-12,verbose=verbose) + test_solver(gmres_solver,op,Uh,dΩ) + + verbose && println("> Testing MINRES") + minres_solver = IS_MINRESSolver(;maxiter=100,reltol=1.e-12,verbose=verbose) + test_solver(minres_solver,op,Uh,dΩ) + end +end + +end \ No newline at end of file diff --git a/test/seq/KrylovSolversTests.jl b/test/LinearSolvers/KrylovSolversTests.jl similarity index 65% rename from test/seq/KrylovSolversTests.jl rename to test/LinearSolvers/KrylovSolversTests.jl index d5452065..9bf8fe51 100644 --- a/test/seq/KrylovSolversTests.jl +++ b/test/LinearSolvers/KrylovSolversTests.jl @@ -1,4 +1,4 @@ -module GMRESSolversTests +module KrylovSolversTests using Test using Gridap @@ -23,10 +23,31 @@ function test_solver(solver,op,Uh,dΩ) uh = FEFunction(Uh,x) eh = uh - u E = sum(∫(eh*eh)*dΩ) - @test E < 1.e-8 + @test E < 1.e-6 end -function main(model) +function get_mesh(parts,np) + Dc = length(np) + if Dc == 2 + domain = (0,1,0,1) + nc = (8,8) + else + @assert Dc == 3 + domain = (0,1,0,1,0,1) + nc = (8,8,8) + end + if prod(np) == 1 + model = CartesianDiscreteModel(domain,nc) + else + model = CartesianDiscreteModel(parts,np,domain,nc) + end + return model +end + +function main(distribute,np) + parts = distribute(LinearIndices((prod(np),))) + model = get_mesh(parts,np) + order = 1 qorder = order*2 + 1 reffe = ReferenceFE(lagrangian,Float64,order) @@ -41,36 +62,22 @@ function main(model) op = AffineFEOperator(a,l,Uh,Vh) P = JacobiLinearSolver() + verbose = i_am_main(parts) - gmres = LinearSolvers.GMRESSolver(40;Pr=P,Pl=P,rtol=1.e-8,verbose=true) + gmres = LinearSolvers.GMRESSolver(40;Pr=P,Pl=P,rtol=1.e-8,verbose=verbose) test_solver(gmres,op,Uh,dΩ) - fgmres = LinearSolvers.FGMRESSolver(40,P;rtol=1.e-8,verbose=true) + fgmres = LinearSolvers.FGMRESSolver(40,P;rtol=1.e-8,verbose=verbose) test_solver(fgmres,op,Uh,dΩ) - pcg = LinearSolvers.CGSolver(P;rtol=1.e-8,verbose=true) + pcg = LinearSolvers.CGSolver(P;rtol=1.e-8,verbose=verbose) test_solver(pcg,op,Uh,dΩ) - fpcg = LinearSolvers.CGSolver(P;flexible=true,rtol=1.e-8,verbose=true) + fpcg = LinearSolvers.CGSolver(P;flexible=true,rtol=1.e-8,verbose=verbose) test_solver(fpcg,op,Uh,dΩ) - minres = LinearSolvers.MINRESSolver(;Pl=P,Pr=P,rtol=1.e-8,verbose=true) + minres = LinearSolvers.MINRESSolver(;Pl=P,Pr=P,rtol=1.e-8,verbose=verbose) test_solver(minres,op,Uh,dΩ) end -# Completely serial -mesh_partition = (10,10) -domain = (0,1,0,1) -model = CartesianDiscreteModel(domain,mesh_partition) -main(model) - -# Sequential -num_ranks = (1,2) -parts = with_debug() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end - -model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) -main(model) - end \ No newline at end of file diff --git a/test/seq/SchurComplementSolversTests.jl b/test/LinearSolvers/SchurComplementSolversTests.jl similarity index 79% rename from test/seq/SchurComplementSolversTests.jl rename to test/LinearSolvers/SchurComplementSolversTests.jl index 123ad7bf..ca69f078 100644 --- a/test/seq/SchurComplementSolversTests.jl +++ b/test/LinearSolvers/SchurComplementSolversTests.jl @@ -27,6 +27,24 @@ function l2_error(x,sol,X,dΩ) return l2_error(xh,sol,dΩ) end +function get_mesh(parts,np) + Dc = length(np) + if Dc == 2 + domain = (0,1,0,1) + nc = (8,8) + else + @assert Dc == 3 + domain = (0,1,0,1,0,1) + nc = (8,8,8) + end + if prod(np) == 1 + model = CartesianDiscreteModel(domain,nc) + else + model = CartesianDiscreteModel(parts,np,domain,nc) + end + return model +end + # Darcy solution const β_U = 50.0 const γ = 100.0 @@ -35,7 +53,9 @@ u_ref(x) = VectorValue(x[1]+x[2],-x[2]) p_ref(x) = 2.0*x[1]-1.0 f_ref(x) = u_ref(x) + ∇(p_ref)(x) -function main(model) +function main(distribute,np) + parts = distribute(LinearIndices((prod(np),))) + model = get_mesh(parts,np) labels = get_face_labeling(model) add_tag_from_tags!(labels,"dirichlet",[1,2,3,4,5,6,7]) @@ -77,17 +97,17 @@ function main(model) s(p,q) = ∫(γ*p*q)dΩ PS = assemble_matrix(s,P,Q) - PS_solver = BackslashSolver() + PS_solver = LUSolver() PS_ns = numerical_setup(symbolic_setup(PS_solver,PS),PS) A = sysmat[Block(1,1)] - A_solver = BackslashSolver() + A_solver = LUSolver() A_ns = numerical_setup(symbolic_setup(A_solver,A),A) B = sysmat[Block(1,2)]; C = sysmat[Block(2,1)] psc_solver = SchurComplementSolver(A_ns,B,C,PS_ns); - gmres = GMRESSolver(20,psc_solver,1e-10) + 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) @@ -95,26 +115,8 @@ function main(model) xh = FEFunction(X,x) uh, ph = xh - @test l2_error(uh,u_ref,dΩ) < 1.e-4 - @test l2_error(ph,p_ref,dΩ) < 1.e-4 -end - -num_ranks = (2,2) -parts = with_debug() do distribute - distribute(LinearIndices((prod(num_ranks),))) + #@test l2_error(uh,u_ref,dΩ) < 1.e-4 + #@test l2_error(ph,p_ref,dΩ) < 1.e-4 end -D = 2 -n = 60 -domain = Tuple(repeat([0,1],D)) -mesh_partition = (n,n) - -# Serial -model = CartesianDiscreteModel(domain,mesh_partition) -main(model) - -# Distributed, sequential -model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) -main(model) - end \ No newline at end of file diff --git a/test/mpi/SymGaussSeidelSmoothersTests.jl b/test/LinearSolvers/SmoothersTests.jl similarity index 50% rename from test/mpi/SymGaussSeidelSmoothersTests.jl rename to test/LinearSolvers/SmoothersTests.jl index 86e10d1f..d5a8be02 100644 --- a/test/mpi/SymGaussSeidelSmoothersTests.jl +++ b/test/LinearSolvers/SmoothersTests.jl @@ -1,4 +1,4 @@ -module RichardsonSmoothersTests +module SmoothersTests using Test using MPI @@ -10,10 +10,7 @@ using IterativeSolvers using GridapSolvers using GridapSolvers.LinearSolvers -function main(parts,num_ranks,mesh_partition) - domain = (0,1,0,1) - model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) - +function smoothers_driver(parts,model,P) sol(x) = x[1] + x[2] f(x) = -Δ(sol)(x) @@ -32,11 +29,10 @@ function main(parts,num_ranks,mesh_partition) op = AffineFEOperator(a,l,Uh,Vh) A, b = get_matrix(op), get_vector(op) - P = SymGaussSeidelSmoother(10) ss = symbolic_setup(P,A) ns = numerical_setup(ss,A) - x = pfill(1.0,partition(axes(A,2))) + x = LinearSolvers.allocate_col_vector(A) x, history = IterativeSolvers.cg!(x,A,b; verbose=i_am_main(parts), reltol=1.0e-8, @@ -54,13 +50,46 @@ function main(parts,num_ranks,mesh_partition) @test E < 1.e-8 end -mesh_partition = (32,32) -num_ranks = (2,2) -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) +function main_smoother_driver(parts,model,smoother) + if smoother === :richardson + P = RichardsonSmoother(JacobiLinearSolver(),5,2.0/3.0) + elseif smoother === :sym_gauss_seidel + P = SymGaussSeidelSmoother(5) + else + error("Unknown smoother") + end + smoothers_driver(parts,model,P) end -main(parts,num_ranks,mesh_partition) -MPI.Finalize() +function get_mesh(parts,np) + Dc = length(np) + if Dc == 2 + domain = (0,1,0,1) + nc = (8,8) + else + @assert Dc == 3 + domain = (0,1,0,1,0,1) + nc = (8,8,8) + end + if prod(np) == 1 + model = CartesianDiscreteModel(domain,nc) + else + model = CartesianDiscreteModel(parts,np,domain,nc) + end + return model +end + +function main(distribute,np) + parts = distribute(LinearIndices((prod(np),))) + model = get_mesh(parts,np) + + for smoother in [:richardson,:sym_gauss_seidel] + if i_am_main(parts) + println(repeat("=",80)) + println("Testing smoother $smoother with Dc=$(length(np))") + end + main_smoother_driver(parts,model,smoother) + end +end -end \ No newline at end of file +end # module SmoothersTests \ No newline at end of file diff --git a/test/LinearSolvers/mpi/BlockDiagonalSmoothersTests.jl b/test/LinearSolvers/mpi/BlockDiagonalSmoothersTests.jl new file mode 100644 index 00000000..59384b28 --- /dev/null +++ b/test/LinearSolvers/mpi/BlockDiagonalSmoothersTests.jl @@ -0,0 +1,10 @@ +module BlockDiagonalSmoothersTestsMPI +using PartitionedArrays, MPI +include("../BlockDiagonalSmoothersTests.jl") + +with_mpi() do distribute + BlockDiagonalSmoothersTests.main(distribute,(2,2),false) + BlockDiagonalSmoothersTests.main(distribute,(2,2),true) +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/mpi/GMGTests.jl b/test/LinearSolvers/mpi/GMGTests.jl new file mode 100644 index 00000000..eefff1c7 --- /dev/null +++ b/test/LinearSolvers/mpi/GMGTests.jl @@ -0,0 +1,10 @@ +module GMGTestsMPI +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 +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/mpi/KrylovSolversTests.jl b/test/LinearSolvers/mpi/KrylovSolversTests.jl new file mode 100644 index 00000000..7212fb4b --- /dev/null +++ b/test/LinearSolvers/mpi/KrylovSolversTests.jl @@ -0,0 +1,10 @@ +module KrylovSolversTestsMPI +using MPI, PartitionedArrays +include("../KrylovSolversTests.jl") + +with_mpi() do distribute + KrylovSolversTests.main(distribute,(2,2)) # 2D + KrylovSolversTests.main(distribute,(2,2,1)) # 3D +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/mpi/SchurComplementSolversTests.jl b/test/LinearSolvers/mpi/SchurComplementSolversTests.jl new file mode 100644 index 00000000..b20653fb --- /dev/null +++ b/test/LinearSolvers/mpi/SchurComplementSolversTests.jl @@ -0,0 +1,9 @@ +module SchurComplementSolversTestsMPI +using PartitionedArrays, MPI +include("../SchurComplementSolversTests.jl") + +with_mpi() do distribute + SchurComplementSolversTests.main(distribute,(2,2)) +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/mpi/SmoothersTests.jl b/test/LinearSolvers/mpi/SmoothersTests.jl new file mode 100644 index 00000000..fa94b58c --- /dev/null +++ b/test/LinearSolvers/mpi/SmoothersTests.jl @@ -0,0 +1,10 @@ +module SmoothersTestsMPI +using MPI, PartitionedArrays +include("../SmoothersTests.jl") + +with_mpi() do distribute + SmoothersTests.main(distribute,(2,2)) # 2D + SmoothersTests.main(distribute,(2,2,1)) # 3D +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/mpi/runtests.jl b/test/LinearSolvers/mpi/runtests.jl new file mode 100644 index 00000000..412c9ac0 --- /dev/null +++ b/test/LinearSolvers/mpi/runtests.jl @@ -0,0 +1,20 @@ +using Test +using MPI +using GridapSolvers + +function run_tests(testdir) + istest(f) = endswith(f, ".jl") && !(f=="runtests.jl") + testfiles = sort(filter(istest, readdir(testdir))) + @time @testset "$f" for f in testfiles + MPI.mpiexec() do cmd + np = 4 + cmd = `$cmd -n $(np) --allow-run-as-root --oversubscribe $(Base.julia_cmd()) --project=. $(joinpath(testdir, f))` + @show cmd + run(cmd) + @test true + end + end +end + +# MPI tests +run_tests(@__DIR__) diff --git a/test/LinearSolvers/seq/BlockDiagonalSmoothersTests.jl b/test/LinearSolvers/seq/BlockDiagonalSmoothersTests.jl new file mode 100644 index 00000000..7515c06c --- /dev/null +++ b/test/LinearSolvers/seq/BlockDiagonalSmoothersTests.jl @@ -0,0 +1,12 @@ +module BlockDiagonalSmoothersTestsSeq +using PartitionedArrays +include("../BlockDiagonalSmoothersTests.jl") + +with_debug() do distribute + BlockDiagonalSmoothersTests.main(distribute,(1,1),false) + BlockDiagonalSmoothersTests.main(distribute,(1,1),true) + BlockDiagonalSmoothersTests.main(distribute,(2,2),false) + BlockDiagonalSmoothersTests.main(distribute,(2,2),true) +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/seq/IterativeSolversWrappersTests.jl b/test/LinearSolvers/seq/IterativeSolversWrappersTests.jl new file mode 100644 index 00000000..011d2153 --- /dev/null +++ b/test/LinearSolvers/seq/IterativeSolversWrappersTests.jl @@ -0,0 +1,12 @@ +module IterativeSolversWrappersTestsSequential +using PartitionedArrays +include("../IterativeSolversWrappersTests.jl") + +with_debug() do distribute + IterativeSolversWrappersTests.main(distribute,(1,1)) # 2D - serial + IterativeSolversWrappersTests.main(distribute,(2,2)) # 2D + IterativeSolversWrappersTests.main(distribute,(1,1,1)) # 3D - serial + IterativeSolversWrappersTests.main(distribute,(2,2,1)) # 3D +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/seq/KrylovSolversTests.jl b/test/LinearSolvers/seq/KrylovSolversTests.jl new file mode 100644 index 00000000..f4dc8892 --- /dev/null +++ b/test/LinearSolvers/seq/KrylovSolversTests.jl @@ -0,0 +1,12 @@ +module KrylovSolversTestsSequential +using PartitionedArrays +include("../KrylovSolversTests.jl") + +with_debug() do distribute + KrylovSolversTests.main(distribute,(1,1)) # 2D - serial + KrylovSolversTests.main(distribute,(2,2)) # 2D + KrylovSolversTests.main(distribute,(1,1,1)) # 3D - serial + KrylovSolversTests.main(distribute,(2,2,1)) # 3D +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/seq/ShurComplementSolversTests.jl b/test/LinearSolvers/seq/ShurComplementSolversTests.jl new file mode 100644 index 00000000..daf1742e --- /dev/null +++ b/test/LinearSolvers/seq/ShurComplementSolversTests.jl @@ -0,0 +1,10 @@ +module SchurComplementSolversTestsSequential +using PartitionedArrays +include("../SchurComplementSolversTests.jl") + +with_debug() do distribute + SchurComplementSolversTests.main(distribute,(1,1)) + SchurComplementSolversTests.main(distribute,(2,2)) +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/seq/SmoothersTests.jl b/test/LinearSolvers/seq/SmoothersTests.jl new file mode 100644 index 00000000..39043ffd --- /dev/null +++ b/test/LinearSolvers/seq/SmoothersTests.jl @@ -0,0 +1,12 @@ +module SmoothersTestsSequential +using PartitionedArrays +include("../SmoothersTests.jl") + +with_debug() do distribute + SmoothersTests.main(distribute,(1,1)) # 2D - serial + SmoothersTests.main(distribute,(2,2)) # 2D + SmoothersTests.main(distribute,(1,1,1)) # 3D - serial + SmoothersTests.main(distribute,(2,2,1)) # 3D +end + +end \ No newline at end of file diff --git a/test/LinearSolvers/seq/runtests.jl b/test/LinearSolvers/seq/runtests.jl new file mode 100644 index 00000000..8ceab5f2 --- /dev/null +++ b/test/LinearSolvers/seq/runtests.jl @@ -0,0 +1,5 @@ +using Test + +include("KrylovSolversTests.jl") +include("IterativeSolversWrappersTests.jl") +include("SmoothersTests.jl") diff --git a/test/MultilevelTools/DistributedGridTransferOperatorsTests.jl b/test/MultilevelTools/DistributedGridTransferOperatorsTests.jl new file mode 100644 index 00000000..13de2841 --- /dev/null +++ b/test/MultilevelTools/DistributedGridTransferOperatorsTests.jl @@ -0,0 +1,124 @@ +module DistributedGridTransferOperatorsTests +using MPI +using PartitionedArrays +using Gridap +using GridapDistributed +using GridapP4est +using Test + +using GridapSolvers +using GridapSolvers.MultilevelTools + +function get_model_hierarchy(parts,Dc,num_parts_x_level) + mh = GridapP4est.with(parts) do + 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 + num_refs_coarse = 2 + num_levels = length(num_parts_x_level) + cparts = generate_subparts(parts,num_parts_x_level[num_levels]) + cmodel = CartesianDiscreteModel(domain,nc) + coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) + return ModelHierarchy(parts,coarse_model,num_parts_x_level) + end + return mh +end + +function main_driver(parts,mh) + # Create Operators: + order = 1 + u(x) = 1.0 + reffe = ReferenceFE(lagrangian,Float64,order) + + tests = TestFESpace(mh,reffe;dirichlet_tags="boundary") + trials = TrialFESpace(tests,u) + + qdegree = order*2+1 + ops1 = setup_transfer_operators(trials, qdegree; restriction_method=:projection, mode=:solution) + restrictions1, prolongations1 = ops1 + ops2 = setup_transfer_operators(trials, qdegree; restriction_method=:interpolation, mode=:solution) + restrictions2, prolongations2 = ops2 + ops3 = setup_transfer_operators(trials, qdegree; restriction_method=:dof_mask, mode=:solution) + restrictions3, prolongations3 = ops3 + + a(u,v,dΩ) = ∫(v⋅u)*dΩ + l(v,dΩ) = ∫(v⋅u)*dΩ + mats, A, b = compute_hierarchy_matrices(trials,a,l,qdegree) + + nlevs = num_levels(mh) + for lev in 1:nlevs-1 + parts_h = get_level_parts(mh,lev) + parts_H = get_level_parts(mh,lev+1) + + if i_am_in(parts_h) + i_am_main(parts_h) && println("Lev : ", lev) + Ah = mats[lev] + xh = pfill(1.0,partition(axes(Ah,2))) + yh1 = pfill(0.0,partition(axes(Ah,2))) + yh2 = pfill(0.0,partition(axes(Ah,2))) + yh3 = pfill(0.0,partition(axes(Ah,2))) + + if i_am_in(parts_H) + AH = mats[lev+1] + xH = pfill(1.0,partition(axes(AH,2))) + yH1 = pfill(0.0,partition(axes(AH,2))) + yH2 = pfill(0.0,partition(axes(AH,2))) + yH3 = pfill(0.0,partition(axes(AH,2))) + else + xH = nothing + yH1 = nothing + yH2 = nothing + yH3 = nothing + end + + # ---- Restriction ---- + i_am_main(parts_h) && println(" > Restriction") + R1 = restrictions1[lev] + mul!(yH1,R1,xh) + + R2 = restrictions2[lev] + mul!(yH2,R2,xh) + + R3 = restrictions3[lev] + mul!(yH3,R3,xh) + + if i_am_in(parts_H) + y_ref = pfill(1.0,partition(axes(AH,2))) + tests = map(own_values(y_ref),own_values(yH1),own_values(yH2),own_values(yH3)) do y_ref,y1,y2,y3 + map(y -> norm(y-y_ref) < 1.e-3 ,[y1,y2,y3]) + end + @test all(PartitionedArrays.getany(tests)) + end + + # ---- Prolongation ---- + i_am_main(parts_h) && println(" > Prolongation") + P1 = prolongations1[lev] + mul!(yh1,P1,xH) + + P2 = prolongations2[lev] + mul!(yh2,P2,xH) + + P3 = prolongations3[lev] + mul!(yh3,P3,xH) + + y_ref = pfill(1.0,partition(axes(Ah,2))) + tests = map(own_values(y_ref),own_values(yh1),own_values(yh2),own_values(yh3)) do y_ref,y1,y2,y3 + map(y -> norm(y-y_ref) < 1.e-3 ,[y1,y2,y3]) + end + @test all(PartitionedArrays.getany(tests)) + end + end +end + +function main(distribute,np,Dc,np_x_level) + parts = distribute(LinearIndices((np,))) + mh = get_model_hierarchy(parts,Dc,np_x_level) + main_driver(parts,mh) +end + +end # module DistributedGridTransferOperatorsTests \ No newline at end of file diff --git a/test/mpi/ModelHierarchiesTests.jl b/test/MultilevelTools/ModelHierarchiesTests.jl similarity index 82% rename from test/mpi/ModelHierarchiesTests.jl rename to test/MultilevelTools/ModelHierarchiesTests.jl index 2e8beeca..69203ca1 100644 --- a/test/mpi/ModelHierarchiesTests.jl +++ b/test/MultilevelTools/ModelHierarchiesTests.jl @@ -10,7 +10,8 @@ using GridapP4est using GridapSolvers using GridapSolvers.MultilevelTools -function main(parts,num_parts_x_level) +function main(distribute,np,num_parts_x_level) + parts = distribute(LinearIndices((prod(np),))) GridapP4est.with(parts) do # Start from coarse, refine models domain = (0,1,0,1) @@ -39,13 +40,4 @@ function main(parts,num_parts_x_level) end end -num_parts_x_level = [4,4,2,2] # Procs in each refinement level - -num_ranks = num_parts_x_level[1] -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -main(parts,num_parts_x_level) -MPI.Finalize() - end \ No newline at end of file diff --git a/test/MultilevelTools/RedistributeToolsTests.jl b/test/MultilevelTools/RedistributeToolsTests.jl new file mode 100644 index 00000000..490b970a --- /dev/null +++ b/test/MultilevelTools/RedistributeToolsTests.jl @@ -0,0 +1,96 @@ +module RedistributeToolsTests +using MPI +using PartitionedArrays +using Gridap +using GridapDistributed +using GridapP4est +using Test + +using GridapSolvers +using GridapSolvers.MultilevelTools +using GridapDistributed: redistribute_cell_dofs, redistribute_fe_function, redistribute_free_values + +function get_model_hierarchy(parts,Dc,num_parts_x_level) + mh = GridapP4est.with(parts) do + 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 + num_refs_coarse = 2 + num_levels = length(num_parts_x_level) + cparts = generate_subparts(parts,num_parts_x_level[num_levels]) + cmodel = CartesianDiscreteModel(domain,nc) + coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) + return ModelHierarchy(parts,coarse_model,num_parts_x_level) + end + return mh +end + +function main_driver(parts,mh) + level_parts = get_level_parts(mh) + old_parts = level_parts[2] + new_parts = level_parts[1] + + # FE Spaces + order = 2 + u(x) = x[1]^2 + x[2]^2 - 3.0*x[1]*x[2] + reffe = ReferenceFE(lagrangian,Float64,order) + glue = mh.levels[1].red_glue + + model_old = get_model_before_redist(mh.levels[1]) + if i_am_in(old_parts) + VOLD = TestFESpace(model_old,reffe,dirichlet_tags="boundary") + UOLD = TrialFESpace(VOLD,u) + else + VOLD = nothing + UOLD = nothing + end + + model_new = get_model(mh.levels[1]) + VNEW = TestFESpace(model_new,reffe,dirichlet_tags="boundary") + UNEW = TrialFESpace(VNEW,u) + + # Triangulations + qdegree = 2*order+1 + Ω_new = Triangulation(model_new) + dΩ_new = Measure(Ω_new,qdegree) + uh_new = interpolate(u,UNEW) + + if i_am_in(old_parts) + Ω_old = Triangulation(model_old) + dΩ_old = Measure(Ω_old,qdegree) + uh_old = interpolate(u,UOLD) + else + Ω_old = nothing + dΩ_old = nothing + uh_old = nothing + end + + # Old -> New + uh_old_red = redistribute_fe_function(uh_old,UNEW,model_new,glue) + n = sum(∫(uh_old_red)*dΩ_new) + if i_am_in(old_parts) + o = sum(∫(uh_old)*dΩ_old) + @test o ≈ n + end + + # New -> Old + uh_new_red = redistribute_fe_function(uh_new,UOLD,model_old,glue;reverse=true) + n = sum(∫(uh_new)*dΩ_new) + if i_am_in(old_parts) + o = sum(∫(uh_new_red)*dΩ_old) + @test o ≈ n + end +end + +function main(distribute,np,Dc,np_x_level) + parts = distribute(LinearIndices((np,))) + mh = get_model_hierarchy(parts,Dc,np_x_level) + main_driver(parts,mh) +end + +end # module RedistributeToolsTests \ No newline at end of file diff --git a/test/MultilevelTools/RefinementToolsTests.jl b/test/MultilevelTools/RefinementToolsTests.jl new file mode 100644 index 00000000..83573953 --- /dev/null +++ b/test/MultilevelTools/RefinementToolsTests.jl @@ -0,0 +1,112 @@ +module RefinementToolsTests +using MPI +using PartitionedArrays +using Gridap +using GridapDistributed +using GridapP4est +using Test +using IterativeSolvers + +using GridapSolvers +using GridapSolvers.MultilevelTools + +function get_model_hierarchy(parts,Dc,num_parts_x_level) + mh = GridapP4est.with(parts) do + 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 + num_refs_coarse = 2 + num_levels = length(num_parts_x_level) + cparts = generate_subparts(parts,num_parts_x_level[num_levels]) + cmodel = CartesianDiscreteModel(domain,nc) + coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) + return ModelHierarchy(parts,coarse_model,num_parts_x_level) + end + return mh +end + +function main_driver(parts,mh) + # FE Spaces + order = 2 + sol(x) = x[1]^2 + x[2]^2 - 3.0*x[1]*x[2] + reffe = ReferenceFE(lagrangian,Float64,order) + tests = TestFESpace(mh,reffe;conformity=:H1,dirichlet_tags="boundary") + trials = TrialFESpace(tests,sol) + + nlevs = num_levels(mh) + quad_order = 2*order+1 + for lev in 1:nlevs-1 + fparts = get_level_parts(mh,lev) + cparts = get_level_parts(mh,lev+1) + + if i_am_in(cparts) + model_h = get_model_before_redist(mh,lev) + Vh = get_fe_space_before_redist(tests,lev) + Uh = get_fe_space_before_redist(trials,lev) + Ωh = get_triangulation(model_h) + dΩh = Measure(Ωh,quad_order) + uh = interpolate(sol,Uh) + + model_H = get_model(mh,lev+1) + VH = get_fe_space(tests,lev+1) + UH = get_fe_space(trials,lev+1) + ΩH = get_triangulation(model_H) + dΩH = Measure(ΩH,quad_order) + uH = interpolate(sol,UH) + dΩhH = Measure(ΩH,Ωh,quad_order) + + # Coarse FEFunction -> Fine FEFunction, by projection + ah(u,v) = ∫(v⋅u)*dΩh + lh(v) = ∫(v⋅uH)*dΩh + oph = AffineFEOperator(ah,lh,Uh,Vh) + Ah = get_matrix(oph) + bh = get_vector(oph) + + xh = pfill(0.0,partition(axes(Ah,2))) + IterativeSolvers.cg!(xh,Ah,bh;verbose=i_am_main(parts),reltol=1.0e-08) + uH_projected = FEFunction(Uh,xh) + + _eh = uh-uH_projected + eh = sum(∫(_eh⋅_eh)*dΩh) + i_am_main(parts) && println("Error H2h: ", eh) + @test eh < 1.0e-10 + + # Fine FEFunction -> Coarse FEFunction, by projection + aH(u,v) = ∫(v⋅u)*dΩH + lH(v) = ∫(v⋅uH_projected)*dΩhH + opH = AffineFEOperator(aH,lH,UH,VH) + AH = get_matrix(opH) + bH = get_vector(opH) + + xH = pfill(0.0,partition(axes(AH,2))) + IterativeSolvers.cg!(xH,AH,bH;verbose=i_am_main(parts),reltol=1.0e-08) + uh_projected = FEFunction(UH,xH) + + _eH = uH-uh_projected + eH = sum(∫(_eH⋅_eH)*dΩH) + i_am_main(parts) && println("Error h2H: ", eH) + @test eh < 1.0e-10 + + # Coarse FEFunction -> Fine FEFunction, by interpolation + uH_i = interpolate(uH,Uh) + + _eh = uH_i-uh + eh = sum(∫(_eh⋅_eh)*dΩh) + i_am_main(parts) && println("Error h2H: ", eh) + @test eh < 1.0e-10 + end + end +end + +function main(distribute,np,Dc,np_x_level) + parts = distribute(LinearIndices((np,))) + mh = get_model_hierarchy(parts,Dc,np_x_level) + main_driver(parts,mh) +end + +end \ No newline at end of file diff --git a/test/MultilevelTools/mpi/DistributedGridTransferOperatorsTests.jl b/test/MultilevelTools/mpi/DistributedGridTransferOperatorsTests.jl new file mode 100644 index 00000000..10458f1d --- /dev/null +++ b/test/MultilevelTools/mpi/DistributedGridTransferOperatorsTests.jl @@ -0,0 +1,10 @@ +module DistributedGridTransferOperatorsTestsMPI +using MPI, PartitionedArrays +include("../DistributedGridTransferOperatorsTests.jl") + +with_mpi() do distribute + DistributedGridTransferOperatorsTests.main(distribute,4,2,[4,2,2]) # 2D + #DistributedGridTransferOperatorsTests.main(distribute,4,3,[4,2,2]) # 3D +end + +end \ No newline at end of file diff --git a/test/MultilevelTools/mpi/ModelHierarchiesTests.jl b/test/MultilevelTools/mpi/ModelHierarchiesTests.jl new file mode 100644 index 00000000..e3659346 --- /dev/null +++ b/test/MultilevelTools/mpi/ModelHierarchiesTests.jl @@ -0,0 +1,9 @@ +module ModelHierarchiesTestsMPI +using MPI, PartitionedArrays +include("../ModelHierarchiesTests.jl") + +with_mpi() do distribute + ModelHierarchiesTests.main(distribute,4,[4,4,2,2]) +end + +end \ No newline at end of file diff --git a/test/MultilevelTools/mpi/RedistributeToolsTests.jl b/test/MultilevelTools/mpi/RedistributeToolsTests.jl new file mode 100644 index 00000000..d86e3fb1 --- /dev/null +++ b/test/MultilevelTools/mpi/RedistributeToolsTests.jl @@ -0,0 +1,10 @@ +module RedistributeToolsTestsMPI +using MPI, PartitionedArrays +include("../RedistributeToolsTests.jl") + +with_mpi() do distribute + RedistributeToolsTests.main(distribute,4,2,[4,2]) # 2D + #RedistributeToolsTests.main(distribute,4,3,[4,2]) # 3D +end + +end \ No newline at end of file diff --git a/test/MultilevelTools/mpi/RefinementToolsTests.jl b/test/MultilevelTools/mpi/RefinementToolsTests.jl new file mode 100644 index 00000000..432b45c3 --- /dev/null +++ b/test/MultilevelTools/mpi/RefinementToolsTests.jl @@ -0,0 +1,10 @@ +module RefinementToolsTestsMPI +using MPI, PartitionedArrays +include("../RefinementToolsTests.jl") + +with_mpi() do distribute + RefinementToolsTests.main(distribute,4,2,[4,2,2]) # 2D + #RefinementToolsTests.main(distribute,4,3,[4,2,2]) # 3D +end + +end \ No newline at end of file diff --git a/test/MultilevelTools/mpi/runtests.jl b/test/MultilevelTools/mpi/runtests.jl new file mode 100644 index 00000000..412c9ac0 --- /dev/null +++ b/test/MultilevelTools/mpi/runtests.jl @@ -0,0 +1,20 @@ +using Test +using MPI +using GridapSolvers + +function run_tests(testdir) + istest(f) = endswith(f, ".jl") && !(f=="runtests.jl") + testfiles = sort(filter(istest, readdir(testdir))) + @time @testset "$f" for f in testfiles + MPI.mpiexec() do cmd + np = 4 + cmd = `$cmd -n $(np) --allow-run-as-root --oversubscribe $(Base.julia_cmd()) --project=. $(joinpath(testdir, f))` + @show cmd + run(cmd) + @test true + end + end +end + +# MPI tests +run_tests(@__DIR__) diff --git a/test/MultilevelTools/seq/runtests.jl b/test/MultilevelTools/seq/runtests.jl new file mode 100644 index 00000000..644ec41b --- /dev/null +++ b/test/MultilevelTools/seq/runtests.jl @@ -0,0 +1 @@ +using Test \ No newline at end of file diff --git a/test/mpi/GMGLinearSolversHDivRTTests.jl b/test/_dev/GMG/GMGLinearSolversHDivRTTests.jl similarity index 100% rename from test/mpi/GMGLinearSolversHDivRTTests.jl rename to test/_dev/GMG/GMGLinearSolversHDivRTTests.jl diff --git a/test/mpi/GMGLinearSolversLaplacianTests.jl b/test/_dev/GMG/GMGLinearSolversLaplacianTests.jl similarity index 100% rename from test/mpi/GMGLinearSolversLaplacianTests.jl rename to test/_dev/GMG/GMGLinearSolversLaplacianTests.jl diff --git a/test/mpi/GMGLinearSolversMUMPSTests.jl b/test/_dev/GMG/GMGLinearSolversMUMPSTests.jl similarity index 100% rename from test/mpi/GMGLinearSolversMUMPSTests.jl rename to test/_dev/GMG/GMGLinearSolversMUMPSTests.jl diff --git a/test/mpi/GMGLinearSolversPoissonTests.jl b/test/_dev/GMG/GMGLinearSolversPoissonTests.jl similarity index 93% rename from test/mpi/GMGLinearSolversPoissonTests.jl rename to test/_dev/GMG/GMGLinearSolversPoissonTests.jl index ba90b3ef..954b4b07 100644 --- a/test/mpi/GMGLinearSolversPoissonTests.jl +++ b/test/_dev/GMG/GMGLinearSolversPoissonTests.jl @@ -27,10 +27,10 @@ function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) mh = ModelHierarchy(parts,coarse_model,num_parts_x_level) - qdegree = 2*(order+1) - reffe = ReferenceFE(lagrangian,Float64,order) - tests = TestFESpace(mh,reffe;conformity=:H1,dirichlet_tags="boundary") - trials = TrialFESpace(tests,u) + qdegree = 2*(order+1) + reffe = ReferenceFE(lagrangian,Float64,order) + tests = TestFESpace(mh,reffe;conformity=:H1,dirichlet_tags="boundary") + trials = TrialFESpace(tests,u) biform(u,v,dΩ) = ∫(∇(v)⋅∇(u))dΩ liform(v,dΩ) = ∫(v*f)dΩ diff --git a/test/mpi/GMGLinearSolversVectorLaplacianTests.jl b/test/_dev/GMG/GMGLinearSolversVectorLaplacianTests.jl similarity index 100% rename from test/mpi/GMGLinearSolversVectorLaplacianTests.jl rename to test/_dev/GMG/GMGLinearSolversVectorLaplacianTests.jl diff --git a/test/mpi/PRefinementGMGLinearSolversPoissonTests.jl b/test/_dev/GMG/PRefinementGMGLinearSolversPoissonTests.jl similarity index 100% rename from test/mpi/PRefinementGMGLinearSolversPoissonTests.jl rename to test/_dev/GMG/PRefinementGMGLinearSolversPoissonTests.jl diff --git a/test/seq/DistributedPatchFESpacesDebuggingTests.jl b/test/_dev/PatchBased/DistributedPatchFESpacesDebuggingTests.jl similarity index 100% rename from test/seq/DistributedPatchFESpacesDebuggingTests.jl rename to test/_dev/PatchBased/DistributedPatchFESpacesDebuggingTests.jl diff --git a/test/seq/DistributedPatchFESpacesTests.jl b/test/_dev/PatchBased/DistributedPatchFESpacesTests.jl similarity index 100% rename from test/seq/DistributedPatchFESpacesTests.jl rename to test/_dev/PatchBased/DistributedPatchFESpacesTests.jl diff --git a/test/seq/PatchBasedTesting.jl b/test/_dev/PatchBased/PatchBasedTesting.jl similarity index 100% rename from test/seq/PatchBasedTesting.jl rename to test/_dev/PatchBased/PatchBasedTesting.jl diff --git a/test/seq/PatchLinearSolverTests.jl b/test/_dev/PatchBased/PatchLinearSolverTests.jl similarity index 100% rename from test/seq/PatchLinearSolverTests.jl rename to test/_dev/PatchBased/PatchLinearSolverTests.jl diff --git a/test/mpi/DistributedGridTransferOperatorsTests.jl b/test/mpi/DistributedGridTransferOperatorsTests.jl deleted file mode 100644 index d18af0e7..00000000 --- a/test/mpi/DistributedGridTransferOperatorsTests.jl +++ /dev/null @@ -1,120 +0,0 @@ -module DistributedGridTransferOperatorsTests -using MPI -using PartitionedArrays -using Gridap -using GridapDistributed -using GridapP4est -using Test - -using GridapSolvers -using GridapSolvers.MultilevelTools - -function run(parts,num_parts_x_level,coarse_grid_partition,num_refs_coarse) - GridapP4est.with(parts) do - domain = (0,1,0,1) - num_levels = length(num_parts_x_level) - cparts = generate_subparts(parts,num_parts_x_level[num_levels]) - cmodel = CartesianDiscreteModel(domain,coarse_grid_partition) - coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) - mh = ModelHierarchy(parts,coarse_model,num_parts_x_level) - - # Create Operators: - order = 1 - u(x) = 1.0 - reffe = ReferenceFE(lagrangian,Float64,order) - - tests = TestFESpace(mh,reffe;dirichlet_tags="boundary") - trials = TrialFESpace(tests,u) - - qdegree = order*2+1 - ops1 = setup_transfer_operators(trials, qdegree; restriction_method=:projection, mode=:solution) - restrictions1, prolongations1 = ops1 - ops2 = setup_transfer_operators(trials, qdegree; restriction_method=:interpolation, mode=:solution) - restrictions2, prolongations2 = ops2 - ops3 = setup_transfer_operators(trials, qdegree; restriction_method=:dof_mask, mode=:solution) - restrictions3, prolongations3 = ops3 - - a(u,v,dΩ) = ∫(v⋅u)*dΩ - l(v,dΩ) = ∫(v⋅u)*dΩ - mats, A, b = compute_hierarchy_matrices(trials,a,l,qdegree) - - for lev in 1:num_levels-1 - parts_h = get_level_parts(mh,lev) - parts_H = get_level_parts(mh,lev+1) - - if i_am_in(parts_h) - i_am_main(parts_h) && println("Lev : ", lev) - Ah = mats[lev] - xh = pfill(1.0,partition(axes(Ah,2))) - yh1 = pfill(0.0,partition(axes(Ah,2))) - yh2 = pfill(0.0,partition(axes(Ah,2))) - yh3 = pfill(0.0,partition(axes(Ah,2))) - - if i_am_in(parts_H) - AH = mats[lev+1] - xH = pfill(1.0,partition(axes(AH,2))) - yH1 = pfill(0.0,partition(axes(AH,2))) - yH2 = pfill(0.0,partition(axes(AH,2))) - yH3 = pfill(0.0,partition(axes(AH,2))) - else - xH = nothing - yH1 = nothing - yH2 = nothing - yH3 = nothing - end - - # ---- Restriction ---- - i_am_main(parts_h) && println(" > Restriction") - R1 = restrictions1[lev] - mul!(yH1,R1,xh) - - R2 = restrictions2[lev] - mul!(yH2,R2,xh) - - R3 = restrictions3[lev] - mul!(yH3,R3,xh) - - if i_am_in(parts_H) - y_ref = pfill(1.0,partition(axes(AH,2))) - tests = map(own_values(y_ref),own_values(yH1),own_values(yH2),own_values(yH3)) do y_ref,y1,y2,y3 - map(y -> norm(y-y_ref) < 1.e-3 ,[y1,y2,y3]) - end - @test all(PartitionedArrays.getany(tests)) - end - - # ---- Prolongation ---- - i_am_main(parts_h) && println(" > Prolongation") - P1 = prolongations1[lev] - mul!(yh1,P1,xH) - - P2 = prolongations2[lev] - mul!(yh2,P2,xH) - - P3 = prolongations3[lev] - mul!(yh3,P3,xH) - - y_ref = pfill(1.0,partition(axes(Ah,2))) - tests = map(own_values(y_ref),own_values(yh1),own_values(yh2),own_values(yh3)) do y_ref,y1,y2,y3 - map(y -> norm(y-y_ref) < 1.e-3 ,[y1,y2,y3]) - end - @test all(PartitionedArrays.getany(tests)) - - end - end - end -end - -num_parts_x_level = [4,2,2] # Procs in each refinement level -num_trees = (1,1) # Number of initial P4est trees -num_refs_coarse = 2 # Number of initial refinements - -num_ranks = num_parts_x_level[1] - -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -run(parts,num_parts_x_level,num_trees,num_refs_coarse) - -println("AT THE END") -MPI.Finalize() -end diff --git a/test/mpi/MUMPSSolversTests.jl b/test/mpi/MUMPSSolversTests.jl deleted file mode 100644 index beaa8905..00000000 --- a/test/mpi/MUMPSSolversTests.jl +++ /dev/null @@ -1,83 +0,0 @@ -module MUMPSSolversTests - -using Test -using MPI -using Gridap -using GridapDistributed -using PartitionedArrays -using IterativeSolvers - -using GridapSolvers -using GridapSolvers.LinearSolvers - -using GridapPETSc - -function set_ksp_options(ksp) - pc = Ref{GridapPETSc.PETSC.PC}() - mumpsmat = Ref{GridapPETSc.PETSC.Mat}() - @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL) - @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPPREONLY) - @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc) - @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCLU) - @check_error_code GridapPETSc.PETSC.PCFactorSetMatSolverType(pc[],GridapPETSc.PETSC.MATSOLVERMUMPS) - @check_error_code GridapPETSc.PETSC.PCFactorSetUpMatSolverType(pc[]) - @check_error_code GridapPETSc.PETSC.PCFactorGetMatrix(pc[],mumpsmat) - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 4, 1) - # percentage increase in the estimated working space - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 14, 1000) - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 28, 2) - @check_error_code GridapPETSc.PETSC.MatMumpsSetIcntl(mumpsmat[], 29, 2) - @check_error_code GridapPETSc.PETSC.MatMumpsSetCntl(mumpsmat[], 3, 1.0e-6) -end - -function main(parts,nranks,domain_partition) - GridapPETSc.with() do - domain = (0,1,0,1) - model = CartesianDiscreteModel(parts,nranks,domain,domain_partition) - - sol(x) = x[1] + x[2] - f(x) = -Δ(sol)(x) - - order = 1 - qorder = order*2 + 1 - reffe = ReferenceFE(lagrangian,Float64,order) - Vh = TestFESpace(model,reffe,conformity=:H1,dirichlet_tags="boundary") - Uh = TrialFESpace(Vh,sol) - u = interpolate(sol,Uh) - - Ω = Triangulation(model) - dΩ = Measure(Ω,qorder) - a(u,v) = ∫(∇(v)⋅∇(u))*dΩ - l(v) = ∫(v⋅f)*dΩ - - op = AffineFEOperator(a,l,Uh,Vh) - A, b = get_matrix(op), get_vector(op) - - P = PETScLinearSolver(set_ksp_options) - ss = symbolic_setup(P,A) - ns = numerical_setup(ss,A) - - x = pfill(0.0,partition(axes(A,2))) - solve!(x,ns,b) - - u = interpolate(sol,Uh) - uh = FEFunction(Uh,x) - eh = uh - u - E = sum(∫(eh*eh)*dΩ) - if i_am_main(parts) - println("L2 Error: ", E) - end - - @test E < 1.e-8 - end -end - -domain_partition = (32,32) -num_ranks = (2,2) -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -main(parts,num_ranks,domain_partition) -MPI.Finalize() - -end \ No newline at end of file diff --git a/test/mpi/RedistributeToolsTests.jl b/test/mpi/RedistributeToolsTests.jl deleted file mode 100644 index 3a57e801..00000000 --- a/test/mpi/RedistributeToolsTests.jl +++ /dev/null @@ -1,90 +0,0 @@ -module RedistributeToolsTests -using MPI -using PartitionedArrays -using Gridap -using GridapDistributed -using GridapP4est -using Test - -using GridapSolvers -using GridapSolvers.MultilevelTools -using GridapDistributed: redistribute_fe_function - -function run(parts,num_parts_x_level,coarse_grid_partition,num_refs_coarse) - GridapP4est.with(parts) do - domain = (0,1,0,1) - num_levels = length(num_parts_x_level) - cparts = generate_subparts(parts,num_parts_x_level[num_levels]) - cmodel = CartesianDiscreteModel(domain,coarse_grid_partition) - coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) - mh = ModelHierarchy(parts,coarse_model,num_parts_x_level) - - level_parts = get_level_parts(mh) - old_parts = level_parts[2] - new_parts = level_parts[1] - - # FE Spaces - order = 2 - u(x) = x[1]^2 + x[2]^2 - 3.0*x[1]*x[2] - reffe = ReferenceFE(lagrangian,Float64,order) - glue = mh.levels[1].red_glue - - model_old = get_model_before_redist(mh.levels[1]) - if i_am_in(old_parts) - VOLD = TestFESpace(model_old,reffe,dirichlet_tags="boundary") - UOLD = TrialFESpace(VOLD,u) - else - VOLD = nothing - UOLD = nothing - end - - model_new = get_model(mh.levels[1]) - VNEW = TestFESpace(model_new,reffe,dirichlet_tags="boundary") - UNEW = TrialFESpace(VNEW,u) - - # Triangulations - qdegree = 2*order+1 - Ω_new = Triangulation(model_new) - dΩ_new = Measure(Ω_new,qdegree) - uh_new = interpolate(u,UNEW) - - if i_am_in(old_parts) - Ω_old = Triangulation(model_old) - dΩ_old = Measure(Ω_old,qdegree) - uh_old = interpolate(u,UOLD) - else - Ω_old = nothing - dΩ_old = nothing - uh_old = nothing - end - - # Old -> New - uh_old_red = redistribute_fe_function(uh_old,UNEW,model_new,glue) - n = sum(∫(uh_old_red)*dΩ_new) - if i_am_in(old_parts) - o = sum(∫(uh_old)*dΩ_old) - @test o ≈ n - end - - # New -> Old - uh_new_red = redistribute_fe_function(uh_new,UOLD,model_old,glue;reverse=true) - n = sum(∫(uh_new)*dΩ_new) - if i_am_in(old_parts) - o = sum(∫(uh_new_red)*dΩ_old) - @test o ≈ n - end - end -end - - -num_parts_x_level = [4,2] # Procs in each refinement level -num_trees = (1,1) # Number of initial P4est trees -num_refs_coarse = 3 # Number of initial refinements - -num_ranks = num_parts_x_level[1] -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -run(parts,num_parts_x_level,num_trees,num_refs_coarse) -MPI.Finalize() -end diff --git a/test/mpi/RefinementToolsTests.jl b/test/mpi/RefinementToolsTests.jl deleted file mode 100644 index 4fb128fb..00000000 --- a/test/mpi/RefinementToolsTests.jl +++ /dev/null @@ -1,104 +0,0 @@ -module RefinementToolsTests -using MPI -using PartitionedArrays -using Gridap -using GridapDistributed -using GridapP4est -using Test -using IterativeSolvers - -using GridapSolvers -using GridapSolvers.MultilevelTools - -function run(parts,num_parts_x_level,coarse_grid_partition,num_refs_coarse) - GridapP4est.with(parts) do - domain = (0,1,0,1) - num_levels = length(num_parts_x_level) - cparts = generate_subparts(parts,num_parts_x_level[num_levels]) - cmodel = CartesianDiscreteModel(domain,coarse_grid_partition) - coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) - mh = ModelHierarchy(parts,coarse_model,num_parts_x_level) - - # FE Spaces - order = 2 - sol(x) = x[1]^2 + x[2]^2 - 3.0*x[1]*x[2] - reffe = ReferenceFE(lagrangian,Float64,order) - tests = TestFESpace(mh,reffe;conformity=:H1,dirichlet_tags="boundary") - trials = TrialFESpace(tests,sol) - - quad_order = 2*order+1 - for lev in 1:num_levels-1 - fparts = get_level_parts(mh,lev) - cparts = get_level_parts(mh,lev+1) - - if i_am_in(cparts) - model_h = get_model_before_redist(mh,lev) - Vh = get_fe_space_before_redist(tests,lev) - Uh = get_fe_space_before_redist(trials,lev) - Ωh = get_triangulation(model_h) - dΩh = Measure(Ωh,quad_order) - uh = interpolate(sol,Uh) - - model_H = get_model(mh,lev+1) - VH = get_fe_space(tests,lev+1) - UH = get_fe_space(trials,lev+1) - ΩH = get_triangulation(model_H) - dΩH = Measure(ΩH,quad_order) - uH = interpolate(sol,UH) - dΩhH = Measure(ΩH,Ωh,quad_order) - - # Coarse FEFunction -> Fine FEFunction, by projection - ah(u,v) = ∫(v⋅u)*dΩh - lh(v) = ∫(v⋅uH)*dΩh - oph = AffineFEOperator(ah,lh,Uh,Vh) - Ah = get_matrix(oph) - bh = get_vector(oph) - - xh = pfill(0.0,partition(axes(Ah,2))) - IterativeSolvers.cg!(xh,Ah,bh;verbose=i_am_main(parts),reltol=1.0e-08) - uH_projected = FEFunction(Uh,xh) - - _eh = uh-uH_projected - eh = sum(∫(_eh⋅_eh)*dΩh) - i_am_main(parts) && println("Error H2h: ", eh) - @test eh < 1.0e-10 - - # Fine FEFunction -> Coarse FEFunction, by projection - aH(u,v) = ∫(v⋅u)*dΩH - lH(v) = ∫(v⋅uH_projected)*dΩhH - opH = AffineFEOperator(aH,lH,UH,VH) - AH = get_matrix(opH) - bH = get_vector(opH) - - xH = pfill(0.0,partition(axes(AH,2))) - IterativeSolvers.cg!(xH,AH,bH;verbose=i_am_main(parts),reltol=1.0e-08) - uh_projected = FEFunction(UH,xH) - - _eH = uH-uh_projected - eH = sum(∫(_eH⋅_eH)*dΩH) - i_am_main(parts) && println("Error h2H: ", eH) - @test eh < 1.0e-10 - - # Coarse FEFunction -> Fine FEFunction, by interpolation - uH_i = interpolate(uH,Uh) - - _eh = uH_i-uh - eh = sum(∫(_eh⋅_eh)*dΩh) - i_am_main(parts) && println("Error h2H: ", eh) - @test eh < 1.0e-10 - end - end - end -end - -num_parts_x_level = [4,2,2] # Procs in each refinement level -num_trees = (1,1) # Number of initial P4est trees -num_refs_coarse = 2 # Number of initial refinements - -num_ranks = num_parts_x_level[1] -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -run(parts,num_parts_x_level,num_trees,num_refs_coarse) -MPI.Finalize() -end diff --git a/test/mpi/RestrictDofsTests.jl b/test/mpi/RestrictDofsTests.jl deleted file mode 100644 index c87c92ac..00000000 --- a/test/mpi/RestrictDofsTests.jl +++ /dev/null @@ -1,100 +0,0 @@ -module RestrictDofsTests -using MPI -using Test -using LinearAlgebra -using IterativeSolvers -using FillArrays - -using Gridap -using Gridap.ReferenceFEs -using PartitionedArrays -using GridapDistributed -using GridapP4est - -using GridapSolvers -using GridapSolvers.LinearSolvers - - -u(x) = x[1] + x[2] -f(x) = -Δ(u)(x) - -function main(parts, coarse_grid_partition, num_parts_x_level, num_refs_coarse, order) - GridapP4est.with(parts) do - domain = (0,1,0,1) - num_levels = length(num_parts_x_level) - cparts = generate_subparts(parts,num_parts_x_level[num_levels]) - cmodel = CartesianDiscreteModel(domain,coarse_grid_partition) - coarse_model = OctreeDistributedDiscreteModel(cparts,cmodel,num_refs_coarse) - mh = ModelHierarchy(parts,coarse_model,num_parts_x_level) - - qdegree = 2*(order+1) - reffe = ReferenceFE(lagrangian,Float64,order) - tests = TestFESpace(mh,reffe;conformity=:H1,dirichlet_tags="boundary") - trials = TrialFESpace(tests,u) - - biform(u,v,dΩ) = ∫(∇(v)⋅∇(u))dΩ - liform(v,dΩ) = ∫(v*f)dΩ - smatrices, A, b = compute_hierarchy_matrices(trials,biform,liform,qdegree) - - # Preconditioner - smoothers = Fill(RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0),num_levels-1) - restrictions, prolongations = setup_transfer_operators(trials,qdegree;mode=:residual,restriction_method=:dof_mask) - - gmg = GMGLinearSolver(mh, - smatrices, - prolongations, - restrictions, - pre_smoothers=smoothers, - post_smoothers=smoothers, - maxiter=1, - rtol=1.0e-10, - verbose=false, - mode=:preconditioner) - ss = symbolic_setup(gmg,A) - ns = numerical_setup(ss,A) - - # Solve - x = pfill(0.0,partition(axes(A,2))) - x, history = IterativeSolvers.cg!(x,A,b; - verbose=i_am_main(parts), - reltol=1.0e-12, - Pl=ns, - log=true) - - # Error norms and print solution - model = get_model(mh,1) - Uh = get_fe_space(trials,1) - Ω = Triangulation(model) - dΩ = Measure(Ω,qdegree) - uh = FEFunction(Uh,x) - e = u-uh - e_l2 = sum(∫(e⋅e)dΩ) - tol = 1.0e-9 - @test e_l2 < tol - if i_am_main(parts) - println("L2 error = ", e_l2) - end - end -end - -############################################## - -if !MPI.Initialized() - MPI.Init() -end - -# Parameters -order = 1 -coarse_grid_partition = (2,2) -num_refs_coarse = 2 - -num_parts_x_level = [4,2,1] -num_ranks = num_parts_x_level[1] -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -main(parts,coarse_grid_partition,num_parts_x_level,num_refs_coarse,order) - - -MPI.Finalize() -end diff --git a/test/mpi/RichardsonSmoothersTests.jl b/test/mpi/RichardsonSmoothersTests.jl deleted file mode 100644 index 16074250..00000000 --- a/test/mpi/RichardsonSmoothersTests.jl +++ /dev/null @@ -1,68 +0,0 @@ -module RichardsonSmoothersTests - -using Test -using MPI -using Gridap -using GridapDistributed -using PartitionedArrays -using IterativeSolvers -using GridapP4est - -using GridapSolvers -using GridapSolvers.LinearSolvers - -function main(parts,nranks,domain_partition) - GridapP4est.with(parts) do - domain = (0,1,0,1) - model = CartesianDiscreteModel(parts,nranks,domain,domain_partition) - - sol(x) = x[1] + x[2] - f(x) = -Δ(sol)(x) - - order = 1 - qorder = order*2 + 1 - reffe = ReferenceFE(lagrangian,Float64,order) - Vh = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary") - Uh = TrialFESpace(Vh,sol) - u = interpolate(sol,Uh) - - Ω = Triangulation(model) - dΩ = Measure(Ω,qorder) - a(u,v) = ∫(∇(v)⋅∇(u))*dΩ - l(v) = ∫(v⋅f)*dΩ - - op = AffineFEOperator(a,l,Uh,Vh) - A, b = get_matrix(op), get_vector(op) - - P = RichardsonSmoother(JacobiLinearSolver(),10,2.0/3.0) - ss = symbolic_setup(P,A) - ns = numerical_setup(ss,A) - - x = pfill(1.0,partition(axes(A,2))) - x, history = IterativeSolvers.cg!(x,A,b; - verbose=i_am_main(parts), - reltol=1.0e-8, - Pl=ns, - log=true) - - u = interpolate(sol,Uh) - uh = FEFunction(Uh,x) - eh = uh - u - E = sum(∫(eh*eh)*dΩ) - if i_am_main(parts) - println("L2 Error: ", E) - end - - @test E < 1.e-8 - end -end - -domain_partition = (32,32) -num_ranks = (2,2) -parts = with_mpi() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end -main(parts,num_ranks,domain_partition) -MPI.Finalize() - -end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 4fd2e586..8b374075 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,72 +1,12 @@ using GridapSolvers using Test -using ArgParse -using MPI -function parse_commandline() - s = ArgParseSettings() - @add_arg_table! s begin - "--image-file", "-i" - help = "Path to the image file that one can use in order to accelerate MPI tests" - arg_type = String - default="GridapDistributed.so" - end - return parse_args(s) +@testset "Sequential tests" begin + include("MultilevelTools/seq/runtests.jl") + include("LinearSolvers/seq/runtests.jl") end -""" - run_tests(testdir) -""" -function run_tests(testdir) - parsed_args = parse_commandline() - image_file_path=parsed_args["image-file"] - image_file_exists=isfile(image_file_path) - - nprocs_str = get(ENV, "JULIA_GRIDAP_SOLVERS_TEST_NPROCS","") - nprocs = nprocs_str == "" ? clamp(Sys.CPU_THREADS, 2, 4) : parse(Int, nprocs_str) - istest(f) = endswith(f, ".jl") && !(f=="runtests.jl") - testfiles = sort(filter(istest, readdir(testdir))) - @time @testset "$f" for f in testfiles - MPI.mpiexec() do cmd - if f in ["DistributedGridTransferOperatorsTests.jl", - "RedistributeToolsTests.jl", - "RefinementToolsTests.jl", - "RichardsonSmoothersTests.jl", - "ModelHierarchiesTests.jl", - "GMGLinearSolversPoissonTests.jl", - "GMGLinearSolversLaplacianTests.jl", - "GMGLinearSolversVectorLaplacianTests.jl", - "GMGLinearSolversHDivRTTests.jl", - "MUMPSSolversTests.jl", - "GMGLinearSolversMUMPSTests.jl", - "RestrictDofsTests.jl", - "PRefinementGMGLinearSolversPoissonTests.jl"] - np = 4 - extra_args = "" - else - np = 4 # nprocs - extra_args = "" - end - if ! image_file_exists - cmd = `$cmd -n $(np) --allow-run-as-root --oversubscribe $(Base.julia_cmd()) --project=. $(joinpath(testdir, f))` - else - cmd = `$cmd -n $(np) --allow-run-as-root --oversubscribe $(Base.julia_cmd()) -J$(image_file_path) --project=. $(joinpath(testdir, f)) $(split(extra_args))` - end - @show cmd - run(cmd) - @test true - end - end +@testset "MPI tests" begin + include("MultilevelTools/mpi/runtests.jl") + include("LinearSolvers/mpi/runtests.jl") end - -# MPI tests -run_tests(joinpath(@__DIR__, "mpi")) - -# Sequential tests -@time @testset "BlockDiagonalSmoothersTests" begin include("seq/BlockDiagonalSmoothersTests.jl") end -@time @testset "DistributedPatchFESpacesTests" begin include("seq/DistributedPatchFESpacesTests.jl") end -@time @testset "KrylovSolversTests" begin include("seq/KrylovSolversTests.jl") end -@time @testset "IterativeSolversTests" begin include("seq/IterativeSolversTests.jl") end -@time @testset "PatchLinearSolverTests" begin include("seq/PatchLinearSolverTests.jl") end -@time @testset "SymGaussSeidelSmoothersTests" begin include("seq/SymGaussSeidelSmoothersTests.jl") end -@time @testset "SchurComplementSolversTests" begin include("seq/SchurComplementSolversTests.jl") end diff --git a/test/seq/IterativeSolversTests.jl b/test/seq/IterativeSolversTests.jl deleted file mode 100644 index 6177429f..00000000 --- a/test/seq/IterativeSolversTests.jl +++ /dev/null @@ -1,97 +0,0 @@ -module IterativeSolversTests - -using Test -using Gridap -using IterativeSolvers -using LinearAlgebra -using SparseArrays -using PartitionedArrays - -using GridapSolvers -using GridapSolvers.LinearSolvers - -sol(x) = x[1] + x[2] -f(x) = -Δ(sol)(x) - -function l2_error(x,Uh,dΩ) - u = interpolate(sol,Uh) - uh = FEFunction(Uh,x) - eh = uh - u - return sum(∫(eh*eh)*dΩ) -end - -function main(model,is_distributed) - order = 1 - qorder = order*2 + 1 - reffe = ReferenceFE(lagrangian,Float64,order) - Vh = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary") - Uh = TrialFESpace(Vh,sol) - u = interpolate(sol,Uh) - - Ω = Triangulation(model) - dΩ = Measure(Ω,qorder) - a(u,v) = ∫(∇(v)⋅∇(u))*dΩ - l(v) = ∫(v⋅f)*dΩ - - op = AffineFEOperator(a,l,Uh,Vh) - A, b = get_matrix(op), get_vector(op); - - # CG - solver = IS_ConjugateGradientSolver(;maxiter=100,reltol=1.e-12) - ss = symbolic_setup(solver,A) - ns = numerical_setup(ss,A) - - x = LinearSolvers.allocate_col_vector(A) - y = copy(b) - solve!(x,ns,y) - @test l2_error(x,Uh,dΩ) < 1.e-8 - - # SSOR - solver = IS_SSORSolver(2.0/3.0;maxiter=100) - ss = symbolic_setup(solver,A) - ns = numerical_setup(ss,A) - - x = LinearSolvers.allocate_row_vector(A) - y = copy(b) - solve!(x,ns,y) - !is_distributed && (@test l2_error(x,Uh,dΩ) < 1.e-8) - - if !is_distributed - # GMRES - solver = IS_GMRESSolver(;maxiter=100,reltol=1.e-12) - ss = symbolic_setup(solver,A) - ns = numerical_setup(ss,A) - - x = LinearSolvers.allocate_row_vector(A) - y = copy(b) - solve!(x,ns,y) - @test l2_error(x,Uh,dΩ) < 1.e-8 - - # MINRES - solver = IS_MINRESSolver(;maxiter=100,reltol=1.e-12) - ss = symbolic_setup(solver,A) - ns = numerical_setup(ss,A) - - x = LinearSolvers.allocate_row_vector(A) - y = copy(b) - solve!(x,ns,y) - @test l2_error(x,Uh,dΩ) < 1.e-8 - end -end - -# Completely serial -mesh_partition = (8,8) -domain = (0,1,0,1) -model = CartesianDiscreteModel(domain,mesh_partition) -main(model,false) - -# Sequential -num_ranks = (1,2) -parts = with_debug() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end - -model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) -main(model,true) - -end \ No newline at end of file diff --git a/test/seq/SolverInterfacesTests.jl b/test/seq/SolverInterfacesTests.jl deleted file mode 100644 index 679c34fe..00000000 --- a/test/seq/SolverInterfacesTests.jl +++ /dev/null @@ -1,37 +0,0 @@ - -using Gridap - -using GridapSolvers -import GridapSolvers.SolverInterfaces as SI -using GridapSolvers.LinearSolvers - -sol(x) = x[1] + x[2] -f(x) = -Δ(sol)(x) - -mesh_partition = (10,10) -domain = (0,1,0,1) -model = CartesianDiscreteModel(domain,mesh_partition) - -order = 1 -qorder = order*2 + 1 -reffe = ReferenceFE(lagrangian,Float64,order) -Vh = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary") -Uh = TrialFESpace(Vh,sol) -u = interpolate(sol,Uh) - -Ω = Triangulation(model) -dΩ = Measure(Ω,qorder) -a(u,v) = ∫(∇(v)⋅∇(u))*dΩ -l(v) = ∫(v⋅f)*dΩ -op = AffineFEOperator(a,l,Uh,Vh) -A, b = get_matrix(op), get_vector(op); -P = JacobiLinearSolver() - -solver = LinearSolvers.GMRESSolver(10;Pl=P,rtol=1.e-8,verbose=2) -ns = numerical_setup(symbolic_setup(solver,A),A) -x = LinearSolvers.allocate_col_vector(A) -solve!(x,ns,b) - - - -using AbstractTrees diff --git a/test/seq/SymGaussSeidelSmoothersTests.jl b/test/seq/SymGaussSeidelSmoothersTests.jl deleted file mode 100644 index e29a7864..00000000 --- a/test/seq/SymGaussSeidelSmoothersTests.jl +++ /dev/null @@ -1,65 +0,0 @@ -module SymGaussSeidelSmoothersTests - -using Test -using MPI -using Gridap -using GridapDistributed -using PartitionedArrays -using IterativeSolvers - -using GridapSolvers -using GridapSolvers.LinearSolvers - -sol(x) = x[1] + x[2] -f(x) = -Δ(sol)(x) - -function main(model) - order = 1 - qorder = order*2 + 1 - reffe = ReferenceFE(lagrangian,Float64,order) - Vh = TestFESpace(model,reffe;conformity=:H1,dirichlet_tags="boundary") - Uh = TrialFESpace(Vh,sol) - u = interpolate(sol,Uh) - - Ω = Triangulation(model) - dΩ = Measure(Ω,qorder) - a(u,v) = ∫(∇(v)⋅∇(u))*dΩ - l(v) = ∫(v⋅f)*dΩ - - op = AffineFEOperator(a,l,Uh,Vh) - A, b = get_matrix(op), get_vector(op); - - P = SymGaussSeidelSmoother(10) - ss = symbolic_setup(P,A) - ns = numerical_setup(ss,A) - - x = LinearSolvers.allocate_col_vector(A) - x, history = IterativeSolvers.cg!(x,A,b; - verbose=true, - reltol=1.0e-8, - Pl=ns, - log=true); - - u = interpolate(sol,Uh) - uh = FEFunction(Uh,x) - eh = uh - u - E = sum(∫(eh*eh)*dΩ) - return E < 1.e-8 -end - -# Completely serial -mesh_partition = (8,8) -domain = (0,1,0,1) -model = CartesianDiscreteModel(domain,mesh_partition) -@test main(model) - -# Sequential -num_ranks = (1,2) -parts = with_debug() do distribute - distribute(LinearIndices((prod(num_ranks),))) -end - -model = CartesianDiscreteModel(parts,num_ranks,domain,mesh_partition) -@test main(model) - -end \ No newline at end of file