From eb05fff0efc838b9597a956f0014a4ac478ba6b7 Mon Sep 17 00:00:00 2001 From: JordiManyer Date: Tue, 12 Dec 2023 16:35:12 +1100 Subject: [PATCH] Added ADS solver --- src/LinearSolvers/LinearSolvers.jl | 6 +- src/LinearSolvers/PETSc/HipmairXuSolvers.jl | 59 ++++++++++++++- src/LinearSolvers/PETSc/PETScUtils.jl | 81 +++++++++++++++++++++ test/_dev/PETSc/HipmairXuHDiv.jl | 28 +++---- 4 files changed, 152 insertions(+), 22 deletions(-) diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl index 40573a67..5eb500ad 100644 --- a/src/LinearSolvers/LinearSolvers.jl +++ b/src/LinearSolvers/LinearSolvers.jl @@ -9,10 +9,7 @@ using BlockArrays using IterativeSolvers using Gridap -using Gridap.Helpers -using Gridap.Algebra -using Gridap.FESpaces -using Gridap.MultiField +using Gridap.Helpers, Gridap.Algebra, Gridap.CellData, Gridap.Arrays, Gridap.FESpaces, Gridap.MultiField using PartitionedArrays using GridapPETSc @@ -45,6 +42,7 @@ include("Krylov/GMRESSolvers.jl") include("Krylov/FGMRESSolvers.jl") include("Krylov/MINRESSolvers.jl") +include("PETSc/PETScUtils.jl") include("PETSc/ElasticitySolvers.jl") include("PETSc/HipmairXuSolvers.jl") diff --git a/src/LinearSolvers/PETSc/HipmairXuSolvers.jl b/src/LinearSolvers/PETSc/HipmairXuSolvers.jl index ca3366a4..7d8b6395 100644 --- a/src/LinearSolvers/PETSc/HipmairXuSolvers.jl +++ b/src/LinearSolvers/PETSc/HipmairXuSolvers.jl @@ -1,9 +1,60 @@ -""" +function ADS_Solver(model,tags,order;rtol=1e-8,maxits=300) + @assert (num_cell_dims(model) == 3) "Not implemented for 2D" + @assert (order == 1) "Only works for linear order" -""" + V_H1_sc, V_H1, V_Hdiv, V_Hcurl = get_ads_spaces(model,order,tags) + G, C, Π_div, Π_curl = get_ads_operators(V_H1_sc,V_H1,V_Hcurl,V_Hdiv) -struct HypreADSSolver end + D = num_cell_dims(model) + ksp_setup(ksp) = ads_ksp_setup(ksp,rtol,maxits,D,G,C,Π_div,Π_curl) + return PETScLinearSolver(ksp_setup) +end -struct HypreAMSSolver end +function get_ads_spaces(model,order,tags) + reffe_H1_sc = ReferenceFE(lagrangian,Float64,order) + V_H1_sc = FESpace(model,reffe_H1_sc;dirichlet_tags=tags) + + reffe_H1 = ReferenceFE(lagrangian,VectorValue{D,Float64},order) + V_H1 = FESpace(model,reffe_H1;dirichlet_tags=tags) + + reffe_Hdiv = ReferenceFE(raviart_thomas,Float64,order-1) + V_Hdiv = FESpace(model,reffe_Hdiv;dirichlet_tags=tags) + + reffe_Hcurl = ReferenceFE(nedelec,Float64,order-1) + V_Hcurl = FESpace(model,reffe_Hcurl;dirichlet_tags=tags) + return V_H1_sc, V_H1, V_Hdiv, V_Hcurl +end + +function get_ads_operators(V_H1_sc,V_H1_vec,V_Hcurl,V_Hdiv) + G = interpolation_operator(u->∇(u),V_H1_sc,V_Hcurl) + C = interpolation_operator(u->cross(∇,u),V_Hcurl,V_Hdiv) + Π_div = interpolation_operator(u->u,V_H1_vec,V_Hdiv) + Π_curl = interpolation_operator(u->u,V_H1_vec,V_Hcurl) + return G, C, Π_div, Π_curl +end + +function ads_ksp_setup(ksp,rtol,maxits,dim,G,C,Π_div,Π_curl) + rtol = PetscScalar(rtol) + atol = GridapPETSc.PETSC.PETSC_DEFAULT + dtol = GridapPETSc.PETSC.PETSC_DEFAULT + maxits = PetscInt(maxits) + + @check_error_code GridapPETSc.PETSC.KSPSetType(ksp[],GridapPETSc.PETSC.KSPGMRES) + @check_error_code GridapPETSc.PETSC.KSPSetTolerances(ksp[], rtol, atol, dtol, maxits) + + pc = Ref{GridapPETSc.PETSC.PC}() + @check_error_code GridapPETSc.PETSC.KSPGetPC(ksp[],pc) + @check_error_code GridapPETSc.PETSC.PCSetType(pc[],GridapPETSc.PETSC.PCHYPRE) + + _G = convert(PETScMatrix,G) + _C = convert(PETScMatrix,C) + _Π_div = convert(PETScMatrix,Π_div) + _Π_curl = convert(PETScMatrix,Π_curl) + @check_error_code GridapPETSc.PETSC.PCHYPRESetDiscreteGradient(pc[],_G.mat[]) + @check_error_code GridapPETSc.PETSC.PCHYPRESetDiscreteCurl(pc[],_C.mat[]) + @check_error_code GridapPETSc.PETSC.PCHYPRESetInterpolations(pc[],dim,_Π_div.mat[],C_NULL,_Π_curl.mat[],C_NULL) + + @check_error_code GridapPETSc.PETSC.KSPView(ksp[],C_NULL) +end \ No newline at end of file diff --git a/src/LinearSolvers/PETSc/PETScUtils.jl b/src/LinearSolvers/PETSc/PETScUtils.jl index eccbbe5f..1b8de887 100644 --- a/src/LinearSolvers/PETSc/PETScUtils.jl +++ b/src/LinearSolvers/PETSc/PETScUtils.jl @@ -54,3 +54,84 @@ end # Interpolation matrices +function interpolation_operator(op,U_in,V_out; + strat=SubAssembledRows(), + Tm=SparseMatrixCSR{0,PetscScalar,PetscInt}, + Tv=Vector{PetscScalar}) + out_dofs = get_fe_dof_basis(V_out) + in_basis = get_fe_basis(U_in) + + cell_interp_mats = out_dofs(op(in_basis)) + local_contr = map(local_views(out_dofs),cell_interp_mats) do dofs, arr + contr = DomainContribution() + add_contribution!(contr,get_triangulation(dofs),arr) + return contr + end + contr = GridapDistributed.DistributedDomainContribution(local_contr) + + matdata = collect_cell_matrix(U_in,V_out,contr) + assem = SparseMatrixAssembler(Tm,Tv,U_in,V_out,strat) + + I = allocate_matrix(assem,matdata) + takelast_matrix!(I,assem,matdata) + return I +end + +function takelast_matrix(a::SparseMatrixAssembler,matdata) + m1 = Gridap.Algebra.nz_counter(get_matrix_builder(a),(get_rows(a),get_cols(a))) + symbolic_loop_matrix!(m1,a,matdata) + m2 = Gridap.Algebra.nz_allocation(m1) + takelast_loop_matrix!(m2,a,matdata) + m3 = Gridap.Algebra.create_from_nz(m2) + return m3 +end + +function takelast_matrix!(mat,a::SparseMatrixAssembler,matdata) + LinearAlgebra.fillstored!(mat,zero(eltype(mat))) + takelast_matrix_add!(mat,a,matdata) +end + +function takelast_matrix_add!(mat,a::SparseMatrixAssembler,matdata) + takelast_loop_matrix!(mat,a,matdata) + Gridap.Algebra.create_from_nz(mat) +end + +function takelast_loop_matrix!(A,a::GridapDistributed.DistributedSparseMatrixAssembler,matdata) + rows = get_rows(a) + cols = get_cols(a) + map(takelast_loop_matrix!,local_views(A,rows,cols),local_views(a),matdata) +end + +function takelast_loop_matrix!(A,a::SparseMatrixAssembler,matdata) + strategy = Gridap.FESpaces.get_assembly_strategy(a) + for (cellmat,_cellidsrows,_cellidscols) in zip(matdata...) + cellidsrows = Gridap.FESpaces.map_cell_rows(strategy,_cellidsrows) + cellidscols = Gridap.FESpaces.map_cell_cols(strategy,_cellidscols) + @assert length(cellidscols) == length(cellidsrows) + @assert length(cellmat) == length(cellidsrows) + if length(cellmat) > 0 + rows_cache = array_cache(cellidsrows) + cols_cache = array_cache(cellidscols) + vals_cache = array_cache(cellmat) + mat1 = getindex!(vals_cache,cellmat,1) + rows1 = getindex!(rows_cache,cellidsrows,1) + cols1 = getindex!(cols_cache,cellidscols,1) + add! = Gridap.Arrays.AddEntriesMap((a,b) -> b) + add_cache = return_cache(add!,A,mat1,rows1,cols1) + caches = add_cache, vals_cache, rows_cache, cols_cache + _takelast_loop_matrix!(A,caches,cellmat,cellidsrows,cellidscols) + end + end + A +end + +@noinline function _takelast_loop_matrix!(mat,caches,cell_vals,cell_rows,cell_cols) + add_cache, vals_cache, rows_cache, cols_cache = caches + add! = Gridap.Arrays.AddEntriesMap((a,b) -> b) + for cell in 1:length(cell_cols) + rows = getindex!(rows_cache,cell_rows,cell) + cols = getindex!(cols_cache,cell_cols,cell) + vals = getindex!(vals_cache,cell_vals,cell) + evaluate!(add_cache,add!,mat,vals,rows,cols) + end +end diff --git a/test/_dev/PETSc/HipmairXuHDiv.jl b/test/_dev/PETSc/HipmairXuHDiv.jl index fc9dcb0b..dd76d834 100644 --- a/test/_dev/PETSc/HipmairXuHDiv.jl +++ b/test/_dev/PETSc/HipmairXuHDiv.jl @@ -8,15 +8,15 @@ using LinearAlgebra using Gridap.Geometry, Gridap.FESpaces, Gridap.CellData, Gridap.Arrays -function get_operators(V_H1_sc,V_H1_vec,V_Hcurl,V_Hdiv,trian) - G = interpolation_operator(u->∇(u),V_H1_sc,V_Hcurl,trian) - C = interpolation_operator(u->cross(∇,u),V_Hcurl,V_Hdiv,trian) - Π_div = interpolation_operator(u->u,V_H1_vec,V_Hdiv,trian) - Π_curl = interpolation_operator(u->u,V_H1_vec,V_Hcurl,trian) +function get_operators(V_H1_sc,V_H1_vec,V_Hcurl,V_Hdiv) + G = interpolation_operator(u->∇(u),V_H1_sc,V_Hcurl) + C = interpolation_operator(u->cross(∇,u),V_Hcurl,V_Hdiv) + Π_div = interpolation_operator(u->u,V_H1_vec,V_Hdiv) + Π_curl = interpolation_operator(u->u,V_H1_vec,V_Hcurl) return G, C, Π_div, Π_curl end -function interpolation_operator(op,U_in,V_out,trian; +function interpolation_operator(op,U_in,V_out; strat=SubAssembledRows(), Tm=SparseMatrixCSR{0,PetscScalar,PetscInt}, Tv=Vector{PetscScalar}) @@ -24,7 +24,7 @@ function interpolation_operator(op,U_in,V_out,trian; in_basis = get_fe_basis(U_in) cell_interp_mats = out_dofs(op(in_basis)) - local_contr = map(local_views(trian),local_views(out_dofs),cell_interp_mats) do trian, dofs, arr + local_contr = map(local_views(out_dofs),cell_interp_mats) do dofs, arr contr = DomainContribution() add_contribution!(contr,get_triangulation(dofs),arr) return contr @@ -126,7 +126,7 @@ end n = 10 D = 3 -order = 2 +order = 1 np = (1,1,1)#Tuple(fill(1,D)) ranks = with_mpi() do distribute distribute(LinearIndices((prod(np),))) @@ -138,25 +138,25 @@ model = CartesianDiscreteModel(ranks,np,domain,ncells) trian = Triangulation(model) reffe_H1_sc = ReferenceFE(lagrangian,Float64,order) -V_H1_sc = FESpace(model,reffe_H1_sc;dirichlet_tags="boundary") +V_H1_sc = FESpace(model,reffe_H1_sc)#;dirichlet_tags="boundary") U_H1_sc = TrialFESpace(V_H1_sc) reffe_H1 = ReferenceFE(lagrangian,VectorValue{D,Float64},order) -V_H1 = FESpace(model,reffe_H1;dirichlet_tags="boundary") +V_H1 = FESpace(model,reffe_H1)#;dirichlet_tags="boundary") U_H1 = TrialFESpace(V_H1) reffe_Hdiv = ReferenceFE(raviart_thomas,Float64,order-1) -V_Hdiv = FESpace(model,reffe_Hdiv;dirichlet_tags="boundary") +V_Hdiv = FESpace(model,reffe_Hdiv)#;dirichlet_tags="boundary") U_Hdiv = TrialFESpace(V_Hdiv) reffe_Hcurl = ReferenceFE(nedelec,Float64,order-1) -V_Hcurl = FESpace(model,reffe_Hcurl;dirichlet_tags="boundary") +V_Hcurl = FESpace(model,reffe_Hcurl)#;dirichlet_tags="boundary") U_Hcurl = TrialFESpace(V_Hcurl) ############################################################################## dΩ = Measure(trian,(order+1)*2) -G, C, Π_div, Π_curl = get_operators(V_H1_sc,V_H1,V_Hcurl,V_Hdiv,trian); +G, C, Π_div, Π_curl = get_operators(V_H1_sc,V_H1,V_Hcurl,V_Hdiv); u(x) = x[1]^3 + x[2]^3 u_h1 = interpolate(u,U_H1_sc) @@ -191,7 +191,7 @@ f(x) = (D==2) ? VectorValue(x[1],x[2]) : VectorValue(x[1],x[2],x[3]) a(u,v) = ∫(u⋅v + α⋅(∇⋅u)⋅(∇⋅v)) * dΩ l(v) = ∫(f⋅v) * dΩ -V = FESpace(model,reffe_Hdiv;dirichlet_tags="boundary") +V = FESpace(model,reffe_Hdiv)#;dirichlet_tags="boundary") U = TrialFESpace(V,sol) op = AffineFEOperator(a,l,V,U) A = get_matrix(op)