diff --git a/src/GridapSolvers.jl b/src/GridapSolvers.jl index a6020006..1135ca96 100644 --- a/src/GridapSolvers.jl +++ b/src/GridapSolvers.jl @@ -1,57 +1,58 @@ module GridapSolvers - include("SolverInterfaces/SolverInterfaces.jl") - include("MultilevelTools/MultilevelTools.jl") - include("BlockSolvers/BlockSolvers.jl") - include("PatchBasedSmoothers/PatchBasedSmoothers.jl") - include("LinearSolvers/LinearSolvers.jl") - include("NonlinearSolvers/NonlinearSolvers.jl") - - using GridapSolvers.SolverInterfaces - using GridapSolvers.MultilevelTools - using GridapSolvers.BlockSolvers - using GridapSolvers.LinearSolvers - using GridapSolvers.PatchBasedSmoothers - using GridapSolvers.NonlinearSolvers - - # MultilevelTools - export get_parts, generate_level_parts, generate_subparts - - export ModelHierarchy, CartesianModelHierarchy, P4estCartesianModelHierarchy - export num_levels, get_level, get_level_parts - export get_model, get_model_before_redist - - export FESpaceHierarchy - export get_fe_space, get_fe_space_before_redist - export compute_hierarchy_matrices - - export DistributedGridTransferOperator - export RestrictionOperator, ProlongationOperator - export setup_transfer_operators - - # BlockSolvers - export BlockDiagonalSolver - - # LinearSolvers - export JacobiLinearSolver - export RichardsonSmoother - export SymGaussSeidelSmoother - export GMGLinearSolver - export BlockDiagonalSmoother - - export ConjugateGradientSolver - export IS_GMRESSolver - export IS_MINRESSolver - export IS_SSORSolver - - export CGSolver - export MINRESSolver - export GMRESSolver - export FGMRESSolver - - # PatchBasedSmoothers - export PatchDecomposition - export PatchFESpace - export PatchBasedLinearSolver +include("SolverInterfaces/SolverInterfaces.jl") +include("MultilevelTools/MultilevelTools.jl") +include("BlockSolvers/BlockSolvers.jl") +include("PatchBasedSmoothers/PatchBasedSmoothers.jl") +include("LinearSolvers/LinearSolvers.jl") +include("NonlinearSolvers/NonlinearSolvers.jl") + +using GridapSolvers.SolverInterfaces +using GridapSolvers.MultilevelTools +using GridapSolvers.BlockSolvers +using GridapSolvers.LinearSolvers +using GridapSolvers.PatchBasedSmoothers +using GridapSolvers.NonlinearSolvers + +# MultilevelTools +export get_parts, generate_level_parts, generate_subparts + +export ModelHierarchy, CartesianModelHierarchy, P4estCartesianModelHierarchy +export num_levels, get_level, get_level_parts +export get_model, get_model_before_redist + +export FESpaceHierarchy +export get_fe_space, get_fe_space_before_redist +export compute_hierarchy_matrices + +export DistributedGridTransferOperator +export RestrictionOperator, ProlongationOperator +export setup_transfer_operators + +# BlockSolvers +export BlockDiagonalSolver + +# LinearSolvers +export JacobiLinearSolver +export RichardsonSmoother +export SymGaussSeidelSmoother +export GMGLinearSolver +export BlockDiagonalSmoother + +export ConjugateGradientSolver +export IS_GMRESSolver +export IS_MINRESSolver +export IS_SSORSolver + +export CGSolver +export MINRESSolver +export GMRESSolver +export FGMRESSolver + +# PatchBasedSmoothers +export PatchDecomposition +export PatchFESpace +export PatchBasedLinearSolver +export Closure end diff --git a/src/LinearSolvers/CallbackSolver.jl b/src/LinearSolvers/CallbackSolver.jl new file mode 100644 index 00000000..87eca8d0 --- /dev/null +++ b/src/LinearSolvers/CallbackSolver.jl @@ -0,0 +1,52 @@ + +struct CallbackSolver{A,B} <: Algebra.LinearSolver + solver :: A + callback :: B + + function CallbackSolver(solver::LinearSolver,callback::Function) + A = typeof(solver) + B = typeof(callback) + new{A,B}(solver,callback) + end +end + +struct CallbackSolverSS{A,B} <: Algebra.SymbolicSetup + solver :: A + ss :: B +end + +function Algebra.symbolic_setup(solver::CallbackSolver,mat::AbstractMatrix) + ss = Algebra.symbolic_setup(solver.solver,mat) + return CallbackSolverSS(solver,ss) +end + +struct CallbackSolverNS{A,B} <: Algebra.NumericalSetup + solver :: A + ns :: B +end + +function Algebra.numerical_setup(ss::CallbackSolverSS,mat::AbstractMatrix) + ns = Algebra.numerical_setup(ss.ss,mat) + return CallbackSolverNS(ss.solver,ns) +end + +function Algebra.numerical_setup(ss::CallbackSolverSS,mat::AbstractMatrix,x::AbstractVector) + ns = Algebra.numerical_setup(ss.ss,mat,x) + return CallbackSolverNS(ss.solver,ns) +end + +function Algebra.numerical_setup!(ns::CallbackSolverNS,mat::AbstractMatrix) + Algebra.numerical_setup!(ns.ns,mat) + return ns +end + +function Algebra.numerical_setup!(ns::CallbackSolverNS,mat::AbstractMatrix,x::AbstractVector) + Algebra.numerical_setup!(ns.ns,mat,x) + return ns +end + +function Algebra.solve!(x::AbstractVector,ns::CallbackSolverNS,b::AbstractVector) + solve!(x,ns.ns,b) + ns.solver.callback(x) + return x +end diff --git a/src/LinearSolvers/GMGLinearSolvers.jl b/src/LinearSolvers/GMGLinearSolvers.jl index cfdf9502..b564e48b 100644 --- a/src/LinearSolvers/GMGLinearSolvers.jl +++ b/src/LinearSolvers/GMGLinearSolvers.jl @@ -415,7 +415,9 @@ function gmg_work_vectors(smatrices::AbstractVector{<:AbstractMatrix}) return work_vectors end -function apply_GMG_level!(lev::Integer,xh::Union{PVector,Nothing},rh::Union{PVector,Nothing},ns::GMGNumericalSetup) +function apply_GMG_level!( + lev::Integer,xh::Union{<:AbstractVector,Nothing},rh::Union{<:AbstractVector,Nothing},ns::GMGNumericalSetup +) mh = ns.solver.mh parts = get_level_parts(mh,lev) if i_am_in(parts) diff --git a/src/LinearSolvers/LinearSolverFromSmoothers.jl b/src/LinearSolvers/LinearSolverFromSmoothers.jl new file mode 100644 index 00000000..1e6264db --- /dev/null +++ b/src/LinearSolvers/LinearSolverFromSmoothers.jl @@ -0,0 +1,50 @@ + +struct LinearSolverFromSmoother{A} <: Algebra.LinearSolver + smoother :: A +end + +struct LinearSolverFromSmootherSS{A,B} <: Algebra.SymbolicSetup + smoother :: A + smoother_ss :: B +end + +struct LinearSolverFromSmootherNS{A,B,C} <: Algebra.NumericalSetup + smoother :: A + smoother_ns :: B + caches :: C +end + +function Gridap.Algebra.symbolic_setup(solver::LinearSolverFromSmoother, mat::AbstractMatrix) + ss = symbolic_setup(solver.smoother,mat) + return LinearSolverFromSmootherSS(solver.smoother,ss) +end + +function Gridap.Algebra.numerical_setup(ss::LinearSolverFromSmootherSS, mat::AbstractMatrix) + ns = numerical_setup(ss.smoother_ss, mat) + caches = allocate_in_domain(mat) + return LinearSolverFromSmootherNS(ss.smoother,ns,caches) +end + +function Gridap.Algebra.numerical_setup(ss::LinearSolverFromSmootherSS, mat::AbstractMatrix, vec::AbstractVector) + ns = numerical_setup(ss.smoother_ss, mat, vec) + caches = allocate_in_domain(mat) + return LinearSolverFromSmootherNS(ss.smoother,ns,caches) +end + +function Gridap.Algebra.numerical_setup!(ns::LinearSolverFromSmootherNS, mat::AbstractMatrix) + numerical_setup!(ns.smoother_ns, mat) + return ns +end + +function Gridap.Algebra.numerical_setup!(ns::LinearSolverFromSmootherNS, mat::AbstractMatrix, vec::AbstractVector) + numerical_setup!(ns.smoother_ns, mat, vec) + return ns +end + +function Gridap.Algebra.solve!(x::AbstractVector,ns::LinearSolverFromSmootherNS, b::AbstractVector) + r = ns.caches + fill!(x,zero(eltype(x))) + copy!(r,b) + solve!(x,ns.smoother_ns,r) + return x +end diff --git a/src/LinearSolvers/LinearSolvers.jl b/src/LinearSolvers/LinearSolvers.jl index 6cf18672..ff806b34 100644 --- a/src/LinearSolvers/LinearSolvers.jl +++ b/src/LinearSolvers/LinearSolvers.jl @@ -18,12 +18,16 @@ using GridapSolvers.MultilevelTools using GridapSolvers.SolverInterfaces using GridapSolvers.PatchBasedSmoothers +export LinearSolverFromSmoother export JacobiLinearSolver export RichardsonSmoother export SymGaussSeidelSmoother export GMGLinearSolver export BlockDiagonalSmoother export SchurComplementSolver +export SchwarzLinearSolver + +export CallbackSolver # Wrappers for IterativeSolvers.jl export IS_ConjugateGradientSolver @@ -49,12 +53,18 @@ include("PETSc/ElasticitySolvers.jl") include("PETSc/HipmairXuSolvers.jl") include("IdentityLinearSolvers.jl") + +include("LinearSolverFromSmoothers.jl") include("JacobiLinearSolvers.jl") include("RichardsonSmoothers.jl") include("SymGaussSeidelSmoothers.jl") + include("GMGLinearSolvers.jl") include("IterativeLinearSolvers.jl") include("SchurComplementSolvers.jl") include("MatrixSolvers.jl") +include("SchwarzLinearSolvers.jl") + +include("CallbackSolver.jl") end \ No newline at end of file diff --git a/src/LinearSolvers/SchwarzLinearSolvers.jl b/src/LinearSolvers/SchwarzLinearSolvers.jl new file mode 100644 index 00000000..cf41701e --- /dev/null +++ b/src/LinearSolvers/SchwarzLinearSolvers.jl @@ -0,0 +1,49 @@ + +# TODO: +# - Implement the multiplicative case +# - Add support for weights/averaging when aggregating the additive case? + +struct SchwarzLinearSolver{T,S,A} <: Algebra.LinearSolver + local_solvers::A + function SchwarzLinearSolver( + solver::Union{S,AbstractVector{<:S}}; + type = :additive + ) where S <: Algebra.LinearSolver + @check type in (:additive,:multiplicative) + @notimplementedif type == :multiplicative # TODO + A = typeof(solver) + new{type,S,A}(solver) + end +end + +struct SchwarzSymbolicSetup{T,S,A,B} <: Algebra.SymbolicSetup + solver::SchwarzLinearSolver{T,S,A} + local_ss::B +end + +function Algebra.symbolic_setup(s::SchwarzLinearSolver,mat::AbstractMatrix) + # TODO: This is where we should compute the comm coloring for the multiplicative case + expand(s) = map(m -> s,partition(mat)) + expand(s::AbstractVector) = s + + local_solvers = expand(s.local_solvers) + local_ss = map(symbolic_setup,local_solvers,partition(mat)) + return SchwarzSymbolicSetup(s,local_ss) +end + +struct SchwarzNumericalSetup{T,S,A,B} <: Algebra.NumericalSetup + solver::SchwarzLinearSolver{T,S,A} + local_ns::B +end + +function Algebra.numerical_setup(ss::SchwarzSymbolicSetup,mat::PSparseMatrix) + local_ns = map(numerical_setup,ss.local_ss,partition(mat)) + return SchwarzNumericalSetup(ss.solver,local_ns) +end + +function Algebra.solve!(x::PVector,ns::SchwarzNumericalSetup{:additive},b::PVector) + map(solve!,partition(x),ns.local_ns,partition(b)) + assemble!(x) |> wait + consistent!(x) |> wait + return x +end diff --git a/src/MultilevelTools/FESpaceHierarchies.jl b/src/MultilevelTools/FESpaceHierarchies.jl index 9f4cd05e..2331af41 100644 --- a/src/MultilevelTools/FESpaceHierarchies.jl +++ b/src/MultilevelTools/FESpaceHierarchies.jl @@ -56,6 +56,16 @@ function _cell_conformity( return CellConformity(cell_reffe,conformity) end +function _cell_conformity( + model::DiscreteModel, + reffe::ReferenceFE; + conformity=nothing, kwargs... +) :: CellConformity + cell_reffe = CompressedArray([reffe],fill(one(Int8),num_cells(model))) + conformity = Conformity(reffe,conformity) + return CellConformity(cell_reffe,conformity) +end + function _cell_conformity( model::GridapDistributed.DistributedDiscreteModel,args...;kwargs... ) :: AbstractVector{<:CellConformity} @@ -134,6 +144,31 @@ function Gridap.MultiField.MultiFieldFESpace(spaces::Vector{<:HierarchicalArray} end end +# Constant FESpaces + +function FESpaces.ConstantFESpace(mh::ModelHierarchy) + map(mh) do mhl + ConstantFESpace(mhl) + end +end + +function FESpaces.ConstantFESpace(mh::MultilevelTools.ModelHierarchyLevel) + reffe = ReferenceFE(lagrangian,Float64,0) + if has_redistribution(mh) + cparts, _ = get_old_and_new_parts(mh.red_glue,Val(false)) + Vh = i_am_in(cparts) ? ConstantFESpace(get_model_before_redist(mh)) : nothing + Vh_red = ConstantFESpace(get_model(mh)) + cell_conformity = i_am_in(cparts) ? _cell_conformity(get_model_before_redist(mh),reffe) : nothing + cell_conformity_red = _cell_conformity(get_model(mh),reffe) + else + Vh = ConstantFESpace(get_model(mh)) + Vh_red = nothing + cell_conformity = _cell_conformity(get_model(mh),reffe) + cell_conformity_red = nothing + end + return FESpaceHierarchyLevel(mh.level,Vh,Vh_red,cell_conformity,cell_conformity_red,mh) +end + # Computing system matrices function compute_hierarchy_matrices( diff --git a/src/MultilevelTools/DistributedGridTransferOperators.jl b/src/MultilevelTools/GridTransferOperators.jl similarity index 91% rename from src/MultilevelTools/DistributedGridTransferOperators.jl rename to src/MultilevelTools/GridTransferOperators.jl index 7d1f0936..8473097f 100644 --- a/src/MultilevelTools/DistributedGridTransferOperators.jl +++ b/src/MultilevelTools/GridTransferOperators.jl @@ -59,17 +59,17 @@ function DistributedGridTransferOperator( return DistributedGridTransferOperator(op_type,redist,restriction_method,sh,cache) end -function _get_interpolation_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode::Symbol) +function _get_interpolation_cache(lev::Int,sh::FESpaceHierarchy,qdegree,mode::Symbol) cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) model_h = get_model_before_redist(sh,lev) Uh = get_fe_space_before_redist(sh,lev) - fv_h = pfill(0.0,partition(Uh.gids)) + fv_h = zero_free_values(Uh) dv_h = (mode == :solution) ? get_dirichlet_dof_values(Uh) : zero_dirichlet_values(Uh) UH = get_fe_space(sh,lev+1) - fv_H = pfill(0.0,partition(UH.gids)) + fv_H = zero_free_values(UH) dv_H = (mode == :solution) ? get_dirichlet_dof_values(UH) : zero_dirichlet_values(UH) cache_refine = model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H @@ -82,7 +82,7 @@ function _get_interpolation_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mod return cache_refine end -function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode::Symbol,solver) +function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree,mode::Symbol,solver) cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) @@ -126,7 +126,7 @@ function _get_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,mode:: return cache_refine end -function _get_dual_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree::Int,solver) +function _get_dual_projection_cache(lev::Int,sh::FESpaceHierarchy,qdegree,solver) cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) @@ -219,10 +219,16 @@ function setup_restriction_operators(sh::FESpaceHierarchy,qdegrees::AbstractArra end end +function update_transfer_operator!( + op::DistributedGridTransferOperator,x::AbstractVector +) + nothing +end + ### Applying the operators: # A) Prolongation, without redistribution -function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:prolongation},Val{false}},x::PVector) +function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:prolongation},Val{false}},x::AbstractVector) cache_refine, cache_redist = A.cache model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine @@ -235,7 +241,7 @@ function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:p end # B.1) Restriction, without redistribution, by interpolation -function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:interpolation}},x::PVector) +function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:interpolation}},x::AbstractVector) cache_refine, cache_redist = A.cache model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine @@ -248,7 +254,7 @@ function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:r end # B.2) Restriction, without redistribution, by projection -function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:projection}},x::PVector) +function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:projection}},x::AbstractVector) cache_refine, cache_redist = A.cache model_h, Uh, fv_h, dv_h, VH, AH_ns, lH, xH, bH, bH0, assem = cache_refine @@ -265,7 +271,7 @@ function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:r end # B.3) Restriction, without redistribution, by dof selection (only nodal dofs) -function LinearAlgebra.mul!(y::PVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dof_mask}},x::PVector) +function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dof_mask}},x::AbstractVector) cache_refine, cache_redist = A.cache model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine @@ -367,6 +373,19 @@ end ############################################################### +function LinearAlgebra.mul!(y::AbstractVector,A::DistributedGridTransferOperator{Val{:restriction},Val{false},Val{:dual_projection}},x::AbstractVector) + 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) + v = get_fe_basis(VH) + assemble_vector!(y,assem,collect_cell_vector(VH,∫(v⋅uh)*dΩhH)) + + 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 diff --git a/src/MultilevelTools/GridapFixes.jl b/src/MultilevelTools/GridapFixes.jl index df058db1..712f6300 100644 --- a/src/MultilevelTools/GridapFixes.jl +++ b/src/MultilevelTools/GridapFixes.jl @@ -1,60 +1,42 @@ -""" -function Base.map(::typeof(Gridap.Arrays.testitem), - a::Tuple{<:AbstractVector{<:AbstractVector{<:VectorValue}},<:AbstractVector{<:Gridap.Fields.LinearCombinationFieldVector}}) - a2=Gridap.Arrays.testitem(a[2]) - a1=Vector{eltype(eltype(a[1]))}(undef,size(a2,1)) - a1.=zero(Gridap.Arrays.testitem(a1)) - (a1,a2) +# Necessary when integrating the result from LocalProjectionMaps onto MultiFieldFEBasisComponents +# in parallel. In general, this would also get triggered when doing a change_domain operation +# between two different Triangulation views. +# It has to do with blocks and how `extend` and `pos_neg_data` are getting dispatched. +function Geometry._similar_empty(val::Fields.LinearCombinationFieldVector) + #Fields.VoidBasis(val,true) + values = zeros(eltype(val.values),size(val.values)) + Gridap.Fields.LinearCombinationFieldVector(values,val.fields) end -""" - -# MultiField/DistributedMultiField missing API +# This below is another attempt to fix the issue. +# I'm close, but doesn't work. """ -function Gridap.FESpaces.zero_dirichlet_values(f::MultiFieldFESpace) - map(zero_dirichlet_values,f.spaces) -end - -function Gridap.FESpaces.interpolate_everywhere!(objects,free_values::AbstractVector,dirichlet_values::Vector,fe::MultiFieldFESpace) - blocks = SingleFieldFEFunction[] - for (field, (U,object)) in enumerate(zip(fe.spaces,objects)) - free_values_i = restrict_to_field(fe,free_values,field) - dirichlet_values_i = dirichlet_values[field] - uhi = interpolate!(object, free_values_i, dirichlet_values_i, U) - push!(blocks,uhi) +function Geometry._similar_empty(val::ArrayBlock{<:AbstractArray{<:Field},N}) where N + T = typeof(Geometry._similar_empty(testitem(val.array))) + array = Array{T,N}(undef,size(val.array)) + touched = copy(val.touched) + a = ArrayBlock(array,touched) + for i in eachindex(a) + if a.touched[i] + a.array[i] = Geometry._similar_empty(val.array[i]) + end end - Gridap.MultiField.MultiFieldFEFunction(free_values,fe,blocks) + a end -function Gridap.FESpaces.interpolate!(objects::GridapDistributed.DistributedMultiFieldFEFunction,free_values::AbstractVector,fe::GridapDistributed.DistributedMultiFieldFESpace) - part_fe_fun = map(local_views(objects),partition(free_values),local_views(fe)) do objects,x,f - interpolate!(objects,x,f) - end - field_fe_fun = GridapDistributed.DistributedSingleFieldFEFunction[] - for i in 1:num_fields(fe) - free_values_i = Gridap.MultiField.restrict_to_field(fe,free_values,i) - fe_space_i = fe.field_fe_space[i] - fe_fun_i = FEFunction(fe_space_i,free_values_i) - push!(field_fe_fun,fe_fun_i) - end - GridapDistributed.DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) +function Gridap.Geometry.pos_neg_data( + ipos_to_val::AbstractArray{<:ArrayBlock{<:AbstractArray{<:Field}}},i_to_iposneg::PosNegPartition +) + nineg = length(i_to_iposneg.ineg_to_i) + void = Geometry._similar_empty(testitem(ipos_to_val)) + ineg_to_val = Fill(void,nineg) + _ipos_to_val = lazy_map(Broadcasting(Fields.VoidBasisMap(false)),ipos_to_val) + _ipos_to_val, ineg_to_val end -function Gridap.FESpaces.FEFunction( - f::GridapDistributed.DistributedMultiFieldFESpace,x::AbstractVector, - dirichlet_values::AbstractArray{<:AbstractVector},isconsistent=false - ) - free_values = GridapDistributed.change_ghost(x,f.gids;is_consistent=isconsistent,make_consistent=true) - part_fe_fun = map(FEFunction,f.part_fe_space,partition(free_values)) - field_fe_fun = GridapDistributed.DistributedSingleFieldFEFunction[] - for i in 1:num_fields(f) - free_values_i = Gridap.MultiField.restrict_to_field(f,free_values,i) - dirichlet_values_i = dirichlet_values[i] - fe_space_i = f.field_fe_space[i] - fe_fun_i = FEFunction(fe_space_i,free_values_i,dirichlet_values_i,true) - push!(field_fe_fun,fe_fun_i) - end - GridapDistributed.DistributedMultiFieldFEFunction(field_fe_fun,part_fe_fun,free_values) +function Arrays.lazy_map(k::Broadcasting{<:Fields.VoidBasisMap},a::Arrays.LazyArray{<:Fill{<:Fields.BlockMap}}) + args = map(ai -> lazy_map(k.f, ai), a.args) + lazy_map(a.maps.value,args) end -""" \ No newline at end of file +""" diff --git a/src/MultilevelTools/LocalProjectionMaps.jl b/src/MultilevelTools/LocalProjectionMaps.jl index 01d8ec8a..c56f5b2b 100644 --- a/src/MultilevelTools/LocalProjectionMaps.jl +++ b/src/MultilevelTools/LocalProjectionMaps.jl @@ -129,17 +129,17 @@ function _compute_local_projections( q = SingleFieldFEBasis(test_shapefuns,Ω,TestBasis(),ReferenceDomain()) op = k.op.op - lhs_data = get_array(∫(p⋅q)dΩ) + lhs_data = get_array(∫(q⋅p)dΩ) rhs_data = get_array(∫(q⋅op(u))dΩ) basis_data = CellData.get_data(q) return lazy_map(k,lhs_data,rhs_data,basis_data) end -# Note on the caches: -# - We CANNOT overwrite `lhs`: In the case of constant cell_maps and `u` being a `FEFunction`, -# the lhs will be a Fill, but the rhs will not be optimized (regular LazyArray). -# In that case, we will see multiple times the same `lhs` being used for different `rhs`. -# - The converse never happens (I think), so we can overwrite `rhs` since it will always be recomputed. +function Arrays.return_value(::LocalProjectionMap,lhs::Matrix{T},rhs::A,basis) where {T,A<:Union{Matrix{T},Vector{T}}} + vec = zeros(T,size(rhs)) + return linear_combination(vec,basis) +end + function Arrays.return_cache(::LocalProjectionMap,lhs::Matrix{T},rhs::A,basis) where {T,A<:Union{Matrix{T},Vector{T}}} return CachedArray(copy(lhs)), CachedArray(copy(rhs)) end @@ -233,6 +233,11 @@ function _compute_local_projections( return lazy_map(k,lhs_data,rhs_data,basis_data,ids) end +function Arrays.return_value(::LocalProjectionMap,lhs::Matrix{T},rhs::A,basis,ids) where {T,A<:Union{Matrix{T},Vector{T}}} + vec = zeros(T,size(rhs)) + return linear_combination(vec,basis) +end + function Arrays.return_cache(::LocalProjectionMap,lhs::Matrix{T},rhs::A,basis,ids) where {T,A<:Union{Matrix{T},Vector{T}}} return CachedArray(copy(lhs)), CachedArray(copy(rhs)), CachedArray(zeros(eltype(rhs),(length(ids),size(rhs,2)))) end diff --git a/src/MultilevelTools/ModelHierarchies.jl b/src/MultilevelTools/ModelHierarchies.jl index 77654e46..ac66a789 100644 --- a/src/MultilevelTools/ModelHierarchies.jl +++ b/src/MultilevelTools/ModelHierarchies.jl @@ -139,6 +139,44 @@ function CartesianModelHierarchy( return HierarchicalArray(meshes,level_parts) end +function ModelHierarchy(models::Vector{<:DiscreteModel}) + nlevs = length(models) + @check all(map(i -> isa(models[i],AdaptedDiscreteModel),1:nlevs-1)) "Hierarchy models are not adapted models." + for lev in 1:nlevs-1 + @check Adaptivity.is_child(models[lev],models[lev+1]) "Incorrect hierarchy of models." + end + + level_parts = fill(DebugArray([1]),nlevs) + meshes = Vector{ModelHierarchyLevel}(undef,nlevs) + for lev in 1:nlevs-1 + glue = Gridap.Adaptivity.get_adaptivity_glue(models[lev]) + model = models[lev] + meshes[lev] = ModelHierarchyLevel(lev,model,glue,nothing,nothing) + end + meshes[nlevs] = ModelHierarchyLevel(nlevs,models[nlevs],nothing,nothing,nothing) + return HierarchicalArray(meshes,level_parts) +end + +function ModelHierarchy(models::Vector{<:GridapDistributed.DistributedDiscreteModel}) + nlevs = length(models) + @check all(map(i -> isa(models[i],GridapDistributed.DistributedAdaptedDiscreteModel),1:nlevs-1)) "Hierarchy models are not adapted models." + for lev in 1:nlevs-1 + @check Adaptivity.is_child(models[lev],models[lev+1]) "Incorrect hierarchy of models." + end + ranks = get_parts(models[1]) + @check all(m -> length(get_parts(m)) === length(ranks), models) "Models have different communicators." + + level_parts = fill(ranks,nlevs) + meshes = Vector{ModelHierarchyLevel}(undef,nlevs) + for lev in 1:nlevs-1 + glue = Gridap.Adaptivity.get_adaptivity_glue(models[lev]) + model = models[lev] + meshes[lev] = ModelHierarchyLevel(lev,model,glue,nothing,nothing) + end + meshes[nlevs] = ModelHierarchyLevel(nlevs,models[nlevs],nothing,nothing,nothing) + return HierarchicalArray(meshes,level_parts) +end + """ P4estCartesianModelHierarchy( ranks,np_per_level,domain,nc::NTuple{D,<:Integer}; diff --git a/src/MultilevelTools/MultiFieldTransferOperators.jl b/src/MultilevelTools/MultiFieldTransferOperators.jl index a1264a65..6b8c69c2 100644 --- a/src/MultilevelTools/MultiFieldTransferOperators.jl +++ b/src/MultilevelTools/MultiFieldTransferOperators.jl @@ -42,11 +42,14 @@ function MultiFieldTransferOperator(sh::FESpaceHierarchy,operators;op_type=:prol return mfops end -function update_transfer_operator!(op::MultiFieldTransferOperator,x::PVector) +function update_transfer_operator!( + op::MultiFieldTransferOperator{Val{:prolongation}},x::Union{AbstractVector,Nothing} +) xh, _ = op.cache - if !isnothing(xh) - copy!(x,xh) + if !isnothing(x) + copy!(xh,x) + isa(xh,PVector) && wait(consistent!(xh)) end for (i,op_i) in enumerate(op.ops) @@ -55,7 +58,27 @@ function update_transfer_operator!(op::MultiFieldTransferOperator,x::PVector) end end -function LinearAlgebra.mul!(x,op::MultiFieldTransferOperator,y) +function update_transfer_operator!( + op::MultiFieldTransferOperator{Val{:restriction}},y::Union{AbstractVector,Nothing} +) + _, yh = op.cache + + if !isnothing(y) + copy!(yh,y) + isa(yh,PVector) && wait(consistent!(yh)) + end + + for (i,op_i) in enumerate(op.ops) + yh_i = isnothing(yh) ? nothing : MultiField.restrict_to_field(op.Vh_in,yh,i) + update_transfer_operator!(op_i,yh_i) + end +end + +function LinearAlgebra.mul!( + x::Union{Nothing,<:AbstractVector}, + op::MultiFieldTransferOperator, + y::Union{Nothing,<:AbstractVector} +) xh, yh = op.cache if !isnothing(yh) @@ -70,7 +93,7 @@ function LinearAlgebra.mul!(x,op::MultiFieldTransferOperator,y) if !isnothing(xh) copy!(x,xh) - consistent!(x) |> fetch + isa(x,PVector) && wait(consistent!(x)) end return x diff --git a/src/MultilevelTools/MultilevelTools.jl b/src/MultilevelTools/MultilevelTools.jl index 5fbb59a5..4ef8c95a 100644 --- a/src/MultilevelTools/MultilevelTools.jl +++ b/src/MultilevelTools/MultilevelTools.jl @@ -51,7 +51,7 @@ include("RefinementTools.jl") include("ModelHierarchies.jl") include("FESpaceHierarchies.jl") include("LocalProjectionMaps.jl") -include("DistributedGridTransferOperators.jl") +include("GridTransferOperators.jl") include("MultiFieldTransferOperators.jl") end \ No newline at end of file diff --git a/src/NonlinearSolvers/ContinuationFEOperators.jl b/src/NonlinearSolvers/ContinuationFEOperators.jl new file mode 100644 index 00000000..4f14314f --- /dev/null +++ b/src/NonlinearSolvers/ContinuationFEOperators.jl @@ -0,0 +1,140 @@ + +""" + mutable struct ContinuationSwitch{A} + callback :: Function + caches :: A + switched :: Bool + end + + ContinuationSwitch(callback::Function, caches) + +Switch that implements the switching logic for a ContinuationFEOperator. + +The callback function must provide the following signatures: + +- Reset: `callback(u::FEFunction,::Nothing,cache) -> (switch::Bool, new_cache)` +- Update: `callback(u::FEFunction,b::AbstractVector,cache) -> (switch::Bool, new_cache)` + +where `u` is the current solution (as a FEFunction) and `b` is the current residual. The first +signature is called by `allocate_residual` and is typically used to reset the switch for a new +nonlinear solve. The second signature is called by `residual!` and is used to update the switch +based on the current residual. + +The cache is (potentially) mutated between each call, and can hold any extra information +needed for the continuation logic. +""" +mutable struct ContinuationSwitch{A} + callback :: Function + caches :: A + switched :: Bool + + function ContinuationSwitch(callback::Function, caches) + A = typeof(caches) + new{A}(callback, caches, false) + end +end + +has_switched(s::ContinuationSwitch) = s.switched + +switch!(s::ContinuationSwitch, u) = switch!(s, u, nothing) + +function switch!(s::ContinuationSwitch, u, b) + has_switched(s) && return s.switched + s.switched, s.caches = s.callback(u, b, s.caches) + if has_switched(s) + @debug "ContinuationFEOperator: Switching operators!" + end + return s.switched +end + +""" + ContinuationSwitch(niter::Int) + +Switch that will change operators after `niter` iterations. +""" +function ContinuationSwitch(niter::Int) + caches = (;it = -1, niter = niter) + callback(u,::Nothing,c) = (false, (;it = -1, niter = c.niter)) + callback(u,b,c) = (c.it+1 >= c.niter, (;it = c.it+1, niter = c.niter)) + return ContinuationSwitch(callback, caches) +end + +""" + struct ContinuationFEOperator <: FESpaces.FEOperator + op1 :: FEOperator + op2 :: FEOperator + switch :: ContinuationSwitch + reuse_caches :: Bool + end + +FEOperator implementing the continuation method for nonlinear solvers. It switches between +its two operators when the switch is triggered. + +Continuation between more that two operators can be achieved by daisy-chaining two or more +ContinuationFEOperators. + +If `reuse_caches` is `true`, the Jacobian of the first operator is reused for the second +operator. This is only possible if the sparsity pattern of the Jacobian does not change. +""" +struct ContinuationFEOperator{A,B,C} <: FESpaces.FEOperator + op1 :: A + op2 :: B + switch :: C + reuse_caches :: Bool + + function ContinuationFEOperator( + op1::FEOperator, + op2::FEOperator, + switch::ContinuationSwitch; + reuse_caches::Bool = true + ) + A, B, C = typeof(op1), typeof(op2), typeof(switch) + new{A,B,C}(op1, op2, switch, reuse_caches) + end +end + +has_switched(op::ContinuationFEOperator) = has_switched(op.switch) + +""" + ContinuationFEOperator(op1::FEOperator, op2::FEOperator, niter::Int; reuse_caches::Bool = true) + +ContinuationFEOperator that switches between `op1` and `op2` after `niter` iterations. +""" +function ContinuationFEOperator( + op1::FEOperator, op2::FEOperator, niter::Int; reuse_caches::Bool = true +) + switch = ContinuationSwitch(niter) + return ContinuationFEOperator(op1, op2, switch; reuse_caches) +end + +# FEOperator API + +function FESpaces.get_test(op::ContinuationFEOperator) + ifelse(!has_switched(op), get_test(op.op1), get_test(op.op2)) +end + +function FESpaces.get_trial(op::ContinuationFEOperator) + ifelse(!has_switched(op), get_trial(op.op1), get_trial(op.op2)) +end + +function Algebra.allocate_residual(op::ContinuationFEOperator,u) + switch!(op.switch, u, nothing) + ifelse(!has_switched(op), allocate_residual(op.op1, u), allocate_residual(op.op2, u)) +end + +function Algebra.residual!(b::AbstractVector,op::ContinuationFEOperator,u) + switch!(op.switch, u, b) + ifelse(!has_switched(op), residual!(b, op.op1, u), residual!(b, op.op2, u)) +end + +function Algebra.allocate_jacobian(op::ContinuationFEOperator,u) + ifelse(!has_switched(op), allocate_jacobian(op.op1, u), allocate_jacobian(op.op2, u)) +end + +function Algebra.jacobian!(A::AbstractMatrix,op::ContinuationFEOperator,u) + if op.reuse_caches + ifelse(!has_switched(op), jacobian!(A, op.op1, u), jacobian!(A, op.op2, u)) + else + ifelse(!has_switched(op), jacobian(op.op1, u), jacobian(op.op2, u)) + end +end diff --git a/src/NonlinearSolvers/NewtonRaphsonSolver.jl b/src/NonlinearSolvers/NewtonRaphsonSolver.jl index 93e0d32f..470366df 100644 --- a/src/NonlinearSolvers/NewtonRaphsonSolver.jl +++ b/src/NonlinearSolvers/NewtonRaphsonSolver.jl @@ -19,6 +19,8 @@ function NewtonSolver(ls;maxiter=100,atol=1e-12,rtol=1.e-6,verbose=0,name="Newto return NewtonSolver(ls,log) end +AbstractTrees.children(s::NewtonSolver) = [s.ls] + struct NewtonCache A::AbstractMatrix b::AbstractVector diff --git a/src/NonlinearSolvers/NonlinearSolvers.jl b/src/NonlinearSolvers/NonlinearSolvers.jl index 4eb24a53..8446ab63 100644 --- a/src/NonlinearSolvers/NonlinearSolvers.jl +++ b/src/NonlinearSolvers/NonlinearSolvers.jl @@ -1,4 +1,5 @@ module NonlinearSolvers + using AbstractTrees using LinearAlgebra using SparseArrays using SparseMatricesCSR @@ -12,7 +13,9 @@ module NonlinearSolvers using GridapSolvers.SolverInterfaces using GridapSolvers.MultilevelTools - using GridapSolvers.SolverInterfaces + + include("ContinuationFEOperators.jl") + export ContinuationFEOperator include("NewtonRaphsonSolver.jl") export NewtonSolver diff --git a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl b/src/PatchBasedSmoothers/DistributedPatchDecompositions.jl similarity index 71% rename from src/PatchBasedSmoothers/mpi/PatchDecompositions.jl rename to src/PatchBasedSmoothers/DistributedPatchDecompositions.jl index 7e0c8055..85897984 100644 --- a/src/PatchBasedSmoothers/mpi/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/DistributedPatchDecompositions.jl @@ -2,6 +2,14 @@ struct DistributedPatchDecomposition{Dr,Dc,Dp,A,B} <: GridapType patch_decompositions::A model::B + function DistributedPatchDecomposition( + patch_decompositions::AbstractVector{<:PatchDecomposition{Dr,Dc,Dp}}, + model::GridapDistributed.DistributedDiscreteModel{Dc,Dp} + ) where {Dr,Dc,Dp} + A = typeof(patch_decompositions) + B = typeof(model) + new{Dr,Dc,Dp,A,B}(patch_decompositions,model) + end end GridapDistributed.local_views(a::DistributedPatchDecomposition) = a.patch_decompositions @@ -12,17 +20,32 @@ function PatchDecomposition( patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude() ) where {Dc,Dp} mark_interface_facets!(model) - patch_decompositions = map(local_views(model)) do lmodel + patch_decompositions = map(local_views(model)) do model PatchDecomposition( - lmodel; + model; Dr=Dr, patch_boundary_style=patch_boundary_style, boundary_tag_names=["boundary","interface"] ) end - A = typeof(patch_decompositions) - B = typeof(model) - return DistributedPatchDecomposition{Dr,Dc,Dp,A,B}(patch_decompositions,model) + return DistributedPatchDecomposition(patch_decompositions,model) +end + +function CoarsePatchDecomposition( + model::GridapDistributed.DistributedAdaptedDiscreteModel{Dc,Dp}; + patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude() +) where {Dc,Dp} + gids = get_cell_gids(Gridap.Adaptivity.get_parent(model)) + mark_interface_facets!(model) + patch_decompositions = map(local_views(model),partition(gids)) do model, cids + own_cells = own_to_local(cids) + glue = Gridap.Adaptivity.get_adaptivity_glue(model) + patch_cells = glue.o2n_faces_map[own_cells] + PatchDecomposition( + model,patch_cells,patch_boundary_style,["boundary","interface"] + ) + end + return DistributedPatchDecomposition(patch_decompositions,model) end function PatchDecomposition(mh::ModelHierarchy;kwargs...) diff --git a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl b/src/PatchBasedSmoothers/DistributedPatchFESpaces.jl similarity index 97% rename from src/PatchBasedSmoothers/mpi/PatchFESpaces.jl rename to src/PatchBasedSmoothers/DistributedPatchFESpaces.jl index 01d3e736..a34150fa 100644 --- a/src/PatchBasedSmoothers/mpi/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/DistributedPatchFESpaces.jl @@ -1,4 +1,5 @@ + const DistributedPatchFESpace = GridapDistributed.DistributedSingleFieldFESpace{<:AbstractVector{<:PatchFESpace}} function PatchFESpace( @@ -44,8 +45,7 @@ function default_patches_mask(patch_decomposition::DistributedPatchDecomposition end function PatchFESpace( - sh::FESpaceHierarchy, - patch_decompositions::AbstractArray{<:DistributedPatchDecomposition} + sh::FESpaceHierarchy,patch_decompositions::AbstractArray ) nlevs = num_levels(sh) psh = map(view(sh,1:nlevs-1),patch_decompositions) do shl,decomp diff --git a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl b/src/PatchBasedSmoothers/PatchBasedLinearSolvers.jl similarity index 63% rename from src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl rename to src/PatchBasedSmoothers/PatchBasedLinearSolvers.jl index 07875b0b..b12801e4 100644 --- a/src/PatchBasedSmoothers/seq/PatchBasedLinearSolvers.jl +++ b/src/PatchBasedSmoothers/PatchBasedLinearSolvers.jl @@ -73,112 +73,92 @@ struct PatchBasedSmootherNumericalSetup{A,B,C,D} <: Gridap.Algebra.NumericalSetu caches :: D end -function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix) - @check !ss.solver.is_nonlinear - solver = ss.solver - Ph, Vh = solver.patch_space, solver.space - weights = solver.weighted ? compute_weight_operators(Ph,Vh) : nothing - - ap(u,v) = solver.biform(u,v) +function assemble_patch_matrices(Ph::FESpace,ap;local_solver=LUSolver()) assem = SparseMatrixAssembler(Ph,Ph) Ap = assemble_matrix(ap,assem,Ph,Ph) - Ap_ns = numerical_setup(symbolic_setup(solver.local_solver,Ap),Ap) - - # Caches - rp = allocate_in_range(Ap); fill!(rp,0.0) - dxp = allocate_in_domain(Ap); fill!(dxp,0.0) - caches = (rp,dxp) + Ap_ns = numerical_setup(symbolic_setup(local_solver,Ap),Ap) + return Ap, Ap_ns +end - return PatchBasedSmootherNumericalSetup(solver,nothing,Ap_ns,weights,caches) +function assemble_patch_matrices(Ph::GridapDistributed.DistributedFESpace,ap;local_solver=LUSolver()) + u, v = get_trial_fe_basis(Ph), get_fe_basis(Ph) + matdata = collect_cell_matrix(Ph,Ph,ap(u,v)) + Ap, Ap_ns = map(local_views(Ph),matdata) do Ph, matdata + assem = SparseMatrixAssembler(Ph,Ph) + Ap = assemble_matrix(assem,matdata) + Ap_ns = numerical_setup(symbolic_setup(local_solver,Ap),Ap) + return Ap, Ap_ns + end |> PartitionedArrays.tuple_of_arrays + return Ap, Ap_ns end -function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix,x::AbstractVector) - @check ss.solver.is_nonlinear - solver = ss.solver - Ph, Vh = solver.patch_space, solver.space - weights = solver.weighted ? compute_weight_operators(Ph,Vh) : nothing - - u0 = FEFunction(Vh,x) - ap(u,v) = solver.biform(u0,u,v) +function update_patch_matrices!(Ap,Ap_ns,Ph::FESpace,ap) assem = SparseMatrixAssembler(Ph,Ph) - Ap = assemble_matrix(ap,assem,Ph,Ph) - Ap_ns = numerical_setup(symbolic_setup(solver.local_solver,Ap),Ap) + assemble_matrix!(Ap,assem,Ph,Ph,ap) + numerical_setup!(Ap_ns,Ap) +end - # Caches - rp = allocate_in_range(Ap); fill!(rp,0.0) - dxp = allocate_in_domain(Ap); fill!(dxp,0.0) - caches = (rp,dxp) +function update_patch_matrices!(Ap,Ap_ns,Ph::GridapDistributed.DistributedFESpace,ap) + u, v = get_trial_fe_basis(Ph), get_fe_basis(Ph) + matdata = collect_cell_matrix(Ph,Ph,ap(u,v)) + map(Ap, Ap_ns, local_views(Ph), matdata) do Ap, Ap_ns, Ph, matdata + assem = SparseMatrixAssembler(Ph,Ph) + assemble_matrix!(Ap,assem,matdata) + numerical_setup!(Ap_ns,Ap) + end +end - return PatchBasedSmootherNumericalSetup(solver,Ap,Ap_ns,weights,caches) +function allocate_patch_workvectors(Ph::FESpace,Vh::FESpace) + rp = zero_free_values(Ph) + dxp = zero_free_values(Ph) + return rp,dxp end -function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::PSparseMatrix) +function allocate_patch_workvectors(Ph::GridapDistributed.DistributedFESpace,Vh::GridapDistributed.DistributedFESpace) + rp = zero_free_values(Ph) + dxp = zero_free_values(Ph) + r = zero_free_values(Vh) + x = zero_free_values(Vh) + return rp,dxp,r,x +end + +function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix) @check !ss.solver.is_nonlinear solver = ss.solver + local_solver = solver.local_solver Ph, Vh = solver.patch_space, solver.space - weights = solver.weighted ? compute_weight_operators(Ph,Vh) : nothing - - # Patch system solver (only local systems need to be solved) - u, v = get_trial_fe_basis(Vh), get_fe_basis(Vh) - contr = solver.biform(u,v) - matdata = collect_cell_matrix(Ph,Ph,contr) - Ap, Ap_ns = map(local_views(Ph),matdata) do Ph, matdata - assem = SparseMatrixAssembler(Ph,Ph) - Ap = assemble_matrix(assem,matdata) - Ap_ns = numerical_setup(symbolic_setup(solver.local_solver,Ap),Ap) - return Ap, Ap_ns - end |> PartitionedArrays.tuple_of_arrays - - # Caches - rp = pfill(0.0,partition(Ph.gids)) - dxp = pfill(0.0,partition(Ph.gids)) - r = pfill(0.0,partition(Vh.gids)) - x = pfill(0.0,partition(Vh.gids)) - caches = (rp,dxp,r,x) + ap(u,v) = solver.biform(u,v) + Ap, Ap_ns = assemble_patch_matrices(Ph,ap;local_solver) + weights = solver.weighted ? compute_weight_operators(Ph,Vh) : nothing + caches = allocate_patch_workvectors(Ph,Vh) return PatchBasedSmootherNumericalSetup(solver,nothing,Ap_ns,weights,caches) end -function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::PSparseMatrix,x::PVector) +function Gridap.Algebra.numerical_setup(ss::PatchBasedSymbolicSetup,A::AbstractMatrix,x::AbstractVector) @check ss.solver.is_nonlinear solver = ss.solver + local_solver = solver.local_solver Ph, Vh = solver.patch_space, solver.space + + u0 = FEFunction(Vh,x) + ap(u,v) = solver.biform(u0,u,v) + Ap, Ap_ns = assemble_patch_matrices(Ph,ap;local_solver) weights = solver.weighted ? compute_weight_operators(Ph,Vh) : nothing - - # Patch system solver (only local systems need to be solved) - u0, u, v = FEFunction(Vh,x), get_trial_fe_basis(Vh), get_fe_basis(Vh) - contr = solver.biform(u0,u,v) - matdata = collect_cell_matrix(Ph,Ph,contr) - Ap, Ap_ns = map(local_views(Ph),matdata) do Ph, matdata - assem = SparseMatrixAssembler(Ph,Ph) - Ap = assemble_matrix(assem,matdata) - Ap_ns = numerical_setup(symbolic_setup(solver.local_solver,Ap),Ap) - return Ap, Ap_ns - end |> PartitionedArrays.tuple_of_arrays - - # Caches - rp = pfill(0.0,partition(Ph.gids)) - dxp = pfill(0.0,partition(Ph.gids)) - r = pfill(0.0,partition(Vh.gids)) - x = pfill(0.0,partition(Vh.gids)) - caches = (rp,dxp,r,x) + caches = allocate_patch_workvectors(Ph,Vh) return PatchBasedSmootherNumericalSetup(solver,Ap,Ap_ns,weights,caches) end -function Gridap.Algebra.numerical_setup!(ns::PatchBasedSmootherNumericalSetup,A::PSparseMatrix,x::PVector) +function Gridap.Algebra.numerical_setup!(ns::PatchBasedSmootherNumericalSetup,A::AbstractMatrix,x::AbstractVector) @check ns.solver.is_nonlinear solver = ns.solver Ph, Vh = solver.patch_space, solver.space Ap, Ap_ns = ns.local_A, ns.local_ns - u0, u, v = FEFunction(Vh,x), get_trial_fe_basis(Vh), get_fe_basis(Vh) - contr = solver.biform(u0,u,v) - matdata = collect_cell_matrix(Ph,Ph,contr) - map(Ap, Ap_ns, local_views(Ph), matdata) do Ap, Ap_ns, Ph, matdata - assem = SparseMatrixAssembler(Ph,Ph) - assemble_matrix!(Ap,assem,matdata) - numerical_setup!(Ap_ns,Ap) - end + u0 = FEFunction(Vh,x) + ap(u,v) = solver.biform(u0,u,v) + update_patch_matrices!(Ap,Ap_ns,Ph,ap) + return ns end function Gridap.Algebra.solve!(x::AbstractVector,ns::PatchBasedSmootherNumericalSetup,r::AbstractVector) diff --git a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl index d00bf981..43053cae 100644 --- a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl +++ b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl @@ -11,25 +11,29 @@ using GridapDistributed using GridapSolvers.MultilevelTools -export PatchDecomposition +export PatchDecomposition, Closure export PatchFESpace -export PatchBasedLinearSolver +export PatchBasedLinearSolver, VankaSolver export PatchProlongationOperator, PatchRestrictionOperator export setup_patch_prolongation_operators, setup_patch_restriction_operators # Geometry -include("seq/PatchDecompositions.jl") -include("mpi/PatchDecompositions.jl") -include("seq/PatchTriangulations.jl") +include("PatchDecompositions.jl") +include("DistributedPatchDecompositions.jl") +include("PatchTriangulations.jl") +include("PatchClosures.jl") # FESpaces -include("seq/PatchFESpaces.jl") -include("mpi/PatchFESpaces.jl") -include("seq/PatchMultiFieldFESpaces.jl") +include("PatchFESpaces.jl") +include("DistributedPatchFESpaces.jl") +include("ZeroMeanPatchFESpaces.jl") +include("PatchMultiFieldFESpaces.jl") # Solvers -include("seq/PatchBasedLinearSolvers.jl") -include("seq/PatchTransferOperators.jl") +include("PatchBasedLinearSolvers.jl") +include("PatchTransferOperators.jl") +include("VankaSolvers.jl") +include("PatchSolvers.jl") end \ No newline at end of file diff --git a/src/PatchBasedSmoothers/PatchClosures.jl b/src/PatchBasedSmoothers/PatchClosures.jl new file mode 100644 index 00000000..723c094d --- /dev/null +++ b/src/PatchBasedSmoothers/PatchClosures.jl @@ -0,0 +1,40 @@ + +struct PatchClosureTriangulation{Dc,Dp} <: Triangulation{Dc,Dp} + trian :: PatchTriangulation{Dc,Dp} +end + +function Closure(PD::PatchDecomposition) + patch_cells = generate_patch_closures(PD) + trian = Triangulation(PD.model) + ptrian = PatchTriangulation(trian,PD,patch_cells,nothing) + return PatchClosureTriangulation(ptrian) +end + +function Geometry.get_background_model(t::PatchClosureTriangulation) + get_background_model(t.trian) +end + +function Geometry.get_grid(t::PatchClosureTriangulation) + get_grid(t.trian) +end + +function Geometry.get_glue(t::PatchClosureTriangulation,::Val{d}) where d + get_glue(t.trian,Val(d)) +end + +function Geometry.get_facet_normal(trian::PatchClosureTriangulation) + get_facet_normal(trian.trian) +end + +function Geometry.is_change_possible(strian::PatchClosureTriangulation,ttrian::PatchClosureTriangulation) + return is_change_possible(strian.trian,ttrian.trian) +end + +function Geometry.is_change_possible(strian::PatchTriangulation,ttrian::PatchClosureTriangulation) + return is_change_possible(strian,ttrian.trian) +end + +function Geometry.move_contributions(scell_to_val::AbstractArray,strian::PatchClosureTriangulation) + vals, _ = move_contributions(scell_to_val,strian.trian) + return vals, strian +end diff --git a/src/PatchBasedSmoothers/seq/PatchDecompositions.jl b/src/PatchBasedSmoothers/PatchDecompositions.jl similarity index 56% rename from src/PatchBasedSmoothers/seq/PatchDecompositions.jl rename to src/PatchBasedSmoothers/PatchDecompositions.jl index 764577e0..8435fe5f 100644 --- a/src/PatchBasedSmoothers/seq/PatchDecompositions.jl +++ b/src/PatchBasedSmoothers/PatchDecompositions.jl @@ -23,21 +23,17 @@ of `Ω` such that `Ω = Σ_i Ω_i`. - `Dr::Integer` : Dimension of the patch root - `model::DiscreteModel{Dc,Dp}` : Underlying discrete model - `patch_cells::Table` : [patch][local cell] -> cell -- `patch_cells_overlapped::Table` : [patch][local cell] -> overlapped cell - `patch_cells_faces_on_boundary::Table` : [d][overlapped cell][local face] -> face is on patch boundary """ struct PatchDecomposition{Dr,Dc,Dp} <: GridapType model :: DiscreteModel{Dc,Dp} - patch_cells :: Gridap.Arrays.Table # [patch][local cell] -> cell - patch_cells_overlapped :: Gridap.Arrays.Table # [patch][local cell] -> overlapped cell - patch_cells_faces_on_boundary :: Vector{Gridap.Arrays.Table} # [d][overlapped cell][local face] -> face is on patch boundary + patch_cells :: Arrays.Table # [patch][local cell] -> cell + patch_cells_faces_on_boundary :: Vector{Arrays.Table} # [d][overlapped cell][local face] -> face is on patch boundary + patch_boundary_style :: PatchBoundaryStyle end -num_patches(a::PatchDecomposition) = length(a.patch_cells) -Gridap.Geometry.num_cells(a::PatchDecomposition) = length(a.patch_cells.data) - -@doc """ +""" function PatchDecomposition( model::DiscreteModel{Dc,Dp}; Dr=0, @@ -49,31 +45,161 @@ Returns an instance of [`PatchDecomposition`](@ref) from a given discrete model. """ function PatchDecomposition( model::DiscreteModel{Dc,Dp}; - Dr=0, + Dr::Integer=0, patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude(), boundary_tag_names::AbstractArray{String}=["boundary"] ) where {Dc,Dp} - Gridap.Helpers.@check 0 <= Dr <= Dc-1 + @assert 0 <= Dr <= Dc-1 topology = get_grid_topology(model) - patch_cells = Gridap.Geometry.get_faces(topology,Dr,Dc) - patch_facets = Gridap.Geometry.get_faces(topology,Dr,Dc-1) - patch_cells_overlapped = compute_patch_overlapped_cells(patch_cells) + patch_cells = Geometry.get_faces(topology,Dr,Dc) + patch_facets = Geometry.get_faces(topology,Dr,Dc-1) patch_cells_faces_on_boundary = compute_patch_cells_faces_on_boundary( - model, patch_cells, patch_cells_overlapped, - patch_facets, patch_boundary_style, boundary_tag_names + model, patch_cells, patch_facets, patch_boundary_style, boundary_tag_names ) return PatchDecomposition{Dr,Dc,Dp}( - model, patch_cells, patch_cells_overlapped, patch_cells_faces_on_boundary + model, patch_cells, patch_cells_faces_on_boundary, patch_boundary_style ) end -function compute_patch_overlapped_cells(patch_cells) - num_overlapped_cells = length(patch_cells.data) - data = Gridap.Arrays.IdentityVector(num_overlapped_cells) - return Gridap.Arrays.Table(data,patch_cells.ptrs) +function PatchDecomposition( + model::DiscreteModel{Dc,Dp}, + patch_cells::AbstractArray{<:AbstractVector{<:Integer}}, + patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude(), + boundary_tag_names::AbstractArray{String}=["boundary"] +) where {Dc,Dp} + patch_facets = generate_patch_facets(model,patch_cells) + patch_cells_faces_on_boundary = compute_patch_cells_faces_on_boundary( + model, patch_cells, patch_facets, patch_boundary_style, boundary_tag_names + ) + return PatchDecomposition{Dc,Dc,Dp}( + model, patch_cells, patch_cells_faces_on_boundary, patch_boundary_style + ) +end + +function CoarsePatchDecomposition( + model::Gridap.Adaptivity.AdaptedDiscreteModel, + patch_boundary_style::PatchBoundaryStyle=PatchBoundaryExclude(), + boundary_tag_names::AbstractArray{String}=["boundary"] +) + glue = Gridap.Adaptivity.get_adaptivity_glue(model) + patch_cells = glue.o2n_faces_map + return PatchDecomposition(model,patch_cells,patch_boundary_style,boundary_tag_names) +end + +""" + num_patches(a::PatchDecomposition) +""" +num_patches(a::PatchDecomposition) = length(a.patch_cells) + +""" + get_patch_cells(PD::PatchDecomposition) -> patch_to_cells +""" +get_patch_cells(PD::PatchDecomposition) = PD.patch_cells + +""" + get_patch_cell_offsets(PD::PatchDecomposition) +""" +get_patch_cell_offsets(PD::PatchDecomposition) = PD.patch_cells.ptrs + +Geometry.num_cells(a::PatchDecomposition) = length(a.patch_cells.data) +Geometry.get_isboundary_face(PD::PatchDecomposition) = PD.patch_cells_faces_on_boundary +Geometry.get_isboundary_face(PD::PatchDecomposition,d::Integer) = PD.patch_cells_faces_on_boundary[d+1] + +""" + get_patch_cells_overlapped(PD::PatchDecomposition) -> patch_to_pcells +""" +function get_patch_cells_overlapped(PD::PatchDecomposition) + patch_cells = get_patch_cells(PD) + n_pcells = length(patch_cells.data) + data = Arrays.IdentityVector(n_pcells) + return Arrays.Table(data,patch_cells.ptrs) +end + +function Base.view(PD::PatchDecomposition{Dr,Dc,Dp},patch_ids) where {Dr,Dc,Dp} + patch_cells = view(get_patch_cells(PD),patch_ids) + patch_faces_on_boundary = [ + view(PD.patch_cells_faces_on_boundary[d+1],patch_ids) for d = 0:Dc-1 + ] + return PatchDecomposition{Dr,Dc,Dp}(PD.model,patch_cells,patch_faces_on_boundary,PD.patch_boundary_style) +end + +""" + patch_view(PD::PatchDecomposition,a::AbstractArray,patch::Integer) + patch_view(PD::PatchDecomposition,a::AbstractArray,patch_ids::AbstractUnitRange{<:Integer}) + +Returns a view of the pcell-wise array `a` restricted to the pcells of the patch `patch` or `patch_ids`. +""" +function patch_view(PD::PatchDecomposition,a::AbstractArray,patch::Integer) + offsets = get_patch_cell_offsets(PD) + return view(a,offsets[patch]:offsets[patch+1]-1) +end + +function patch_view(PD::PatchDecomposition,a::AbstractArray,patch_ids::AbstractUnitRange{<:Integer}) + offsets = get_patch_cell_offsets(PD) + start = offsets[first(patch_ids)] + stop = offsets[last(patch_ids)+1]-1 + return view(a,start:stop) +end + +""" + patch_reindex(PD::PatchDecomposition,cell_to_data) -> pcell_to_data +""" +function patch_reindex(PD::PatchDecomposition,cell_to_data) + patch_cells = get_patch_cells(PD) + pcell_to_data = lazy_map(Reindex(cell_to_data),patch_cells.data) + return pcell_to_data +end + +function patch_extend(PD::PatchDecomposition,patch_to_data) + pcell_to_patch = fill(0,num_cells(PD)) + patch_offsets = get_patch_cell_offsets(PD) + for patch in 1:num_patches(PD) + pcell_to_patch[patch_offsets[patch]:patch_offsets[patch+1]-1] .= patch + end + + pcell_to_data = lazy_map(Reindex(patch_to_data),pcell_to_patch) + return pcell_to_data +end + +""" + allocate_patch_cell_array(PD::PatchDecomposition,cell_to_data::Table{T};init=zero(T)) + +Allocates a patch-cell-wise array from a cell-wise array. +""" +function allocate_patch_cell_array( + PD::PatchDecomposition, cell_to_data::Table{T}; init=zero(T) +) where T + patch_cells = get_patch_cells(PD) + return allocate_patch_cell_array(patch_cells,cell_to_data;init) +end + +function allocate_patch_cell_array( + patch_cells::Table, cell_to_data::Table{T}; init = zero(T) +) where T + ptrs = zeros(Int,length(patch_cells.data)+1) + ptrs[1] = 1 + for (pcell,cell) in enumerate(patch_cells.data) + n = cell_to_data.ptrs[cell+1] - cell_to_data.ptrs[cell] + ptrs[pcell+1] = ptrs[pcell] + n + end + data = fill(init,ptrs[end]-1) + return Arrays.Table(data,ptrs) +end + +function allocate_patch_cell_array( + patch_cells::Table, cell_to_data::AbstractVector{<:AbstractVector{T}}; init = zero(T) +) where T + ptrs = zeros(Int,length(patch_cells.data)+1) + ptrs[1] = 1 + for (pcell,cell) in enumerate(patch_cells.data) + n = length(cell_to_data[cell]) + ptrs[pcell+1] = ptrs[pcell] + n + end + data = fill(init,ptrs[end]-1) + return Arrays.Table(data,ptrs) end # patch_cell_faces_on_boundary :: @@ -81,17 +207,16 @@ end function compute_patch_cells_faces_on_boundary( model::DiscreteModel, patch_cells, - patch_cells_overlapped, patch_facets, patch_boundary_style, boundary_tag_names ) patch_cell_faces_on_boundary = _allocate_patch_cells_faces_on_boundary(model,patch_cells) - if !isa(patch_boundary_style,PatchBoundaryInclude) + if isa(patch_boundary_style,PatchBoundaryExclude) _compute_patch_cells_faces_on_boundary!( patch_cell_faces_on_boundary, - model, patch_cells, patch_cells_overlapped, - patch_facets, patch_boundary_style, boundary_tag_names + model, patch_cells,patch_facets, + patch_boundary_style, boundary_tag_names ) end return patch_cell_faces_on_boundary @@ -103,7 +228,7 @@ function _allocate_patch_cells_faces_on_boundary(model::DiscreteModel{Dc},patch_ patch_cells_faces_on_boundary = Vector{Gridap.Arrays.Table}(undef,Dc) for d = 0:Dc-1 - ctype_to_num_dfaces = map(reffe -> num_faces(reffe,d),ctype_to_reffe) + ctype_to_num_dfaces = map(reffe -> num_faces(reffe,d), ctype_to_reffe) patch_cells_faces_on_boundary[d+1] = _allocate_ocell_to_dface(Bool, patch_cells,cell_to_ctype, ctype_to_num_dfaces) end @@ -111,11 +236,11 @@ function _allocate_patch_cells_faces_on_boundary(model::DiscreteModel{Dc},patch_ end function _allocate_ocell_to_dface(::Type{T},patch_cells,cell_to_ctype,ctype_to_num_dfaces) where T<:Number - num_overlapped_cells = length(patch_cells.data) - ptrs = Vector{Int}(undef,num_overlapped_cells+1) + n_pcells = length(patch_cells.data) + ptrs = Vector{Int}(undef,n_pcells+1) ptrs[1] = 1 - for i = 1:num_overlapped_cells + for i = 1:n_pcells cell = patch_cells.data[i] ctype = cell_to_ctype[cell] ptrs[i+1] = ptrs[i] + ctype_to_num_dfaces[ctype] @@ -128,7 +253,6 @@ function _compute_patch_cells_faces_on_boundary!( patch_cells_faces_on_boundary, model::DiscreteModel, patch_cells, - patch_cells_overlapped, patch_facets, patch_boundary_style, boundary_tag_names @@ -137,14 +261,14 @@ function _compute_patch_cells_faces_on_boundary!( cache_patch_cells = array_cache(patch_cells) cache_patch_facets = array_cache(patch_facets) for patch = 1:num_patches + first_patch_cell = patch_cells.ptrs[patch] current_patch_cells = getindex!(cache_patch_cells,patch_cells,patch) current_patch_facets = getindex!(cache_patch_facets,patch_facets,patch) _compute_patch_cells_faces_on_boundary!( patch_cells_faces_on_boundary, model, - patch, + first_patch_cell, current_patch_cells, - patch_cells_overlapped, current_patch_facets, patch_boundary_style, boundary_tag_names @@ -155,9 +279,8 @@ end function _compute_patch_cells_faces_on_boundary!( patch_cells_faces_on_boundary, model::DiscreteModel{Dc}, - patch, + first_patch_cell, patch_cells, - patch_cells_overlapped, patch_facets, patch_boundary_style, boundary_tag_names @@ -166,7 +289,7 @@ function _compute_patch_cells_faces_on_boundary!( topology = get_grid_topology(model) boundary_tags = findall(x -> (x ∈ boundary_tag_names), face_labeling.tag_to_name) - Gridap.Helpers.@check !isempty(boundary_tags) + @assert !isempty(boundary_tags) boundary_entities = vcat(face_labeling.tag_to_entities[boundary_tags]...) # Cells facets @@ -182,7 +305,7 @@ function _compute_patch_cells_faces_on_boundary!( # Go over all cells in the current patch for (lcell,cell) in enumerate(patch_cells) - overlapped_cell = patch_cells_overlapped.data[patch_cells_overlapped.ptrs[patch]+lcell-1] + overlapped_cell = first_patch_cell+lcell-1 cell_facets = getindex!(cache_cell_to_facets,cell_to_facets,cell) # Go over the facets (i.e., faces of dim Dc-1) in the current cell for (lfacet,facet) in enumerate(cell_facets) @@ -216,7 +339,7 @@ function _compute_patch_cells_faces_on_boundary!( lface = findfirst(x -> x==facet_face, cell_dfaces) lcell2 = findfirst(x -> x==cell_around_face, patch_cells) - overlapped_cell2 = patch_cells_overlapped.data[patch_cells_overlapped.ptrs[patch]+lcell2-1] + overlapped_cell2 = first_patch_cell+lcell2-1 position = patch_cells_faces_on_boundary[d+1].ptrs[overlapped_cell2]+lface-1 patch_cells_faces_on_boundary[d+1].data[position] = true end @@ -228,11 +351,19 @@ function _compute_patch_cells_faces_on_boundary!( end end -# Patch cell faces: -# patch_faces[pcell] = [face1, face2, ...] -# where face1, face2, ... are the faces of the overlapped cell `pcell` such that -# - they are NOT on the boundary of the patch -# - they are flagged `true` in faces_mask +""" + get_patch_cell_faces(PD::PatchDecomposition,Df::Integer) + get_patch_cell_faces(PD::PatchDecomposition,Df::Integer,faces_mask::AbstractVector{Bool}) + +Returns a patch-wise Table containing the faces on each patch cell, i.e + + patch_faces[pcell] = [face1, face2, ...] + +where face1, face2, ... are the faces on the overlapped cell `pcell` such that + + - they are NOT on the boundary of the patch + - they are flagged `true` in `faces_mask` +""" function get_patch_cell_faces(PD::PatchDecomposition,Df::Integer) model = PD.model topo = get_grid_topology(model) @@ -240,18 +371,19 @@ function get_patch_cell_faces(PD::PatchDecomposition,Df::Integer) return get_patch_cell_faces(PD,Df,faces_mask) end -function get_patch_cell_faces(PD::PatchDecomposition{Dr,Dc},Df::Integer,faces_mask) where {Dr,Dc} +function get_patch_cell_faces( + PD::PatchDecomposition{Dr,Dc},Df::Integer,faces_mask::AbstractVector{Bool} +) where {Dr,Dc} model = PD.model topo = get_grid_topology(model) c2e_map = Gridap.Geometry.get_faces(topo,Dc,Df) patch_cells = Gridap.Arrays.Table(PD.patch_cells) - patch_cell_faces = map(Reindex(c2e_map),patch_cells.data) + patch_cell_faces = lazy_map(Reindex(c2e_map),patch_cells.data) faces_on_boundary = PD.patch_cells_faces_on_boundary[Df+1] patch_faces = _allocate_patch_cell_faces(patch_cell_faces,faces_on_boundary,faces_mask) _generate_patch_cell_faces!(patch_faces,patch_cell_faces,faces_on_boundary,faces_mask) - return patch_faces end @@ -300,21 +432,29 @@ function _generate_patch_cell_faces!(patch_faces,patch_cell_faces,faces_on_bound return patch_faces end -# Patch faces: -# patch_faces[patch] = [face1, face2, ...] -# where face1, face2, ... are the faces of the patch such that -# - they are NOT on the boundary of the patch -# - they are flagged `true` in faces_mask -# If `reverse` is `true`, the faces are the ones ON the boundary of the patch. +""" + get_patch_faces(PD::PatchDecomposition,Df::Integer,faces_mask::AbstractVector{Bool};reverse=false) + +Returns a patch-wise Table containing the faces on each patch, i.e + + patch_faces[patch] = [face1, face2, ...] + +where face1, face2, ... are the faces on the patch such that + + - they are NOT on the boundary of the patch + - they are flagged `true` in `faces_mask` + +If `reverse` is `true`, the faces are the ones ON the boundary of the patch. +""" function get_patch_faces( - PD::PatchDecomposition{Dr,Dc},Df::Integer,faces_mask;reverse=false + PD::PatchDecomposition{Dr,Dc},Df::Integer,faces_mask::AbstractVector{Bool};reverse=false ) where {Dr,Dc} model = PD.model topo = get_grid_topology(model) c2e_map = Gridap.Geometry.get_faces(topo,Dc,Df) patch_cells = Gridap.Arrays.Table(PD.patch_cells) - patch_cell_faces = map(Reindex(c2e_map),patch_cells.data) + patch_cell_faces = lazy_map(Reindex(c2e_map),patch_cells.data) faces_on_boundary = PD.patch_cells_faces_on_boundary[Df+1] patch_faces = _allocate_patch_faces(patch_cells,patch_cell_faces,faces_on_boundary,faces_mask,reverse) @@ -392,12 +532,24 @@ function _generate_patch_faces!( return patch_faces end -# Face connectivity for the patches -# pface_to_pcell[pface] = [pcell1, pcell2, ...] -# This would be the Gridap equivalent to `get_faces(patch_topology,Df,Dc)`. -# On top of this, we also return the local cell index (within the face connectivity) of cell, -# which is needed when a pface touches different cells depending on the patch. -# The argument `patch_faces` allows to select only some pfaces (i.e boundary/skeleton/etc...). +""" + get_pface_to_pcell(PD::PatchDecomposition{Dr,Dc},Df::Integer,patch_faces) + +Returns two pface-wise Tables containing + + 1) the patch cells touched by each patch face and + 2) the local cell index (within the face connectivity) of the cell touched by the patch face, + which is needed when a pface touches different cells depending on the patch + +i.e + + pface_to_pcell[pface] = [pcell1, pcell2, ...] + pface_to_lcell[pface] = [lcell1, lcell2, ...] + +where pcell1, pcell2, ... are the patch cells touched by the patch face `pface`. + +This would be the Gridap equivalent to `get_faces(patch_topology,Df,Dc)`. +""" function get_pface_to_pcell(PD::PatchDecomposition{Dr,Dc},Df::Integer,patch_faces) where {Dr,Dc} model = PD.model topo = get_grid_topology(model) @@ -405,7 +557,7 @@ function get_pface_to_pcell(PD::PatchDecomposition{Dr,Dc},Df::Integer,patch_face faces_to_cells = Gridap.Geometry.get_faces(topo,Df,Dc) pfaces_to_cells = lazy_map(Reindex(faces_to_cells),patch_faces.data) patch_cells = Gridap.Arrays.Table(PD.patch_cells) - patch_cells_overlapped = PD.patch_cells_overlapped + patch_cells_overlapped = get_patch_cells_overlapped(PD) num_patches = length(patch_cells) num_pfaces = length(pfaces_to_cells) @@ -447,3 +599,109 @@ function get_pface_to_pcell(PD::PatchDecomposition{Dr,Dc},Df::Integer,patch_face pface_to_pcell = Gridap.Arrays.Table(data,ptrs) return pface_to_pcell, pface_to_lcell end + +""" + generate_patch_closures(PD::PatchDecomposition{Dr,Dc}) + +Returns a patch-wise Table containing the closure of each patch. +""" +function generate_patch_closures(PD::PatchDecomposition{Dr,Dc}) where {Dr,Dc} + topo = get_grid_topology(PD.model) + nodes_to_cells = Geometry.get_faces(topo,0,Dc) + cells_to_nodes = Geometry.get_faces(topo,Dc,0) + + patch_cells = get_patch_cells(PD) + + n_patches = length(patch_cells) + ptrs = zeros(Int,n_patches+1) + for patch in 1:n_patches + cells = view(patch_cells,patch) + closure = Set(cells) + for cell in cells + nodes = view(cells_to_nodes,cell) + for node in nodes + nbors = view(nodes_to_cells,node) + push!(closure,nbors...) + end + end + ptrs[patch+1] = length(closure) + end + Arrays.length_to_ptrs!(ptrs) + + data = zeros(Int,ptrs[end]-1) + for patch in 1:n_patches + cells = view(patch_cells,patch) + + # First we push the interior patch cells + for cell in cells + data[ptrs[patch]] = cell + ptrs[patch] += 1 + end + + # Then we push the extra cells in the closure + closure = Set(cells) + for cell in cells + nodes = view(cells_to_nodes,cell) + for node in nodes + nbors = view(nodes_to_cells,node) + for nbor in nbors + if nbor ∉ closure + data[ptrs[patch]] = nbor + push!(closure,nbor) + ptrs[patch] += 1 + end + end + end + end + end + Arrays.rewind_ptrs!(ptrs) + + return Table(data,ptrs) +end + +""" + generate_patch_facets(model::DiscreteModel{Dc,Dp},patch_cells) + +Given a model and the cells within each patch, returns a patch-wise Table containing the facets on each patch. +""" +function generate_patch_facets( + model::DiscreteModel{Dc,Dp},patch_cells +) where {Dc,Dp} + topo = get_grid_topology(model) + c2f_map = Geometry.get_faces(topo,Dc,Dc-1) + f2c_map = Geometry.get_faces(topo,Dc-1,Dc) + + touched = fill(false,num_faces(topo,Dc-1)) + ptrs = fill(0,length(patch_cells)+1) + for (patch,cells) in enumerate(patch_cells) + for c in cells + for f in c2f_map[c] + nbors = f2c_map[f] + if !touched[f] && (length(nbors) == 2) && all(n -> n ∈ cells, nbors) + touched[f] = true + ptrs[patch+1] += 1 + end + end + end + end + Arrays.length_to_ptrs!(ptrs) + + data = fill(0,ptrs[end]-1) + fill!(touched,false) + for (patch,cells) in enumerate(patch_cells) + for c in cells + for f in c2f_map[c] + nbors = f2c_map[f] + if !touched[f] && (length(nbors) == 2) && all(n -> n ∈ cells, nbors) + touched[f] = true + data[ptrs[patch]] = f + ptrs[patch] += 1 + end + end + end + end + Arrays.rewind_ptrs!(ptrs) + + patch_facets = Table(data,ptrs) + return patch_facets +end diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/PatchFESpaces.jl similarity index 63% rename from src/PatchBasedSmoothers/seq/PatchFESpaces.jl rename to src/PatchBasedSmoothers/PatchFESpaces.jl index d019307c..a95d992d 100644 --- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/PatchFESpaces.jl @@ -38,12 +38,24 @@ FESpace representing a patch-based subspace decomposition `V = Σ_i V_i` of a global space `V`. """ -struct PatchFESpace <: FESpaces.SingleFieldFESpace - Vh :: FESpaces.SingleFieldFESpace +struct PatchFESpace{A,B} <: FESpaces.SingleFieldFESpace + Vh :: A patch_decomposition :: PatchDecomposition num_dofs :: Int patch_cell_dofs_ids :: Arrays.Table dof_to_pdof :: Arrays.Table + metadata :: B + + function PatchFESpace( + space::SingleFieldFESpace, + patch_decomposition::PatchDecomposition, + num_dofs,patch_cell_dofs_ids,dof_to_pdof, + matadata=nothing + ) + A = typeof(space) + B = typeof(matadata) + new{A,B}(space,patch_decomposition,num_dofs,patch_cell_dofs_ids,dof_to_pdof,matadata) + end end @doc """ @@ -92,24 +104,24 @@ function PatchFESpace( patches_mask = Fill(false,num_patches(patch_decomposition)) ) cell_dofs_ids = get_cell_dof_ids(space) + patch_cells_overlapped = get_patch_cells_overlapped(patch_decomposition) patch_cell_dofs_ids, num_dofs = generate_patch_cell_dofs_ids( - get_grid_topology(patch_decomposition.model), patch_decomposition.patch_cells, - patch_decomposition.patch_cells_overlapped, + patch_cells_overlapped, patch_decomposition.patch_cells_faces_on_boundary, - cell_dofs_ids,cell_conformity,patches_mask + cell_dofs_ids,cell_conformity;patches_mask ) dof_to_pdof = generate_dof_to_pdof(space,patch_decomposition,patch_cell_dofs_ids) return PatchFESpace(space,patch_decomposition,num_dofs,patch_cell_dofs_ids,dof_to_pdof) end -FESpaces.get_dof_value_type(a::PatchFESpace) = Gridap.FESpaces.get_dof_value_type(a.Vh) -FESpaces.get_free_dof_ids(a::PatchFESpace) = Base.OneTo(a.num_dofs) -FESpaces.get_fe_basis(a::PatchFESpace) = get_fe_basis(a.Vh) -FESpaces.ConstraintStyle(::PatchFESpace) = Gridap.FESpaces.UnConstrained() -FESpaces.ConstraintStyle(::Type{PatchFESpace}) = Gridap.FESpaces.UnConstrained() -FESpaces.get_vector_type(a::PatchFESpace) = get_vector_type(a.Vh) -FESpaces.get_fe_dof_basis(a::PatchFESpace) = get_fe_dof_basis(a.Vh) +FESpaces.get_dof_value_type(a::PatchFESpace) = Gridap.FESpaces.get_dof_value_type(a.Vh) +FESpaces.get_free_dof_ids(a::PatchFESpace) = Base.OneTo(a.num_dofs) +FESpaces.get_fe_basis(a::PatchFESpace) = get_fe_basis(a.Vh) +FESpaces.ConstraintStyle(::PatchFESpace) = Gridap.FESpaces.UnConstrained() +FESpaces.ConstraintStyle(::Type{<:PatchFESpace}) = Gridap.FESpaces.UnConstrained() +FESpaces.get_vector_type(a::PatchFESpace) = get_vector_type(a.Vh) +FESpaces.get_fe_dof_basis(a::PatchFESpace) = get_fe_dof_basis(a.Vh) function Gridap.CellData.get_triangulation(a::PatchFESpace) PD = a.patch_decomposition @@ -136,14 +148,14 @@ function FESpaces.get_cell_dof_ids(::Triangulation,a::PatchFESpace,trian::PatchT end function FESpaces.get_cell_dof_ids(::BoundaryTriangulation,a::PatchFESpace,trian::PatchTriangulation) - cell_dof_ids = get_cell_dof_ids(a) + cell_dof_ids = get_cell_dof_ids(a) pface_to_pcell = trian.pface_to_pcell pcells = isempty(pface_to_pcell) ? Int[] : lazy_map(x->x[1],pface_to_pcell) return lazy_map(Reindex(cell_dof_ids),pcells) end function FESpaces.get_cell_dof_ids(::SkeletonTriangulation,a::PatchFESpace,trian::PatchTriangulation) - cell_dof_ids = get_cell_dof_ids(a) + cell_dof_ids = get_cell_dof_ids(a) pface_to_pcell = trian.pface_to_pcell pcells_plus = isempty(pface_to_pcell) ? Int[] : lazy_map(x->x[1],pface_to_pcell) @@ -154,6 +166,11 @@ function FESpaces.get_cell_dof_ids(::SkeletonTriangulation,a::PatchFESpace,trian return lazy_map(Fields.BlockMap(2,[1,2]),plus,minus) end +function FESpaces.get_cell_dof_ids(a::PatchFESpace,trian::PatchClosureTriangulation) + patch_cells = trian.trian.patch_faces + return propagate_patch_dof_ids(a,patch_cells) +end + # scatter dof values function FESpaces.scatter_free_and_dirichlet_values(f::PatchFESpace,free_values,dirichlet_values) @@ -164,98 +181,71 @@ end # Construction of the patch cell dofs ids function generate_patch_cell_dofs_ids( - topology, patch_cells, patch_cells_overlapped, patch_cells_faces_on_boundary, cell_dofs_ids, - cell_conformity, - patches_mask + cell_conformity::CellConformity; + patches_mask = Fill(false,length(patch_cells)), + numbering = :global ) - patch_cell_dofs_ids = allocate_patch_cell_dofs_ids(patch_cells,cell_dofs_ids) + @assert numbering in [:global,:local] + patch_cell_dofs_ids = allocate_patch_cell_array(patch_cells,cell_dofs_ids;init=Int32(-1)) num_dofs = generate_patch_cell_dofs_ids!( - patch_cell_dofs_ids,topology, + patch_cell_dofs_ids, patch_cells,patch_cells_overlapped, patch_cells_faces_on_boundary, - cell_dofs_ids,cell_conformity,patches_mask + cell_dofs_ids,cell_conformity; + patches_mask,numbering ) return patch_cell_dofs_ids, num_dofs end -function allocate_patch_cell_dofs_ids(patch_cells,cell_dofs_ids) - cache_cells = array_cache(patch_cells) - cache_cdofs = array_cache(cell_dofs_ids) - - num_overlapped_cells = length(patch_cells.data) - ptrs = Vector{Int}(undef,num_overlapped_cells+1) - ptrs[1] = 1; ncells = 1 - for patch = 1:length(patch_cells) - cells = getindex!(cache_cells,patch_cells,patch) - for cell in cells - current_cell_dof_ids = getindex!(cache_cdofs,cell_dofs_ids,cell) - ptrs[ncells+1] = ptrs[ncells]+length(current_cell_dof_ids) - ncells += 1 - end - end - - @check num_overlapped_cells+1 == ncells - data = Vector{Int}(undef,ptrs[end]-1) - return Gridap.Arrays.Table(data,ptrs) -end - function generate_patch_cell_dofs_ids!( patch_cell_dofs_ids, - topology, patch_cells, patch_cells_overlapped, patch_cells_faces_on_boundary, cell_dofs_ids, - cell_conformity, - patches_mask + cell_conformity::CellConformity; + patches_mask = Fill(false,length(patch_cells)), + numbering = :global ) - cache = array_cache(patch_cells) - num_patches = length(patch_cells) - current_dof = 1 - for patch = 1:num_patches - current_patch_cells = getindex!(cache,patch_cells,patch) - current_dof = generate_patch_cell_dofs_ids!( + @assert numbering in [:global,:local] + dof_offset = 1 + for patch = 1:length(patch_cells) + current_patch_cells = view(patch_cells,patch) + dof_offset = generate_patch_cell_dofs_ids!( patch_cell_dofs_ids, - topology, - patch, - current_patch_cells, - patch_cells_overlapped, + patch,current_patch_cells,patch_cells_overlapped, patch_cells_faces_on_boundary, - cell_dofs_ids, - cell_conformity; - free_dofs_offset=current_dof, + cell_dofs_ids,cell_conformity; + dof_offset=dof_offset, mask=patches_mask[patch] ) + if numbering == :local + dof_offset = 1 + end end - return current_dof-1 + return dof_offset-1 end -# TO-THINK/STRESS: -# 1. MultiFieldFESpace case? -# 2. FESpaces which are directly defined on physical space? We think this case is covered by -# the fact that we are using a CellConformity instance to rely on ownership info. -# free_dofs_offset : the ID from which we start to assign free DoF IDs upwards -# Note: we do not actually need to generate a global numbering for Dirichlet DoFs. We can -# tag all as them with -1, as we are always imposing homogenous Dirichlet boundary -# conditions, and thus there is no need to address the result of interpolating Dirichlet -# Data into the FE space. +# Note: We do not actually need to generate a global numbering for Dirichlet DoFs. We can +# tag all as them with -1, as we are always imposing homogenous Dirichlet boundary +# conditions, and thus there is no need to address the result of interpolating Dirichlet +# Data into the FE space. function generate_patch_cell_dofs_ids!( patch_cell_dofs_ids, - topology, patch::Integer, patch_cells::AbstractVector{<:Integer}, patch_cells_overlapped::Gridap.Arrays.Table, patch_cells_faces_on_boundary, global_space_cell_dofs_ids, - cell_conformity; - free_dofs_offset=1, + cell_conformity::CellConformity; + dof_offset=1, mask=false ) - o = patch_cells_overlapped.ptrs[patch] + o = patch_cells_overlapped.ptrs[patch] if mask for lpatch_cell = 1:length(patch_cells) cell_overlapped_mesh = patch_cells_overlapped.data[o+lpatch_cell-1] @@ -266,7 +256,6 @@ function generate_patch_cell_dofs_ids!( else g2l = Dict{Int,Int}() Dc = length(patch_cells_faces_on_boundary) - d_to_cell_to_dface = [Gridap.Geometry.get_faces(topology,Dc,d) for d in 0:Dc-1] # Loop over cells of the patch (local_cell_id_within_patch) for (lpatch_cell,patch_cell) in enumerate(patch_cells) @@ -276,38 +265,31 @@ function generate_patch_cell_dofs_ids!( current_patch_cell_dofs_ids = view(patch_cell_dofs_ids.data,s:e) ctype = cell_conformity.cell_ctype[patch_cell] - # 1) DoFs belonging to faces (Df < Dc) face_offset = 0 - for d = 0:Dc-1 - num_cell_faces = length(d_to_cell_to_dface[d+1][patch_cell]) - for lface in 1:num_cell_faces + for d = 0:Dc + n_dfaces = cell_conformity.d_ctype_num_dfaces[d+1][ctype] + for lface in 1:n_dfaces for ldof in cell_conformity.ctype_lface_own_ldofs[ctype][face_offset+lface] gdof = global_space_cell_dofs_ids[patch_cell][ldof] - face_in_patch_boundary = patch_cells_faces_on_boundary[d+1][cell_overlapped_mesh][lface] + face_in_patch_boundary = (d != Dc) && patch_cells_faces_on_boundary[d+1][cell_overlapped_mesh][lface] dof_is_dirichlet = (gdof < 0) if face_in_patch_boundary || dof_is_dirichlet current_patch_cell_dofs_ids[ldof] = -1 elseif gdof in keys(g2l) current_patch_cell_dofs_ids[ldof] = g2l[gdof] else - g2l[gdof] = free_dofs_offset - current_patch_cell_dofs_ids[ldof] = free_dofs_offset - free_dofs_offset += 1 + g2l[gdof] = dof_offset + current_patch_cell_dofs_ids[ldof] = dof_offset + dof_offset += 1 end end end face_offset += cell_conformity.d_ctype_num_dfaces[d+1][ctype] end - - # 2) Interior DoFs - for ldof in cell_conformity.ctype_lface_own_ldofs[ctype][face_offset+1] - current_patch_cell_dofs_ids[ldof] = free_dofs_offset - free_dofs_offset += 1 - end end end - return free_dofs_offset + return dof_offset end function generate_dof_to_pdof(Vh,PD,patch_cell_dofs_ids) @@ -379,9 +361,117 @@ function _generate_dof_to_pdof!(dof_to_pdof,Vh,PD,patch_cell_dofs_ids) end end +function generate_pdof_to_dof( + patch_decomposition::PatchDecomposition, + cell_dof_ids::Table{Ti}, + patch_cell_lids::Table{Ti} +) where Ti <: Integer + + n_patches = num_patches(patch_decomposition) + patch_cells = get_patch_cells(patch_decomposition) + + ptrs = fill(0,n_patches+1) + for patch in 1:n_patches + cell_lids = patch_view(patch_decomposition,patch_cell_lids,patch) + ptrs[patch+1] = maximum(map(maximum,cell_lids)) + end + PartitionedArrays.length_to_ptrs!(ptrs) + + data = fill(0,ptrs[end]-1) + for (patch,cells) in enumerate(patch_cells) + cell_lids = patch_view(patch_decomposition,patch_cell_lids,patch) + for (lcell,cell) in enumerate(cells) + dofs = view(cell_dof_ids,cell) + pdofs = view(cell_lids,lcell) + for (ldof,dof) in zip(pdofs,dofs) + if ldof > 0 + data[ptrs[patch]+ldof-1] = dof + end + end + end + end + + return Arrays.Table(data,ptrs) +end + +# TODO: Just For lagrange multipliers, fiz this better +function generate_pdof_to_dof( + patch_decomposition::PatchDecomposition, + cell_dof_ids::Fill, + patch_cell_lids::Table{Ti} +) where Ti <: Integer + ptrs = collect(1:num_patches(patch_decomposition)+1) + data = collect(Fill(1,num_patches(patch_decomposition))) + return Arrays.Table(data,ptrs) +end + +""" + propagate_patch_dof_ids(patch_space::PatchFESpace,new_patch_cells::Table) + +Propagates the DoF ids of the patch_space to a new set of patch cells given by +a patch-wise Table `new_patch_cells`. +""" +function propagate_patch_dof_ids( + patch_space::PatchFESpace, + new_patch_cells::Table +) + space = patch_space.Vh + patch_decomposition = patch_space.patch_decomposition + patch_cells = get_patch_cells(patch_decomposition) + patch_pcells = get_patch_cells_overlapped(patch_decomposition) + + cell_dof_ids = get_cell_dof_ids(space) + patch_cell_dof_ids = patch_space.patch_cell_dofs_ids + ext_dof_ids = allocate_patch_cell_array(new_patch_cells,cell_dof_ids;init=Int32(-1)) + + new_pcell = 1 + n_patches = length(patch_cells) + for patch in 1:n_patches + cells = view(patch_cells,patch) + pcells = view(patch_pcells,patch) + + # Create local dof to pdof maps + dof_to_pdof = Dict{Int,Int}() + for (cell,pcell) in zip(cells,pcells) + dofs = view(cell_dof_ids,cell) + pdofs = view(patch_cell_dof_ids,pcell) + for (dof,pdof) in zip(dofs,pdofs) + if pdof != -1 + dof_to_pdof[dof] = pdof + end + end + end + + # Propagate dofs to patch extensions + for new_cell in view(new_patch_cells,patch) + dofs = view(cell_dof_ids,new_cell) + pdofs = view(ext_dof_ids,new_pcell) + + pos = findfirst(c -> c == new_cell, cells) + if !isnothing(pos) # Cell is in the patch + pcell = pcells[pos] + pdofs .= view(patch_cell_dof_ids,pcell) + else # Cell is new (not in the patch) + for (ldof,dof) in enumerate(dofs) + if haskey(dof_to_pdof,dof) + pdofs[ldof] = dof_to_pdof[dof] + else + pdofs[ldof] = -1 + end + end + end + new_pcell += 1 + end + end + + return ext_dof_ids +end + # x \in PatchFESpace # y \in SingleFESpace -function prolongate!(x,Ph::PatchFESpace,y;dof_ids=LinearIndices(y)) +function prolongate!( + x::AbstractVector,Ph::PatchFESpace,y::AbstractVector;dof_ids=LinearIndices(y) +) dof_to_pdof = Ph.dof_to_pdof ptrs = dof_to_pdof.ptrs @@ -396,7 +486,9 @@ end # x \in SingleFESpace # y \in PatchFESpace -function inject!(x,Ph::PatchFESpace,y) +function inject!( + x::AbstractVector,Ph::PatchFESpace,y::AbstractVector +) dof_to_pdof = Ph.dof_to_pdof ptrs = dof_to_pdof.ptrs @@ -410,7 +502,9 @@ function inject!(x,Ph::PatchFESpace,y) end end -function inject!(x,Ph::PatchFESpace,y,w,w_sums) +function inject!( + x::AbstractVector,Ph::PatchFESpace,y::AbstractVector,w,w_sums +) dof_to_pdof = Ph.dof_to_pdof ptrs = dof_to_pdof.ptrs diff --git a/src/PatchBasedSmoothers/seq/PatchMultiFieldFESpaces.jl b/src/PatchBasedSmoothers/PatchMultiFieldFESpaces.jl similarity index 89% rename from src/PatchBasedSmoothers/seq/PatchMultiFieldFESpaces.jl rename to src/PatchBasedSmoothers/PatchMultiFieldFESpaces.jl index efc56486..3d4fe41d 100644 --- a/src/PatchBasedSmoothers/seq/PatchMultiFieldFESpaces.jl +++ b/src/PatchBasedSmoothers/PatchMultiFieldFESpaces.jl @@ -76,21 +76,24 @@ function inject!(x,Ph::MultiFieldFESpace,y) end end -function prolongate!(x::PVector, - Ph::GridapDistributed.DistributedMultiFieldFESpace, - y::PVector; - is_consistent::Bool=false) +function prolongate!( + x::PVector, + Ph::GridapDistributed.DistributedMultiFieldFESpace, + y::PVector; + is_consistent::Bool=false +) if !is_consistent consistent!(y) |> fetch end map(prolongate!,partition(x),local_views(Ph),partition(y)) end -function inject!(x::PVector, - Ph::GridapDistributed.DistributedMultiFieldFESpace, - y::PVector; - make_consistent::Bool=true) - +function inject!( + x::PVector, + Ph::GridapDistributed.DistributedMultiFieldFESpace, + y::PVector; + make_consistent::Bool=true +) map(partition(x),local_views(Ph),partition(y)) do x,Ph,y inject!(x,Ph,y) end diff --git a/src/PatchBasedSmoothers/PatchSolvers.jl b/src/PatchBasedSmoothers/PatchSolvers.jl new file mode 100644 index 00000000..39960241 --- /dev/null +++ b/src/PatchBasedSmoothers/PatchSolvers.jl @@ -0,0 +1,164 @@ + +struct PatchSolver{A,B,C,D,E} <: Algebra.LinearSolver + space::A + patch_decomposition::B + patch_ids::C + patch_cell_lids::D + biform::E +end + +function PatchSolver( + space::FESpace,patch_decomposition::PatchDecomposition,biform,reffe;conformity=nothing +) + cell_conformity = MultilevelTools._cell_conformity( + get_background_model(get_triangulation(space)),reffe;conformity + ) + return PatchSolver(space,patch_decomposition,biform,cell_conformity) +end + +function PatchSolver( + space::FESpace,patch_decomposition::PatchDecomposition,biform,cell_conformity::CellConformity +) + cell_dof_ids = get_cell_dof_ids(space) + patch_cells = get_patch_cells(patch_decomposition) + patch_cells_overlapped = get_patch_cells_overlapped(patch_decomposition) + + patch_cell_lids, _ = generate_patch_cell_dofs_ids( + patch_cells,patch_cells_overlapped, + patch_decomposition.patch_cells_faces_on_boundary, + cell_dof_ids,cell_conformity;numbering=:local + ) + patch_ids = generate_pdof_to_dof( + patch_decomposition,cell_dof_ids,patch_cell_lids + ) + return PatchSolver(space,patch_decomposition,patch_ids,patch_cell_lids,biform) +end + +function PatchSolver( + space::MultiFieldFESpace,patch_decomposition::PatchDecomposition,biform,reffes;kwargs... +) + solvers = map((Si,reffe)->PatchSolver(Si,patch_decomposition,biform,reffe;kwargs...),space,reffes) + field_to_patch_cell_lids = map(s -> s.patch_cell_lids ,solvers) + field_to_patch_ids = map(s -> s.patch_ids ,solvers) + field_to_ndofs = map(s -> num_free_dofs(s), space) + + patch_cell_lids, patch_ids = block_patch_ids(patch_decomposition,field_to_patch_cell_lids,field_to_patch_ids,field_to_ndofs) + return PatchSolver(space,patch_decomposition,patch_ids,patch_cell_lids,biform) +end + +function compute_patch_offsets( + patch_decomposition::PatchDecomposition, + field_to_patch_ids::Vector +) + nfields = length(field_to_patch_ids) + npatches = num_patches(patch_decomposition) + + offsets = zeros(Int,(nfields,npatches)) + for i in 1:(nfields-1) + offsets[i+1,:] .= offsets[i,:] .+ map(length,field_to_patch_ids[i]) + end + + return offsets +end + +function block_patch_ids( + patch_decomposition::PatchDecomposition, + field_to_patch_cell_lids::Vector, + field_to_patch_ids::Vector, + field_to_ndofs::Vector +) + offsets = compute_patch_offsets(patch_decomposition,field_to_patch_ids) + nfields = length(field_to_patch_ids) + npatches = num_patches(patch_decomposition) + + field_offsets = zeros(Int,nfields) + for i in 1:(nfields-1) + field_offsets[i+1] = field_offsets[i] + field_to_ndofs[i] + end + + block_cell_lids = Any[] + block_ids = Any[] + for i in 1:nfields + patch_cell_lids_i = field_to_patch_cell_lids[i] + patch_ids_i = field_to_patch_ids[i] + if i == 1 + push!(block_cell_lids,patch_cell_lids_i) + push!(block_ids,patch_ids_i) + else + offsets_i = Int32.(offsets[i,:]) + pcell_to_offsets = patch_extend(patch_decomposition,offsets_i) + patch_cell_lids_i_b = lazy_map(Broadcasting(Gridap.MultiField._sum_if_first_positive),patch_cell_lids_i,pcell_to_offsets) + patch_ids_i_b = lazy_map(Broadcasting(Gridap.MultiField._sum_if_first_positive),patch_ids_i,Fill(field_offsets[i],npatches)) + push!(block_cell_lids,patch_cell_lids_i_b) + push!(block_ids,patch_ids_i_b) + end + end + patch_cell_lids = lazy_map(BlockMap(nfields,collect(1:nfields)),block_cell_lids...) + patch_ids = map(vcat,block_ids...) + return patch_cell_lids, patch_ids +end + +struct PatchSS{A} <: Algebra.SymbolicSetup + solver::PatchSolver{A} +end + +Algebra.symbolic_setup(s::PatchSolver,mat::AbstractMatrix) = PatchSS(s) + +struct PatchNS{A,B} <: Algebra.NumericalSetup + solver::A + cache ::B +end + +function Algebra.numerical_setup(ss::PatchSS,mat::AbstractMatrix) + T = eltype(mat) + + space = ss.solver.space + biform = ss.solver.biform + cellmat, _ = move_contributions( # TODO: Generalise + get_array(biform(get_trial_fe_basis(space),get_fe_basis(space))), + Triangulation(ss.solver.patch_decomposition) + ) + cache = (CachedArray(zeros(T,1)), CachedArray(zeros(T,1,1)), cellmat) + return PatchNS(ss.solver,cache) +end + +function Algebra.solve!(x::AbstractVector,ns::PatchNS,b::AbstractVector) + PD = ns.solver.patch_decomposition + patch_ids, patch_cell_lids = ns.solver.patch_ids, ns.solver.patch_cell_lids + c_xk, c_Ak, cellmat = ns.cache + fill!(x,0.0) + + rows_cache = array_cache(patch_cell_lids) + cols_cache = array_cache(patch_cell_lids) + vals_cache = array_cache(cellmat) + + add! = AddEntriesMap(+) + add_cache = return_cache(add!,c_Ak.array,first(cellmat),first(patch_cell_lids),first(patch_cell_lids)) + caches = add_cache, vals_cache, rows_cache, cols_cache + + n_patches = length(patch_ids) + for patch in 1:n_patches + ids = patch_ids[patch] + + n = length(ids) + setsize!(c_xk,(n,)) + setsize!(c_Ak,(n,n)) + Ak = c_Ak.array + bk = c_xk.array + + patch_cellmat = patch_view(PD,cellmat,patch) + patch_rows = patch_view(PD,patch_cell_lids,patch) + patch_cols = patch_view(PD,patch_cell_lids,patch) + fill!(Ak,0.0) + FESpaces._numeric_loop_matrix!(Ak,caches,patch_cellmat,patch_rows,patch_cols) + + copyto!(bk,view(b,ids)) + f = lu!(Ak;check=false) + @check issuccess(f) "Factorization failed" + ldiv!(f,bk) + + x[ids] .+= bk + end + + return x +end diff --git a/src/PatchBasedSmoothers/seq/PatchTransferOperators.jl b/src/PatchBasedSmoothers/PatchTransferOperators.jl similarity index 85% rename from src/PatchBasedSmoothers/seq/PatchTransferOperators.jl rename to src/PatchBasedSmoothers/PatchTransferOperators.jl index 55e1a41b..36025ecc 100644 --- a/src/PatchBasedSmoothers/seq/PatchTransferOperators.jl +++ b/src/PatchBasedSmoothers/PatchTransferOperators.jl @@ -66,22 +66,12 @@ function _get_patch_cache(lev,sh,PD,lhs,rhs,is_nonlinear,cache_refine) cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) # Patch-based correction fespace - glue = sh[lev].mh_level.ref_glue - patches_mask = get_coarse_node_mask(model_h,glue) cell_conformity = MultilevelTools.get_cell_conformity_before_redist(sh,lev) - Ph = PatchFESpace(Uh,PD,cell_conformity;patches_mask) + Ph = PatchFESpace(Uh,PD,cell_conformity) # Solver caches - u, v = get_trial_fe_basis(Uh), get_fe_basis(Uh) - contr = is_nonlinear ? lhs(zero(Uh),u,v) : lhs(u,v) - matdata = collect_cell_matrix(Ph,Ph,contr) - Ap_ns, Ap = map(local_views(Ph),matdata) do Ph, matdata - assem = SparseMatrixAssembler(Ph,Ph) - Ap = assemble_matrix(assem,matdata) - Ap_ns = numerical_setup(symbolic_setup(LUSolver(),Ap),Ap) - return Ap_ns, Ap - end |> tuple_of_arrays - Ap = is_nonlinear ? Ap : nothing + ap(u,v) = is_nonlinear ? lhs(zero(Uh),u,v) : lhs(u,v) + Ap, Ap_ns = assemble_patch_matrices(Ph,ap) # TODO: This is using an LUSolver duh = zero(Uh) dxp, rp = zero_free_values(Ph), zero_free_values(Ph) @@ -107,18 +97,12 @@ function MultilevelTools.update_transfer_operator!(op::PatchProlongationOperator end if !isa(fv_h,Nothing) - u, v = get_trial_fe_basis(Uh), get_fe_basis(Uh) - contr = op.is_nonlinear ? op.lhs(FEFunction(Uh,fv_h),u,v) : op.lhs(u,v) - matdata = collect_cell_matrix(Ph,Ph,contr) - map(Ap_ns,Ap,local_views(Ph),matdata) do Ap_ns, Ap, Ph, matdata - assem = SparseMatrixAssembler(Ph,Ph) - assemble_matrix!(Ap,assem,matdata) - numerical_setup!(Ap_ns,Ap) - end + ap(u,v) = op.is_nonlinear ? op.lhs(FEFunction(Uh,fv_h),u,v) : op.lhs(u,v) + update_patch_matrices!(Ap,Ap_ns,Ph,ap) end end -function LinearAlgebra.mul!(y::PVector,A::PatchProlongationOperator{Val{false}},x::PVector) +function LinearAlgebra.mul!(y::AbstractVector,A::PatchProlongationOperator{Val{false}},x::AbstractVector) cache_refine, cache_patch, cache_redist = A.caches model_h, Uh, fv_h, dv_h, UH, fv_H, dv_H = cache_refine Ph, Ap_ns, Ap, duh, dxp, rp = cache_patch @@ -130,7 +114,11 @@ function LinearAlgebra.mul!(y::PVector,A::PatchProlongationOperator{Val{false}}, uh = FEFunction(Uh,fv_h,dv_h) assemble_vector!(v->A.rhs(uh,v),rp,Ph) - map(solve!,partition(dxp),Ap_ns,partition(rp)) + if isa(y,PVector) + map(solve!,partition(dxp),Ap_ns,partition(rp)) + else + solve!(dxp,Ap_ns,rp) + end inject!(dxh,Ph,dxp) fv_h .= fv_h .- dxh copy!(y,fv_h) @@ -167,11 +155,11 @@ end function setup_patch_prolongation_operators(sh,lhs,rhs,qdegrees;is_nonlinear=false) map(view(linear_indices(sh),1:num_levels(sh)-1)) do lev - qdegree = isa(qdegrees,Number) ? qdegrees : qdegrees[lev] + qdegree = isa(qdegrees,Vector) ? qdegrees[lev] : qdegrees cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) model = get_model_before_redist(sh,lev) - PD = PatchDecomposition(model) + PD = CoarsePatchDecomposition(model) Ω = Triangulation(PD) dΩ = Measure(Ω,qdegree) lhs_i = is_nonlinear ? (u,du,dv) -> lhs(u,du,dv,dΩ) : (u,v) -> lhs(u,v,dΩ) @@ -230,7 +218,6 @@ struct PatchRestrictionOperator{R,A,B} end function PatchRestrictionOperator(lev,sh,Ip,rhs,qdegree;solver=LUSolver()) - cache_refine = MultilevelTools._get_dual_projection_cache(lev,sh,qdegree,solver) cache_redist = MultilevelTools._get_redistribution_cache(lev,sh,:residual,:restriction,:dual_projection,cache_refine) cache_patch = Ip.caches[2] @@ -242,12 +229,13 @@ function PatchRestrictionOperator(lev,sh,Ip,rhs,qdegree;solver=LUSolver()) end function MultilevelTools.update_transfer_operator!(op::PatchRestrictionOperator,x::Union{PVector,Nothing}) + # Note: Update is done in the prolongation operator, with which we share the cache nothing end function setup_patch_restriction_operators(sh,patch_prolongations,rhs,qdegrees;kwargs...) map(view(linear_indices(sh),1:num_levels(sh)-1)) do lev - qdegree = isa(qdegrees,Number) ? qdegrees : qdegrees[lev] + qdegree = isa(qdegrees,Vector) ? qdegrees[lev] : qdegrees cparts = get_level_parts(sh,lev+1) if i_am_in(cparts) model = get_model_before_redist(sh,lev) @@ -262,7 +250,7 @@ function setup_patch_restriction_operators(sh,patch_prolongations,rhs,qdegrees;k end end -function LinearAlgebra.mul!(y::PVector,A::PatchRestrictionOperator{Val{false}},x::PVector) +function LinearAlgebra.mul!(y::AbstractVector,A::PatchRestrictionOperator{Val{false}},x::AbstractVector) cache_refine, cache_patch, _ = A.caches model_h, Uh, VH, Mh_ns, rh, uh, assem, dΩhH = cache_refine Ph, Ap_ns, Ap, duh, dxp, rp = cache_patch @@ -272,17 +260,19 @@ function LinearAlgebra.mul!(y::PVector,A::PatchRestrictionOperator{Val{false}},x copy!(fv_h,x) fill!(rp,0.0) prolongate!(rp,Ph,fv_h) - map(solve!,partition(dxp),Ap_ns,partition(rp)) + if isa(y,PVector) + map(solve!,partition(dxp),Ap_ns,partition(rp)) + else + solve!(dxp,Ap_ns,rp) + end inject!(dxh,Ph,dxp) - consistent!(dxh) |> fetch assemble_vector!(v->A.rhs(duh,v),rh,Uh) fv_h .= fv_h .- rh - consistent!(fv_h) |> fetch solve!(rh,Mh_ns,fv_h) copy!(fv_h,rh) - consistent!(fv_h) |> fetch + isa(y,PVector) && wait(consistent!(fv_h)) v = get_fe_basis(VH) assemble_vector!(y,assem,collect_cell_vector(VH,∫(v⋅uh)*dΩhH)) diff --git a/src/PatchBasedSmoothers/seq/PatchTriangulations.jl b/src/PatchBasedSmoothers/PatchTriangulations.jl similarity index 81% rename from src/PatchBasedSmoothers/seq/PatchTriangulations.jl rename to src/PatchBasedSmoothers/PatchTriangulations.jl index 5c5ab3e8..9b1f1be1 100644 --- a/src/PatchBasedSmoothers/seq/PatchTriangulations.jl +++ b/src/PatchBasedSmoothers/PatchTriangulations.jl @@ -5,7 +5,7 @@ Wrapper around a Triangulation, for patch-based assembly. """ -struct PatchTriangulation{Dc,Dp,A,B,C,D} <: Gridap.Geometry.Triangulation{Dc,Dp} +struct PatchTriangulation{Dc,Dp,A,B,C,D} <: Triangulation{Dc,Dp} trian :: A PD :: B patch_faces :: C @@ -45,18 +45,18 @@ end # For now, I am disabling changes from PatchTriangulations to other Triangulations. # Reason: When the tface_to_mface map is not injective (i.e when we have overlapping), -# the glue is not well defined. Gridap will throe an error when trying to create +# the glue is not well defined. Gridap will throw an error when trying to create # the inverse map mface_to_tface. # I believe this could technically be relaxed in the future, but for now I don't see a # scenario where we would need this. -function Geometry.is_change_possible(strian::PatchTriangulation,ttrian::Triangulation) - return strian === ttrian +function Geometry.is_change_possible(strian::PatchTriangulation,ttrian::PatchTriangulation) + return (strian === ttrian) || (strian.PD === ttrian.PD) end # Constructors function Geometry.Triangulation(PD::PatchDecomposition) - patch_cells = Gridap.Arrays.Table(PD.patch_cells) + patch_cells = PD.patch_cells trian = Triangulation(PD.model) return PatchTriangulation(trian,PD,patch_cells,nothing) end @@ -75,7 +75,6 @@ function Geometry.BoundaryTriangulation( pface_to_pcell, pface_to_lcell = get_pface_to_pcell(PD,Df,patch_faces) trian = OverlappingBoundaryTriangulation(model,patch_faces.data,pface_to_lcell) - return PatchTriangulation(trian,PD,patch_faces,pface_to_pcell) end @@ -88,11 +87,7 @@ function Geometry.SkeletonTriangulation(PD::PatchDecomposition{Dr,Dc}) where {Dr patch_faces = get_patch_faces(PD,Df,is_interior) pface_to_pcell, _ = get_pface_to_pcell(PD,Df,patch_faces) - nfaces = length(patch_faces.data) - plus = OverlappingBoundaryTriangulation(model,patch_faces.data,fill(Int8(1),nfaces)) - minus = OverlappingBoundaryTriangulation(model,patch_faces.data,fill(Int8(2),nfaces)) - trian = SkeletonTriangulation(plus,minus) - + trian = OverlappingSkeletonTriangulation(model,patch_faces.data) return PatchTriangulation(trian,PD,patch_faces,pface_to_pcell) end @@ -124,11 +119,9 @@ function Geometry.move_contributions( scell_to_val::AbstractArray, strian::PatchTriangulation ) - display(ndims(eltype(scell_to_val))) return scell_to_val, strian end - # Overlapping BoundaryTriangulation # # This is the situation: I do not see any show-stopper for us to have an overlapping @@ -140,11 +133,23 @@ end # mostly copied from Gridap/Geometry/BoundaryTriangulations.jl function OverlappingBoundaryTriangulation( - model::DiscreteModel, + model::Gridap.Adaptivity.AdaptedDiscreteModel, face_to_bgface::AbstractVector{<:Integer}, face_to_lcell::AbstractVector{<:Integer} ) + trian = OverlappingBoundaryTriangulation( + Gridap.Adaptivity.get_model(model), + face_to_bgface, + face_to_lcell, + ) + return Gridap.Adaptivity.AdaptedTriangulation(trian,model) +end +function OverlappingBoundaryTriangulation( + model::DiscreteModel, + face_to_bgface::AbstractVector{<:Integer}, + face_to_lcell::AbstractVector{<:Integer} +) D = num_cell_dims(model) topo = get_grid_topology(model) bgface_grid = Grid(ReferenceFE{D-1},model) @@ -195,7 +200,8 @@ function OverlappingFaceToCellGlue( face_to_ftype, cell_to_ctype, cell_to_lface_to_pindex, - ctype_to_lface_to_ftype) + ctype_to_lface_to_ftype + ) end function overlapped_find_local_index( @@ -217,3 +223,27 @@ function overlapped_find_local_index( end return c_to_lc end + + +# Overlapping SkeletonTriangulation + +function OverlappingSkeletonTriangulation( + model::Gridap.Adaptivity.AdaptedDiscreteModel, + face_to_bgface::AbstractVector{<:Integer}, +) + trian = OverlappingSkeletonTriangulation( + Gridap.Adaptivity.get_model(model), + face_to_bgface, + ) + return Gridap.Adaptivity.AdaptedTriangulation(trian,model) +end + +function OverlappingSkeletonTriangulation( + model::DiscreteModel, + face_to_bgface::AbstractVector{<:Integer}, +) + nfaces = length(face_to_bgface) + plus = OverlappingBoundaryTriangulation(model,face_to_bgface,fill(Int8(1),nfaces)) + minus = OverlappingBoundaryTriangulation(model,face_to_bgface,fill(Int8(2),nfaces)) + return SkeletonTriangulation(plus,minus) +end diff --git a/src/PatchBasedSmoothers/VankaSolvers.jl b/src/PatchBasedSmoothers/VankaSolvers.jl new file mode 100644 index 00000000..6e05fdcc --- /dev/null +++ b/src/PatchBasedSmoothers/VankaSolvers.jl @@ -0,0 +1,95 @@ + +struct VankaSolver{A} <: Algebra.LinearSolver + patch_ids::A + function VankaSolver( + patch_ids::Union{T,AbstractVector{<:T}} # Serial / Distributed + ) where T <: AbstractVector{<:AbstractVector{<:Integer}} + A = typeof(patch_ids) + new{A}(patch_ids) + end +end + +function VankaSolver(space::FESpace) + trian = get_triangulation(space) + ncells = num_cells(trian) + patch_cells = Table(1:ncells,1:ncells+1) + return VankaSolver(space,patch_cells) +end + +function VankaSolver(space::FESpace,patch_decomposition::PatchDecomposition) + patch_cells = patch_decomposition.patch_cells + return VankaSolver(space,patch_cells) +end + +function VankaSolver(space::FESpace,patch_cells::Table{<:Integer}) + collect_ids(ids::AbstractArray) = ids + collect_ids(ids::ArrayBlock) = vcat(ids.array...) + + cell_ids = get_cell_dof_ids(space) + patch_ids = map(patch_cells) do cells + ids = vcat([collect_ids(cell_ids[cell]) for cell in cells]...) + filter!(x->x>0,ids) + sort!(ids) + unique!(ids) + return ids + end + return VankaSolver(patch_ids) +end + +function VankaSolver(space::GridapDistributed.DistributedMultiFieldFESpace) + local_solvers = map(VankaSolver,local_views(space)) + return SchwarzLinearSolver(local_solvers) +end + +function VankaSolver( + space::GridapDistributed.DistributedMultiFieldFESpace, + patch_decomposition::DistributedPatchDecomposition +) + local_solvers = map(VankaSolver,local_views(space),local_views(patch_decomposition)) + return SchwarzLinearSolver(local_solvers) +end + +struct VankaSS{A} <: Algebra.SymbolicSetup + solver::VankaSolver{A} +end + +Algebra.symbolic_setup(s::VankaSolver,mat::AbstractMatrix) = VankaSS(s) + +struct VankaNS{A,B,C} <: Algebra.NumericalSetup + solver::VankaSolver{A} + matrix::B + cache ::C +end + +function Algebra.numerical_setup(ss::VankaSS,mat::AbstractMatrix) + T = eltype(mat) + cache = (CachedArray(zeros(T,1)), CachedArray(zeros(T,1,1))) + return VankaNS(ss.solver,mat,cache) +end + +function Algebra.solve!(x::AbstractVector,ns::VankaNS,b::AbstractVector) + A, patch_ids = ns.matrix, ns.solver.patch_ids + c_xk, c_Ak = ns.cache + fill!(x,0.0) + + n_patches = length(patch_ids) + for patch in 1:n_patches + ids = patch_ids[patch] + + n = length(ids) + setsize!(c_xk,(n,)) + setsize!(c_Ak,(n,n)) + Ak = c_Ak.array + bk = c_xk.array + + copyto!(Ak,view(A,ids,ids)) + copyto!(bk,view(b,ids)) + f = lu!(Ak,NoPivot();check=false) + @check issuccess(f) "Factorization failed" + ldiv!(f,bk) + + x[ids] .+= bk + end + + return x +end diff --git a/src/PatchBasedSmoothers/ZeroMeanPatchFESpaces.jl b/src/PatchBasedSmoothers/ZeroMeanPatchFESpaces.jl new file mode 100644 index 00000000..0785dd9a --- /dev/null +++ b/src/PatchBasedSmoothers/ZeroMeanPatchFESpaces.jl @@ -0,0 +1,143 @@ + +const ZeroMeanPatchFESpace{CA,S,B} = PatchFESpace{ZeroMeanFESpace{CA,S},B} + +function PatchFESpace( + space::FESpaces.ZeroMeanFESpace, + patch_decomposition::PatchDecomposition, + reffe::Union{ReferenceFE,Tuple{<:ReferenceFEs.ReferenceFEName,Any,Any}}; + conformity=nothing, + patches_mask = Fill(false,num_patches(patch_decomposition)) +) + pspace = PatchFESpace(space.space.space,patch_decomposition,reffe;conformity,patches_mask) + + n_patches = num_patches(patch_decomposition) + n_pdofs = pspace.num_dofs + patch_cells = patch_decomposition.patch_cells + pcell_to_pdofs = pspace.patch_cell_dofs_ids + dof_to_pdof = pspace.dof_to_pdof + + dof_to_dvol = space.vol_i + pdof_to_dvol = fill(0.0,n_pdofs) + for (dof,pdofs) in enumerate(dof_to_pdof) + pdof_to_dvol[pdofs] .= dof_to_dvol[dof] + end + + patch_vol = fill(0.0,n_patches) + pdof_to_new_pdof = fill(0,n_pdofs) + n_pdofs_free = 0 + n_pdofs_fixed = 1 # -1 reserved for homogeneous dirichlet on filtered patches + for patch in 1:n_patches + if patches_mask[patch] + continue + end + cell_s = patch_cells.ptrs[patch] + cell_e = patch_cells.ptrs[patch+1]-1 + pdof_s = pcell_to_pdofs.ptrs[cell_s] + pdof_e = pcell_to_pdofs.ptrs[cell_e+1]-1 + + # Patch volume + patch_vol[patch] = sum(pdof_to_dvol[pcell_to_pdofs.data[pdof_s:pdof_e]]) + + # First pdof per patch is fixed + pdof = pcell_to_pdofs.data[pdof_s] + n_pdofs_fixed += 1 + pdof_to_new_pdof[pdof] = -n_pdofs_fixed + pcell_to_pdofs.data[pdof_s] = -n_pdofs_fixed + + # The rest is free + for k in pdof_s+1:pdof_e + pdof = pcell_to_pdofs.data[k] + n_pdofs_free += 1 + pdof_to_new_pdof[pdof] = n_pdofs_free + pcell_to_pdofs.data[k] = n_pdofs_free + end + end + + dof_to_new_pdof = Table(collect(pdof_to_new_pdof[dof_to_pdof.data]), dof_to_pdof.ptrs) + + metadata = (patch_vol, pdof_to_dvol, n_pdofs, n_pdofs_free, n_pdofs_fixed, patches_mask) + return PatchFESpace(space,patch_decomposition,n_pdofs_free,pcell_to_pdofs,dof_to_new_pdof,metadata) +end + +# x \in PatchFESpace +# y \in SingleFESpace +function prolongate!(x,Ph::ZeroMeanPatchFESpace,y;dof_ids=LinearIndices(y)) + dof_to_pdof = Ph.dof_to_pdof + fixed_dof = Ph.Vh.space.dof_to_fix + + ptrs = dof_to_pdof.ptrs + data = dof_to_pdof.data + z = VectorWithEntryInserted(y,fixed_dof,zero(eltype(x))) + for dof in dof_ids + for k in ptrs[dof]:ptrs[dof+1]-1 + pdof = data[k] + if pdof > 0 # Is this correct? Should we correct the values in each patch? + x[pdof] = z[dof] + end + end + end +end + +# x \in SingleFESpace +# y \in PatchFESpace +function inject!(x,Ph::ZeroMeanPatchFESpace,y) + dof_to_pdof = Ph.dof_to_pdof + fixed_vals, _ = _compute_new_patch_fixedvals(y,Ph) + fixed_dof = Ph.Vh.space.dof_to_fix + + ptrs = dof_to_pdof.ptrs + data = dof_to_pdof.data + z = VectorWithEntryInserted(x,fixed_dof,zero(eltype(x))) + for dof in 1:length(dof_to_pdof) + z[dof] = 0.0 + for k in ptrs[dof]:ptrs[dof+1]-1 + pdof = data[k] + if pdof > 0 + z[dof] += y[pdof] + elseif pdof != -1 + z[dof] += fixed_vals[-pdof-1] + end + end + end + + x .-= z.value + return x +end + +function _compute_new_patch_fixedvals(x,Ph::ZeroMeanPatchFESpace) + patch_vol, pdof_to_dvol, n_pdofs, n_pdofs_free, n_pdofs_fixed, patches_mask = Ph.metadata + PD = Ph.patch_decomposition + patch_cells = PD.patch_cells + pcell_to_pdofs = Ph.patch_cell_dofs_ids + + n_patches = num_patches(PD) + k = 0 + fixed_vals = fill(0.0,n_pdofs_fixed) + fixed_pdofs = fill(0,n_pdofs_fixed) + for patch in 1:n_patches + if patches_mask[patch] + continue + end + pcell_s = patch_cells.ptrs[patch] + pcell_e = patch_cells.ptrs[patch+1]-1 + + pdof_s = pcell_to_pdofs.ptrs[pcell_s] + pdof_e = pcell_to_pdofs.ptrs[pcell_e+1]-1 + + fixed_pdof = -pcell_to_pdofs.data[pdof_s] + free_pdofs = pcell_to_pdofs.data[pdof_s+1]:pcell_to_pdofs.data[pdof_e] + pdofs = (pcell_to_pdofs.data[pdof_s+1]:pcell_to_pdofs.data[pdof_e]+1) .+ k + + fv = view(x,free_pdofs) + dv = [zero(eltype(fv))] + vol_i = view(pdof_to_dvol,pdofs) + c = FESpaces._compute_new_fixedval(fv,dv,vol_i,patch_vol[patch],1) + + k += 1 + x[free_pdofs] .+= c + fixed_vals[k] = c + fixed_pdofs[k] = fixed_pdof + end + + return fixed_vals, fixed_pdofs +end diff --git a/src/SolverInterfaces/ConvergenceLogs.jl b/src/SolverInterfaces/ConvergenceLogs.jl index 87603b5c..17839d9e 100644 --- a/src/SolverInterfaces/ConvergenceLogs.jl +++ b/src/SolverInterfaces/ConvergenceLogs.jl @@ -62,6 +62,26 @@ end @inline get_tabulation(log::ConvergenceLog) = get_tabulation(log,2) @inline get_tabulation(log::ConvergenceLog,n::Int) = repeat(' ', n + 2*log.depth) +""" + set_depth!(log::ConvergenceLog,depth::Int) + set_depth!(log::LinearSolver,depth::Int) + +Sets the tabulation depth of the convergence log `log` to `depth`. +""" +function set_depth!(log::ConvergenceLog,depth::Int) + log.depth = depth + return log +end + +function set_depth!(solver::Algebra.LinearSolver,depth::Int) + if hasproperty(solver,:log) + set_depth!(solver.log,depth) + end + map(children(solver)) do child + set_depth!(child,depth+2) + end +end + """ reset!(log::ConvergenceLog{T}) @@ -122,7 +142,7 @@ function finalize!(log::ConvergenceLog{T},r::T) where T msg = @sprintf("Iterations: %3i - Residuals: %.2e, %.2e ", log.num_iters, r, r_rel) println(t,msg) if log.verbose > SOLVER_VERBOSE_LOW - footer = " Exiting $(log.name) solver " + footer = " Exiting $(log.name) solver " println(t,rpad(string(repeat('-',15),footer),55,'-')) end end diff --git a/src/SolverInterfaces/SolverInfos.jl b/src/SolverInterfaces/SolverInfos.jl index 9ceae76b..49cb7b98 100644 --- a/src/SolverInterfaces/SolverInfos.jl +++ b/src/SolverInterfaces/SolverInfos.jl @@ -4,10 +4,15 @@ struct SolverInfo data :: Dict{Symbol, Any} end -SolverInfo(name::String) = SolverInfo(name,Dict{Symbol, Any}()) +SolverInfo(name::String;kwargs...) = SolverInfo(name,Dict{Symbol,Any}(kwargs)) function get_solver_info(solver::Gridap.Algebra.LinearSolver) - return SolverInfo(string(typeof(solver))) + info = SolverInfo(string(typeof(solver))) + if hasproperty(solver,:log) + add_tolerance_info!(info,solver.log) + add_convergence_info!(info,solver.log) + end + return info end function merge_info!(a::SolverInfo,b::SolverInfo;prefix=b.name) diff --git a/src/SolverInterfaces/SolverInterfaces.jl b/src/SolverInterfaces/SolverInterfaces.jl index 24563838..8b9aa724 100644 --- a/src/SolverInterfaces/SolverInterfaces.jl +++ b/src/SolverInterfaces/SolverInterfaces.jl @@ -14,7 +14,7 @@ include("SolverInfos.jl") export SolverVerboseLevel, SolverConvergenceFlag export SolverTolerances, get_solver_tolerances, set_solver_tolerances! -export ConvergenceLog, init!, update!, finalize!, reset!, print_message +export ConvergenceLog, init!, update!, finalize!, reset!, print_message, set_depth! export SolverInfo diff --git a/src/SolverInterfaces/SolverTolerances.jl b/src/SolverInterfaces/SolverTolerances.jl index dfe4aaf2..f3ac835d 100644 --- a/src/SolverInterfaces/SolverTolerances.jl +++ b/src/SolverInterfaces/SolverTolerances.jl @@ -114,7 +114,7 @@ end Returns `true` if the solver has finished, `false` otherwise. """ -function finished(tols::SolverTolerances,niter,e_a,e_r) :: Bool +@inline function finished(tols::SolverTolerances,niter,e_a,e_r) :: Bool return (niter >= tols.maxiter) || converged(tols,niter,e_a,e_r) end @@ -123,7 +123,7 @@ end Returns `true` if the solver has converged, `false` otherwise. """ -function converged(tols::SolverTolerances,niter,e_a,e_r) :: Bool +@inline function converged(tols::SolverTolerances,niter,e_a,e_r) :: Bool return (e_r < tols.rtol) || (e_a < tols.atol) end diff --git a/test/LinearSolvers/SchwarzSolversTests.jl b/test/LinearSolvers/SchwarzSolversTests.jl new file mode 100644 index 00000000..b2fd0605 --- /dev/null +++ b/test/LinearSolvers/SchwarzSolversTests.jl @@ -0,0 +1,51 @@ +module SchwarzSolversTests + +using Test +using MPI +using Gridap, Gridap.Algebra +using GridapDistributed +using PartitionedArrays + +using GridapSolvers +using GridapSolvers.LinearSolvers + +function main(distribute,np) + sol(x) = x[1] + x[2] + f(x) = -Δ(sol)(x) + + parts = distribute(LinearIndices((prod(np),))) + model = CartesianDiscreteModel(parts,np,(0,1,0,1),(8,8)) + + 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⋅sol)*dΩ + + op = AffineFEOperator(a,l,Uh,Vh) + A, b = get_matrix(op), get_vector(op) + + P = SchwarzLinearSolver(LUSolver()) + solver = CGSolver(P;rtol=1.0e-8,verbose=i_am_main(parts)) + ns = numerical_setup(symbolic_setup(solver,A),A) + x = allocate_in_domain(A); fill!(x,zero(eltype(x))) + 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 # module \ No newline at end of file diff --git a/test/LinearSolvers/SmoothersTests.jl b/test/LinearSolvers/SmoothersTests.jl index 530029f0..8261c7fd 100644 --- a/test/LinearSolvers/SmoothersTests.jl +++ b/test/LinearSolvers/SmoothersTests.jl @@ -5,7 +5,6 @@ using MPI using Gridap, Gridap.Algebra using GridapDistributed using PartitionedArrays -using IterativeSolvers using GridapSolvers using GridapSolvers.LinearSolvers @@ -19,7 +18,6 @@ function smoothers_driver(parts,model,P) 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) @@ -29,15 +27,10 @@ function smoothers_driver(parts,model,P) op = AffineFEOperator(a,l,Uh,Vh) A, b = get_matrix(op), get_vector(op) - ss = symbolic_setup(P,A) - ns = numerical_setup(ss,A) - - x = allocate_in_domain(A); ; fill!(x,zero(eltype(x))) - x, history = IterativeSolvers.cg!(x,A,b; - verbose=i_am_main(parts), - reltol=1.0e-8, - Pl=ns, - log=true) + solver = CGSolver(P;rtol=1.0e-8,verbose=i_am_main(parts)) + ns = numerical_setup(symbolic_setup(solver,A),A) + x = allocate_in_domain(A); fill!(x,zero(eltype(x))) + solve!(x,ns,b) u = interpolate(sol,Uh) uh = FEFunction(Uh,x) @@ -52,12 +45,13 @@ end function main_smoother_driver(parts,model,smoother) if smoother === :richardson - P = RichardsonSmoother(JacobiLinearSolver(),5,2.0/3.0) + S = RichardsonSmoother(JacobiLinearSolver(),5,2.0/3.0) elseif smoother === :sym_gauss_seidel - P = SymGaussSeidelSmoother(5) + S = SymGaussSeidelSmoother(5) else error("Unknown smoother") end + P = LinearSolverFromSmoother(S) smoothers_driver(parts,model,P) end diff --git a/test/LinearSolvers/mpi/SchwarzSolversTests.jl b/test/LinearSolvers/mpi/SchwarzSolversTests.jl new file mode 100644 index 00000000..12e2466a --- /dev/null +++ b/test/LinearSolvers/mpi/SchwarzSolversTests.jl @@ -0,0 +1,10 @@ +module SmoothersTestsMPI +using MPI, PartitionedArrays +include("../SchwarzSolversTests.jl") + +with_mpi() do distribute + SchwarzSolversTests.main(distribute,(2,2)) # 2D + SchwarzSolversTests.main(distribute,(2,2,1)) # 3D +end + +end \ No newline at end of file diff --git a/test/NonlinearSolvers/NonlinearSolversTests.jl b/test/NonlinearSolvers/NonlinearSolversTests.jl index 74360760..871851fa 100644 --- a/test/NonlinearSolvers/NonlinearSolversTests.jl +++ b/test/NonlinearSolvers/NonlinearSolversTests.jl @@ -13,7 +13,9 @@ using Gridap.Algebra using GridapSolvers using GridapSolvers.NonlinearSolvers, GridapSolvers.LinearSolvers -function main(ranks,model,solver) +const ModelTypes = Union{<:DiscreteModel,<:GridapDistributed.DistributedDiscreteModel} + +function main(ranks,model::ModelTypes,solver::Symbol) u_sol(x) = sum(x) order = 1 @@ -40,6 +42,9 @@ function main(ranks,model,solver) nls = NLsolveNonlinearSolver(ls; show_trace=true, method=:anderson, linesearch=BackTracking(), iterations=40, m=10) elseif solver == :newton nls = NewtonSolver(ls,maxiter=20,verbose=true) + elseif solver == :newton_continuation + op = ContinuationFEOperator(op,op,2;reuse_caches=true) + nls = NewtonSolver(ls,maxiter=20,verbose=true) else @error "Unknown solver" end @@ -52,15 +57,15 @@ function main(ranks,model,solver) end # Serial -function main(solver) +function main(solver::Symbol) model = CartesianDiscreteModel((0,1,0,1),(8,8)) ranks = DebugArray([1]) main(ranks,model,solver) end # Distributed -function main(distribute,solver) - ranks = distribute(LinearIndices((4,))) +function main(distribute,np,solver::Symbol) + ranks = distribute(LinearIndices((prod(np),))) model = CartesianDiscreteModel((0,1,0,1),(8,8)) main(ranks,model,solver) end diff --git a/test/NonlinearSolvers/mpi/NonlinearSolversTests.jl b/test/NonlinearSolvers/mpi/NonlinearSolversTests.jl new file mode 100644 index 00000000..24cf4fcf --- /dev/null +++ b/test/NonlinearSolvers/mpi/NonlinearSolversTests.jl @@ -0,0 +1,10 @@ +module NonlinearSolversTestsMPI +using MPI, PartitionedArrays +include("../NonlinearSolversTests.jl") + +with_mpi() do distribute + NonlinearSolversTests.main(distribute,4,:newton) + NonlinearSolversTests.main(distribute,4,:newton_continuation) +end + +end \ No newline at end of file diff --git a/test/NonlinearSolvers/mpi/runtests.jl b/test/NonlinearSolvers/mpi/runtests.jl new file mode 100644 index 00000000..b43e871a --- /dev/null +++ b/test/NonlinearSolvers/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) --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/NonlinearSolvers/seq/runtests.jl b/test/NonlinearSolvers/seq/runtests.jl index c9c1b0d1..7600c8a4 100644 --- a/test/NonlinearSolvers/seq/runtests.jl +++ b/test/NonlinearSolvers/seq/runtests.jl @@ -1,9 +1,10 @@ using Test +using PartitionedArrays include("../NonlinearSolversTests.jl") -@testset "NonlinearSolvers" begin +@testset "NonlinearSolvers - Serial" begin @testset "NLSolvers" begin NonlinearSolversTests.main(:nlsolvers_newton) NonlinearSolversTests.main(:nlsolvers_trust_region) @@ -11,5 +12,18 @@ include("../NonlinearSolversTests.jl") end @testset "Newton" begin NonlinearSolversTests.main(:newton) + NonlinearSolversTests.main(:newton_continuation) + end +end + +@testset "NonlinearSolvers - Sequential" begin + @testset "NLSolvers" begin + NonlinearSolversTests.main(DebugArray,4,:nlsolvers_newton) + NonlinearSolversTests.main(DebugArray,4,:nlsolvers_trust_region) + NonlinearSolversTests.main(DebugArray,4,:nlsolvers_anderson) + end + @testset "Newton" begin + NonlinearSolversTests.main(DebugArray,4,:newton) + NonlinearSolversTests.main(DebugArray,4,:newton_continuation) end end diff --git a/test/_dev/MacroFEs/block_stokes.jl b/test/_dev/MacroFEs/block_stokes.jl new file mode 100644 index 00000000..10dae068 --- /dev/null +++ b/test/_dev/MacroFEs/block_stokes.jl @@ -0,0 +1,81 @@ + +using Gridap +using Gridap.Adaptivity, Gridap.Geometry, Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Arrays +using Gridap.MultiField, Gridap.Algebra + +using GridapSolvers +using GridapSolvers.BlockSolvers, GridapSolvers.LinearSolvers + +using FillArrays + +# Parameters + +Dc = 3 +reftype = 1 + +u_sol(x) = (Dc == 2) ? VectorValue(x[1]^2*x[2], -x[1]*x[2]^2) : VectorValue(x[1]^2*x[2], -x[1]*x[2]^2,0.0) +p_sol(x) = (x[1] - 1.0/2.0) + +domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) +nc = (Dc == 2) ? (20,20) : (4,4,4) +model = simplexify(CartesianDiscreteModel(domain,nc)) + +min_order = (reftype == 1) ? Dc : Dc-1 +order = max(2,min_order) +poly = (Dc == 2) ? TRI : TET +rrule = (reftype == 1) ? Adaptivity.BarycentricRefinementRule(poly) : Adaptivity.PowellSabinRefinementRule(poly) +reffes = Fill(LagrangianRefFE(VectorValue{Dc,Float64},poly,order),Adaptivity.num_subcells(rrule)) +reffe_u = Adaptivity.MacroReferenceFE(rrule,reffes) +reffe_p = LagrangianRefFE(Float64,poly,order-1) + +qdegree = 2*order +quad = Quadrature(poly,Adaptivity.CompositeQuadrature(),rrule,qdegree) + +V = FESpace(model,reffe_u,dirichlet_tags=["boundary"]) +Q = FESpace(model,reffe_p,conformity=:L2,constraint=:zeromean) +U = TrialFESpace(V,u_sol) + +Ω = Triangulation(model) +dΩ = Measure(Ω,quad) +@assert abs(sum(∫(p_sol)dΩ)) < 1.e-15 +@assert abs(sum(∫(divergence(u_sol))dΩ)) < 1.e-15 + +mfs = BlockMultiFieldStyle() +X = MultiFieldFESpace([U,Q];style=mfs) +Y = MultiFieldFESpace([V,Q];style=mfs) + +α = 1.e2 +f(x) = -Δ(u_sol)(x) + ∇(p_sol)(x) +graddiv(u,v) = ∫(α*(∇⋅v)⋅(∇⋅u))dΩ +lap(u,v) = ∫(∇(u)⊙∇(v))dΩ +a((u,p),(v,q)) = lap(u,v) + ∫(divergence(u)*q - divergence(v)*p)dΩ + graddiv(u,v) +l((v,q)) = ∫(f⋅v)dΩ + +op = AffineFEOperator(a,l,X,Y) +A = get_matrix(op) +b = get_vector(op) + +solver_u = LUSolver() +solver_p = CGSolver(JacobiLinearSolver();maxiter=20,atol=1e-14,rtol=1.e-6) +solver_p.log.depth = 2 + +bblocks = [LinearSystemBlock() LinearSystemBlock(); + LinearSystemBlock() BiformBlock((p,q) -> ∫(-(1.0/α)*p*q)dΩ,Q,Q)] +coeffs = [1.0 1.0; + 0.0 1.0] +P = BlockTriangularSolver(bblocks,[solver_u,solver_p],coeffs,:upper) +solver = FGMRESSolver(20,P;atol=1e-10,rtol=1.e-12,verbose=true) +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = allocate_in_domain(A); fill!(x,0.0) +solve!(x,ns,b) +xh = FEFunction(X,x) + +uh, ph = xh +eh_u = uh - u_sol +eh_p = ph - p_sol +sum(∫(eh_u⋅eh_u)dΩ) +sum(∫(eh_p*eh_p)dΩ) + +norm(assemble_vector( q -> ∫(divergence(uh)*q)dΩ, Q)) +abs(sum(∫(divergence(uh))dΩ)) diff --git a/test/_dev/MacroFEs/riesz_GMG.jl b/test/_dev/MacroFEs/riesz_GMG.jl new file mode 100644 index 00000000..2a00e7f9 --- /dev/null +++ b/test/_dev/MacroFEs/riesz_GMG.jl @@ -0,0 +1,201 @@ +using Gridap +using GridapSolvers +using LinearAlgebra +using FillArrays + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays, Gridap.Adaptivity +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +function mock_model() + coords = [VectorValue(0.0,0.0),VectorValue(1.0,0.0),VectorValue(0.5,sqrt(3.0)/2.0)] + polys = [TRI] + cell_types = Int8[1] + conn = Table([[1,2,3]]) + reffes = map(x->LagrangianRefFE(Float64,x,1),polys) + + ref_grid = UnstructuredGrid(coords,conn,reffes,cell_types;has_affine_map=true) + model = UnstructuredDiscreteModel(ref_grid) + + labels = get_face_labeling(model) + labels.d_to_dface_to_entity[1] .= [1,2,3] + labels.d_to_dface_to_entity[2] .= [4,5,6] + labels.d_to_dface_to_entity[3] .= [7] + for i in 1:7 + push!(labels.tag_to_entities,Int32[i]) + push!(labels.tag_to_name,"tag_$(i)") + end + push!(labels.tag_to_entities,Int32[1,2,3,4,5,6]) + push!(labels.tag_to_name,"boundary") + push!(labels.tag_to_entities,Int32[7]) + push!(labels.tag_to_name,"interior") + + return model +end + +u_exact(x) = VectorValue(-x[1],x[2]) + +Dc = 2 +#cmodel = simplexify(CartesianDiscreteModel((0,1,0,1),(2,2))) +cmodel = refine(mock_model()).model + +labels = get_face_labeling(cmodel) +#add_tag_from_tags!(labels,"newman",[5]) +#add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:3)...,4,5]) +add_tag_from_tags!(labels,"newman",[6]) +model = refine(cmodel) + +order = 2 +rrule = Adaptivity.BarycentricRefinementRule(TRI) +reffes = Fill(LagrangianRefFE(VectorValue{2,Float64},TRI,order),Adaptivity.num_subcells(rrule)) +reffe_u = Adaptivity.MacroReferenceFE(rrule,reffes) + +VH = TestFESpace(cmodel,reffe_u;dirichlet_tags="dirichlet") +UH = TrialFESpace(VH,u_exact) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) + +qdegree = 2*order +quad = Quadrature(rrule,qdegree) +Ωh = Triangulation(model) +dΩh = Measure(Ωh,quad) +ΩH = Triangulation(cmodel) +dΩH = Measure(ΩH,quad) +dΩHh = Measure(ΩH,Ωh,quad) + +Γh = Boundary(model.model) +dΓh = Measure(Γh,4*order) +nh = get_normal_vector(Γh) + +ν = 1.e-2 +α = 10.0 +f(x) = -ν*Δ(u_exact)(x) +graddiv(u,v,dΩ) = ∫(α*(∇⋅u)⋅(∇⋅v))*dΩ +function a(u,v,dΩ) + c = ∫(ν*∇(u)⊙∇(v))*dΩ + if α > 0.0 + c += graddiv(u,v,dΩ) + end + return c +end +l(v,dΩ) = ∫(v⋅f)dΩ + +∇u_exact(x) = ∇(u_exact)(x) +ah(x,y) = a(x,y,dΩh) +aH(x,y) = a(x,y,dΩH) +lh(v) = l(v,dΩh) + ∫(ν*v⋅(∇u_exact⋅nh))dΓh + +op_h = AffineFEOperator(ah,lh,Uh,Vh) +Ah = get_matrix(op_h) +bh = get_vector(op_h) +xh_exact = Ah\bh + +AH = assemble_matrix(aH,UH,VH) +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Uh,Vh) + +uh_exact = FEFunction(Uh,xh_exact) +l2_error(uh_exact,u_exact,dΩh) +l2_norm(∇⋅uh_exact,dΩh) + +topo = get_grid_topology(model) +patches_mask = get_isboundary_face(topo,0) +PD = PatchDecomposition(model) +Ph = PatchFESpace(Vh,PD,reffe_u;conformity=H1Conformity())#,patches_mask) +Ωp = Triangulation(PD) +dΩp = Measure(Ωp,qdegree) +ap(x,y) = a(x,y,dΩp) +smoother = RichardsonSmoother(PatchBasedLinearSolver(ap,Ph,Vh),20,0.01) +smoother_ns = numerical_setup(symbolic_setup(smoother,Ah),Ah) + +function project_f2c(rh) + Qrh = Mhh\rh + uh = FEFunction(Vh,Qrh) + ll(v) = ∫(v⋅uh)*dΩHh + assemble_vector(ll,VH) +end +function interp_c2f(xH) + get_free_dof_values(interpolate(FEFunction(VH,xH),Vh)) +end + +PDi = PatchBasedSmoothers.CoarsePatchDecomposition(model) +Ih = PatchFESpace(Vh,PDi,reffe_u;conformity=H1Conformity()) + +Ωpi = Triangulation(PDi) +dΩpi = Measure(Ωpi,qdegree) +api(x,y) = a(x,y,dΩpi) +Ai = assemble_matrix(api,Ih,Ih) + +function P1(dxH) + interp_c2f(dxH) +end +function R1(rh) + project_f2c(rh) +end + +function P2(dxH) + dxh = interp_c2f(dxH) + uh = FEFunction(Vh,dxh) + r̃h = assemble_vector(v -> graddiv(uh,v,dΩpi),Ih) + dx̃ = Ai\r̃h + Pdxh = fill(0.0,length(dxh)) + PatchBasedSmoothers.inject!(Pdxh,Ih,dx̃) + y = dxh - Pdxh + return y +end +function R2(rh) + r̃h = zero_free_values(Ih) + PatchBasedSmoothers.prolongate!(r̃h,Ih,rh) + dr̃h = Ai\r̃h + dxh = zero_free_values(Vh) + PatchBasedSmoothers.inject!(dxh,Ih,dr̃h) + drh = assemble_vector(v -> graddiv(FEFunction(Vh,dxh),v,dΩh),Vh) + rH = project_f2c(rh - drh) + return rH +end + + +xh = randn(size(Ah,2)) +rh = bh - Ah*xh +niters = 5 + +iter = 0 +error0 = norm(rh) +error = error0 +e_rel = error/error0 +while iter < niters && e_rel > 1.0e-8 + println("Iter $iter:") + println(" > Initial: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + println(" > Pre-smoother: ", norm(rh)) + + rH = R2(rh) + qH = AH\rH + qh = P2(qH) + + rh = rh - Ah*qh + xh = xh + qh + println(" > Post-correction: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + + iter += 1 + error = norm(rh) + e_rel = error/error0 + println(" > Final: ",error, " - ", e_rel) +end + +uh = FEFunction(Xh,xh) +l2_error(uh,u_exact,dΩh) +l2_norm(∇⋅uh,dΩh) diff --git a/test/_dev/MacroFEs/stokes.jl b/test/_dev/MacroFEs/stokes.jl index ea5887b1..d575fc53 100644 --- a/test/_dev/MacroFEs/stokes.jl +++ b/test/_dev/MacroFEs/stokes.jl @@ -18,7 +18,6 @@ reffes = Fill(LagrangianRefFE(VectorValue{2,Float64},TRI,order),Adaptivity.num_s reffe_u = Adaptivity.MacroReferenceFE(rrule,reffes) reffe_p = LagrangianRefFE(Float64,TRI,order-1) -#reffe_p = Adaptivity.MacroReferenceFE(rrule,reffes;conformity=L2Conformity()) qdegree = 2*order quad = Quadrature(TRI,Adaptivity.CompositeQuadrature(),rrule,qdegree) diff --git a/test/_dev/MultilevelTools/LocalProjectionMaps.jl b/test/_dev/MultilevelTools/LocalProjectionMaps.jl index eb820ad0..8af39925 100644 --- a/test/_dev/MultilevelTools/LocalProjectionMaps.jl +++ b/test/_dev/MultilevelTools/LocalProjectionMaps.jl @@ -3,33 +3,46 @@ using Gridap using GridapDistributed, PartitionedArrays using GridapSolvers using GridapSolvers.MultilevelTools -using Gridap.Arrays +using Gridap.Arrays, Gridap.CellData, Gridap.FESpaces -np = (2,2) +np = (2,2,1) ranks = with_debug() do distribute distribute(LinearIndices((prod(np),))) end +function half_empty_trian(ranks,model) + cell_ids = get_cell_gids(model) + trians = map(ranks,local_views(model),partition(cell_ids)) do rank, model, ids + cell_mask = zeros(Bool, num_cells(model)) + if rank ∈ (3,4) + cell_mask[own_to_local(ids)] .= true + end + Triangulation(model,cell_mask) + end + GridapDistributed.DistributedTriangulation(trians,model) +end + layer(x) = sign(x)*abs(x)^(1/3) -cmap(x) = VectorValue(layer(x[1]),layer(x[2])) -model = CartesianDiscreteModel((0,-1,0,-1),(10,10),map=cmap) +cmap(x) = VectorValue(layer(x[1]),layer(x[2]),x[3]) +model = CartesianDiscreteModel(ranks,np,(0,1,0,1,0,1),(5,5,5),map=cmap) -Ω = Triangulation(model) +#Ω = Triangulation(model) +Ω = half_empty_trian(ranks,model) dΩ = Measure(Ω,2) order = 2 -qdegree = 2*(order+1) -reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +qdegree = 2*(order) +reffe_u = ReferenceFE(lagrangian,VectorValue{3,Float64},order) reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) -V = TestFESpace(model,reffe_u,dirichlet_tags="boundary"); -Q = TestFESpace(model,reffe_p;conformity=:L2,constraint=:zeromean) +V = TestFESpace(Ω,reffe_u,dirichlet_tags="boundary"); +Q = TestFESpace(Ω,reffe_p;conformity=:L2) mfs = Gridap.MultiField.BlockMultiFieldStyle() X = MultiFieldFESpace([V,Q];style=mfs) -#Π_Qh = LocalProjectionMap(divergence,reffe_p,qdegree) -Π_Qh = LocalProjectionMap(divergence,Q,qdegree) +Π_Qh = LocalProjectionMap(divergence,reffe_p,qdegree) +#Π_Qh = LocalProjectionMap(divergence,Q,qdegree) A = assemble_matrix((u,v) -> ∫(∇(v)⊙∇(u))dΩ, V, V) D = assemble_matrix((u,v) -> ∫(Π_Qh(u)⋅(∇⋅v))dΩ, V, V) @@ -44,3 +57,68 @@ G = Matrix(D) F-G norm(F-G) maximum(abs.(F-G)) + +u,p = get_trial_fe_basis(X) +v,q = get_fe_basis(X) + +Πu = Π_Qh(u) +∇v = ∇⋅(v) +Πu ⋅ ∇v + +∫(Π_Qh(u)⋅(∇⋅v))dΩ + +cf = Π_Qh(u)⋅(∇⋅v) + +i = 1 +integrate(cf.fields.items[i],dΩ.measures.items[i]) + +map(num_cells,local_views(Ω)) + +data = Πu.fields.items[1].cell_field.args[1].args[1] +testitem(data) + +length(u.fields.items[1].cell_basis) + +####################################################################### + +Ωf = Triangulation(model,tags="interior") +Vf = TestFESpace(Ωf,reffe_u); +Qf = TestFESpace(Ωf,reffe_p;conformity=:L2) +Xf = MultiFieldFESpace([Vf,Qf];style=mfs) + +Ωv = Triangulation(model,[1,2,3]) +dΩv = Measure(Ωv,2) + +u1 = get_trial_fe_basis(Vf) +cf1 = Π_Qh(u1) +∫(cf1)dΩv # OK + +u2 = get_trial_fe_basis(Xf)[1] +cf2 = Π_Qh(u2) +testitem(cf2.cell_field.args[1].args[1]) +∫(cf2)dΩv # Not OK + +cf2_bis = GenericCellField(cf2.cell_field.args[1].args[1],Ωf,ReferenceDomain()) +∫(cf2_bis)dΩv # OK + +u3 = get_trial_fe_basis(Xf)[2] +∫(u3)dΩv # OK + +u3_bis = ∇⋅(u3) +u3_bisbis = change_domain(u3_bis,Ωv,ReferenceDomain()) + +################################### + +eltype(cf2.cell_field.args[1].args[1]) +eltype(cf2.cell_field.args[1]) + +val = testitem(cf2.cell_field.args[1]) +Gridap.Geometry._similar_empty(val) + +val_copy = deepcopy(val) + +Gridap.Geometry._similar_empty(val.array[1]) + +val1 = val.array[1] +zs = 0 .* size(val1) +void = similar(val,eltype(val1),zs) diff --git a/test/_dev/MultilevelTools/SerialModelHierarchies.jl b/test/_dev/MultilevelTools/SerialModelHierarchies.jl new file mode 100644 index 00000000..3febf2a7 --- /dev/null +++ b/test/_dev/MultilevelTools/SerialModelHierarchies.jl @@ -0,0 +1,61 @@ +using Gridap +using GridapDistributed, PartitionedArrays +using GridapSolvers + +using Gridap.MultiField, Gridap.Adaptivity +using GridapSolvers.LinearSolvers, GridapSolvers.MultilevelTools, GridapSolvers.PatchBasedSmoothers + +function get_bilinear_form(mh_lev,biform,qdegree) + model = GridapSolvers.get_model(mh_lev) + Ω = Triangulation(model) + dΩ = Measure(Ω,qdegree) + return (u,v) -> biform(u,v,dΩ) +end + +cmodel = CartesianDiscreteModel((0,1,0,1),(10,10)) +fmodel = refine(cmodel,2) +mh = ModelHierarchy([fmodel,cmodel]) + +order = 1 +reffe = ReferenceFE(lagrangian,Float64,order) +sh = FESpace(mh,reffe;dirichlet_tags="boundary") + +qdegree = 2*(order+1) +biform(u,v,dΩ) = ∫(∇(u)⋅∇(v))dΩ + +biforms = map(mhl -> get_bilinear_form(mhl,biform,qdegree),mh) + +smoothers = [RichardsonSmoother(JacobiLinearSolver(),10,0.2)] +prolongations = setup_prolongation_operators( + sh,qdegree;mode=:residual +) +restrictions = setup_restriction_operators( + sh,qdegree;mode=:residual +) + +gmg = GMGLinearSolver( + mh,sh,sh,biforms, + prolongations,restrictions, + pre_smoothers=smoothers, + post_smoothers=smoothers, + coarsest_solver=LUSolver(), + maxiter=3,mode=:preconditioner,verbose=true +) + +solver = CGSolver(gmg;verbose=true) + +# Finest level +model = GridapSolvers.get_model(mh,1) +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) +V = get_fe_space(sh,1) +a(u,v) = biform(u,v,dΩ) +l(v) = ∫(v)dΩ + +op = AffineFEOperator(a,l,V,V) +A, b = get_matrix(op), get_vector(op) + +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = zeros(size(b)) +solve!(x,ns,b) diff --git a/test/_dev/PatchBased/BoundaryTrians.jl b/test/_dev/PatchBased/BoundaryTrians.jl index bcc6a484..2ebf4513 100644 --- a/test/_dev/PatchBased/BoundaryTrians.jl +++ b/test/_dev/PatchBased/BoundaryTrians.jl @@ -14,6 +14,12 @@ using GridapSolvers.LinearSolvers using GridapSolvers.MultilevelTools using GridapSolvers.PatchBasedSmoothers +function get_cell_size(Ω::Triangulation{Dc}) where Dc + dΩ = Measure(Ω,1) + h_e = Gridap.CellData.get_array(∫(1)dΩ) .^ (1/Dc) + return CellField(h_e,Ω) +end + function weakforms(model) Ω = Triangulation(model) Γ = BoundaryTriangulation(model) @@ -26,17 +32,26 @@ function weakforms(model) n_Γ = get_normal_vector(Γ) n_Λ = get_normal_vector(Λ) + h_e = get_cell_size(Λ) + a1(u,v) = ∫(∇(u)⊙∇(v))dΩ a2(u,v) = ∫((∇(v)⋅n_Γ)⋅u)dΓ a3(u,v) = ∫(jump(u⋅n_Λ)⋅jump(v⋅n_Λ))dΛ - return a1, a2, a3 + a4(u,v) = ∫(h_e*jump(u⋅n_Λ)⋅jump(v⋅n_Λ))dΛ + return a1, a2, a3, a4 end model = CartesianDiscreteModel((0,1,0,1),(2,2)) +model = Gridap.Adaptivity.refine(model) -order = 1 +#order = 1 +#qorder = 2*order+1 +#reffe = ReferenceFE(raviart_thomas,Float64,order-1) +#Vh = FESpace(model,reffe) + +order = 2 qorder = 2*order+1 -reffe = ReferenceFE(raviart_thomas,Float64,order-1) +reffe = ReferenceFE(lagrangian,VectorValue{2,Float64},order) Vh = FESpace(model,reffe) Ω = Triangulation(model) @@ -52,8 +67,8 @@ Ph = PatchFESpace(Vh,PD,reffe) Γp_virtual = BoundaryTriangulation(PD,tags=["boundary","interior"],reverse=true) -a1, a2, a3 = weakforms(model) -ap1, ap2, ap3 = weakforms(PD) +a1, a2, a3, a4 = weakforms(model) +ap1, ap2, ap3, ap4 = weakforms(PD) A1 = assemble_matrix(a1,Vh,Vh) Ap1 = assemble_matrix(ap1,Ph,Ph) @@ -63,3 +78,19 @@ Ap2 = assemble_matrix(ap2,Ph,Ph) A3 = assemble_matrix(a3,Vh,Vh) Ap3 = assemble_matrix(ap3,Ph,Ph) + +A4 = assemble_matrix(a4,Vh,Vh) +Ap4 = assemble_matrix(ap4,Ph,Ph) + + +Λ = SkeletonTriangulation(PD) +dΛ = Measure(Λ, qorder) +n_Λ = get_normal_vector(Λ) + +up = get_trial_fe_basis(Ph) +ap(u,v) = ∫(jump(u⋅n_Λ)⋅jump(v⋅n_Λ))dΛ + +A = assemble_matrix(ap,Ph,Ph) + + +@which is_change_possible(Ω,Λ) diff --git a/test/_dev/PatchBased/CoarsePatchDecompositions.jl b/test/_dev/PatchBased/CoarsePatchDecompositions.jl new file mode 100644 index 00000000..c8b4abaf --- /dev/null +++ b/test/_dev/PatchBased/CoarsePatchDecompositions.jl @@ -0,0 +1,16 @@ +using Gridap +using Gridap.Geometry, Gridap.Arrays + +using GridapSolvers +using GridapSolvers: PatchBasedSmoothers + +cmodel = CartesianDiscreteModel((0,1,0,1),(2,2)) +model = Gridap.Adaptivity.refine(cmodel) + +glue = Gridap.Adaptivity.get_adaptivity_glue(model) + +PD = PatchBasedSmoothers.CoarsePatchDecomposition(model) + +reffe = ReferenceFE(lagrangian,Float64,2) +Vh = FESpace(model,reffe) +Ph = PatchFESpace(Vh,PD,reffe) diff --git a/test/_dev/PatchBased/DistributedPatchFESpacesDebuggingTests.jl b/test/_dev/PatchBased/DistributedPatchFESpacesDebuggingTests.jl index a88f390b..6c937d33 100644 --- a/test/_dev/PatchBased/DistributedPatchFESpacesDebuggingTests.jl +++ b/test/_dev/PatchBased/DistributedPatchFESpacesDebuggingTests.jl @@ -72,7 +72,7 @@ function run(ranks) ap(u,v) = biform(u,v,dΩₚ) lp(v) = liform(v,dΩₚ) - assembler_P = SparseMatrixAssembler(Ph,Ph,FullyAssembledRows()) + assembler_P = SparseMatrixAssembler(Ph,Ph,LocallyAssembled()) Ahp = assemble_matrix(ap,assembler_P,Ph,Ph) fhp = assemble_vector(lp,assembler_P,Ph) diff --git a/test/_dev/PatchBased/PatchSolvers.jl b/test/_dev/PatchBased/PatchSolvers.jl new file mode 100644 index 00000000..e8b17239 --- /dev/null +++ b/test/_dev/PatchBased/PatchSolvers.jl @@ -0,0 +1,80 @@ +using Gridap +using Gridap.Geometry, Gridap.FESpaces, Gridap.ReferenceFEs, Gridap.MultiField + +using GridapSolvers +using GridapSolvers.MultilevelTools +using GridapSolvers.LinearSolvers +using GridapSolvers.PatchBasedSmoothers + +using LinearAlgebra + +function mean_value(p,dΩ) + sum(∫(p)dΩ) +end + +model = UnstructuredDiscreteModel(CartesianDiscreteModel((0,1,0,1),(4,4))) + +order = 2 +qdegree = 2*order+1 +reffe_u = LagrangianRefFE(VectorValue{2,Float64},QUAD,order) +reffe_p = LagrangianRefFE(Float64,QUAD,order-1;space=:P) +reffe_λ = LagrangianRefFE(Float64,QUAD,0) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="boundary") +Qh = TestFESpace(model,reffe_p,conformity=:L2) +Λh = ConstantFESpace(model) +Xh = MultiFieldFESpace([Vh,Qh,Λh]) + +qdegree = 2*order+1 +f(x) = VectorValue(1.0,1.0) +Π = LocalProjectionMap(divergence,Qh,qdegree) +a(u,v,dΩ) = ∫(∇(u)⊙∇(v))dΩ +b((u,p),(v,q),dΩ) = ∫(p*(∇⋅v) + q*(∇⋅u))dΩ +d(u,v,dΩ) = ∫(Π(v)*Π(u))dΩ +c((p,λ),(q,μ),dΩ) = ∫(λ*q + μ*p)dΩ + +β = 1.0 +biform((u,p,λ),(v,q,μ),dΩ) = a(u,v,dΩ) + b((u,p),(v,q),dΩ) + β*c((p,λ),(q,μ),dΩ) #+ d(u,v,dΩ) +liform((v,q),dΩ) = ∫(f⋅v)dΩ + +PD = PatchDecomposition(model,patch_boundary_style=PatchBasedSmoothers.PatchBoundaryInclude()) + +Vp = PatchFESpace(Vh,PD,reffe_u) +Qp = PatchFESpace(Qh,PD,reffe_p) +Λp = PatchFESpace(Λh,PD,reffe_λ) +Zh = MultiFieldFESpace([Vp,Qp,Λp]) + +Ω = Triangulation(model) +dΩ = Measure(Ω,4) + +Ωp = Triangulation(PD) +dΩp = Measure(Ωp,4) + +ap(u,v) = a(u,v,dΩp) +ps = PatchBasedSmoothers.PatchSolver(Vh,PD,ap,reffe_u) + +patch_ids = ps.patch_ids +patch_cell_lids = ps.patch_cell_lids + +A = assemble_matrix((u,v)->a(u,v,dΩ),Vh,Vh) +ss = symbolic_setup(ps,A) +ns = numerical_setup(ss,A) + +rhs = assemble_vector(v -> ∫(f⋅v)dΩ,Vh) +x = zeros(num_free_dofs(Vh)) +solve!(x,ns,rhs) + +biformp((u,p,λ),(v,q,μ)) = biform((u,p,λ),(v,q,μ),dΩp) +ps_mf = PatchBasedSmoothers.PatchSolver(Xh,PD,biformp,[reffe_u,reffe_p,reffe_λ]) + +patch_ids = ps_mf.patch_ids +patch_cell_lids = ps_mf.patch_cell_lids + +A = assemble_matrix((u,v)->biform(u,v,dΩ),Xh,Xh) +ss = symbolic_setup(ps_mf,A) +ns = numerical_setup(ss,A) + +rhs = assemble_vector(y -> liform(y,dΩ),Xh) +x = zeros(num_free_dofs(Xh)) +solve!(x,ns,rhs) + diff --git a/test/_dev/PatchBased/ZeroMeanPatchBasedFESpaces.jl b/test/_dev/PatchBased/ZeroMeanPatchBasedFESpaces.jl new file mode 100644 index 00000000..63e97732 --- /dev/null +++ b/test/_dev/PatchBased/ZeroMeanPatchBasedFESpaces.jl @@ -0,0 +1,96 @@ + +using Gridap +using GridapSolvers + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] + x[2] + +model = CartesianDiscreteModel((0,1,0,1),(2,2)) +PD = PatchDecomposition(model) + +Ω = Triangulation(model) +Ωp = Triangulation(PD) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +Vh = TestFESpace(model,reffe_u,dirichlet_tags="boundary") +Uh = TrialFESpace(Vh,u_exact) + +Qh = TestFESpace(model,reffe_p,conformity=:L2,constraint=:zeromean) + +topo = get_grid_topology(model) +#patches_mask = fill(false,PatchBasedSmoothers.num_patches(PD)) +patches_mask = get_isboundary_face(topo,0) +Ph = PatchFESpace(Vh,PD,reffe_u;conformity=H1Conformity(),patches_mask=patches_mask) +Lh = PatchFESpace(Qh,PD,reffe_p;conformity=L2Conformity(),patches_mask=patches_mask) + +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) +Zh = MultiFieldFESpace([Ph,Lh]) + +qdegree = 2*order +dΩ = Measure(Ω,qdegree) +dΩp = Measure(Ωp,qdegree) + +p_mean = sum(∫(p_exact)dΩ)/sum(∫(1)dΩ) +p_exact_zm(x) = p_exact(x) - p_mean +f(x) = Δ(u_exact)(x) + ∇(p_exact_zm)(x) +liform((v,q),dΩ) = ∫(v⋅f)dΩ +biform((u,p),(v,q),dΩ) = ∫(∇(u)⊙∇(v) - (∇⋅v)*p - (∇⋅u)*q)dΩ +a(x,y) = biform(x,y,dΩ) +l(y) = liform(y,dΩ) +ap(x,y) = biform(x,y,dΩp) + +op = AffineFEOperator(a,l,Xh,Yh) +A = get_matrix(op) +b = get_vector(op) +x_exact = A\b + +P = PatchBasedSmoothers.PatchBasedLinearSolver(ap,Zh,Yh) +P_ns = numerical_setup(symbolic_setup(P,A),A) +x = zeros(num_free_dofs(Yh)) +solve!(x,P_ns,b) + +solver = FGMRESSolver(10,P;verbose=true) +#solver = GMRESSolver(10,verbose=true) +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = zeros(num_free_dofs(Yh)) +solve!(x,ns,b) + +norm(A*x - b) + +############################################################################################ + +patch_pcells = get_patch_cells_overlapped(PD) +pcell_to_pdofs_u = Ph.patch_cell_dofs_ids +pcell_to_pdofs_p = Lh.patch_cell_dofs_ids + +patch_dofs_u = map(pcells -> unique(filter(x -> x > 0, vcat(pcell_to_pdofs_u[pcells]...))),patch_pcells) +patch_dofs_p = map(pcells -> unique(filter(x -> x > 0, vcat(pcell_to_pdofs_p[pcells]...))),patch_pcells) + +o = num_free_dofs(Ph) +patch_dofs = map((du,dp) -> vcat(du,dp .+ o),patch_dofs_u,patch_dofs_p) + +patch_cells = PD.patch_cells +patch_models = map(patch_cells) do cells + pmodel = DiscreteModelPortion(model,cells) + pgrid = UnstructuredGrid(get_grid(pmodel)) + ptopo = get_grid_topology(pmodel) + plabels = FaceLabeling(ptopo) + UnstructuredDiscreteModel(pgrid,ptopo,plabels) +end + +patch_spaces = map(patch_models) do pmodel + Vhi = TestFESpace(pmodel,reffe_u;dirichlet_tags="boundary") + Qhi = TestFESpace(pmodel,reffe_p,conformity=:L2,constraint=:zeromean) + Yhi = MultiFieldFESpace([Vhi,Qhi]) +end + +Ap = assemble_matrix(ap,Zh,Zh) +patch_mats = map(dofs -> Ap[dofs,dofs],patch_dofs) diff --git a/test/_dev/PatchBased/patch_closures.jl b/test/_dev/PatchBased/patch_closures.jl new file mode 100644 index 00000000..82014750 --- /dev/null +++ b/test/_dev/PatchBased/patch_closures.jl @@ -0,0 +1,121 @@ +using Gridap +using Gridap.Geometry, Gridap.Arrays + +using GridapSolvers +using GridapSolvers: PatchBasedSmoothers + +function generate_dual_graph( + topo::GridTopology{Dc}, D::Integer = Dc +) where Dc + @assert 0 < D <= Dc + edge_to_face = Geometry.get_faces(topo,D-1,D) + n_faces = Geometry.num_faces(topo,D) + return generate_dual_graph(edge_to_face,n_faces) +end + +# Given a table `edge_to_face`, creates the dual graph `face_to_face`. +function generate_dual_graph( + edge_to_face::Table, + n_faces = maximum(edge_to_face.data) +) + n_edges = length(edge_to_face) + + ptrs = zeros(Int,n_faces+1) + for e in 1:n_edges + faces = view(edge_to_face,e) + if length(faces) > 1 + @assert length(faces) == 2 + f1, f2 = faces + ptrs[f1+1] += 1 + ptrs[f2+1] += 1 + end + end + Arrays.length_to_ptrs!(ptrs) + + data = zeros(Int,ptrs[end]-1) + for e in 1:n_edges + faces = view(edge_to_face,e) + if length(faces) > 1 + f1, f2 = faces + data[ptrs[f1]] = f2 + data[ptrs[f2]] = f1 + ptrs[f1] += 1 + ptrs[f2] += 1 + end + end + Arrays.rewind_ptrs!(ptrs) + + return Table(data,ptrs) +end + +model = CartesianDiscreteModel((0,1,0,1),(2,2)) +topo = get_grid_topology(model) + +PD = PatchDecomposition(model;patch_boundary_style=PatchBasedSmoothers.PatchBoundaryInclude()) + +order = 2 +reffe = ReferenceFE(lagrangian,Float64,order) +Vh = FESpace(model,reffe) +Ph = PatchFESpace(Vh,PD,reffe) + +Ω = Triangulation(model) +Ωp = Triangulation(PD) +Ωc = Closure(PD) + +dΩ = Measure(Ω,2*order) +dΩp = Measure(Ωp,2*order) +dΩc = Measure(Ωc,2*order) + +biform(u,v,dΩ) = ∫(u⋅v)dΩ +a(u,v) = biform(u,v,dΩ) +ap(u,v) = biform(u,v,dΩp) +ac(u,v) = biform(u,v,dΩc) + +A = assemble_matrix(a,Vh,Vh) +Ap = assemble_matrix(ap,Ph,Ph) +Ac = assemble_matrix(ac,Ph,Ph) + +vanka_PD = PatchDecomposition(model;patch_boundary_style=PatchBasedSmoothers.PatchBoundaryExclude()) +vanka = PatchBasedSmoothers.VankaSolver(Vh,vanka_PD) +vanka_ids = vanka.patch_ids + +dof_to_pdof = Ph.dof_to_pdof +patch_cell_ids = get_cell_dof_ids(Ph) +pdof_to_dof = flatten_partition(dof_to_pdof) + +for (patch,ids) in enumerate(vanka_ids) + patch_ids = unique(vcat(PatchBasedSmoothers.patch_view(PD,patch_cell_ids,patch)...)) + + perm = sortperm(pdof_to_dof[patch_ids]) + patch_ids = patch_ids[perm] + + println("> Patch $patch") + println(" > Vanka: ",ids) + println(" > Space: ",pdof_to_dof[patch_ids]) + @assert ids == pdof_to_dof[patch_ids] + + A_vanka = A[ids,ids] + A_vanka_bis = Ac[patch_ids,patch_ids] + @assert A_vanka ≈ A_vanka_bis +end + +b(u,v) = ∫(u⋅v)dΩc + ∫(u⋅v)dΩp +B = assemble_matrix(b,Ph,Ph) +@assert B ≈ Ac + Ap + +# Multifield + +Xh = MultiFieldFESpace([Vh,Vh]) +Zh = MultiFieldFESpace([Ph,Ph]) + +biform_mf((u1,u2),(v1,v2),dΩ) = ∫(u1⋅v1 + u2⋅v2)dΩ +a_mf(u,v) = biform_mf(u,v,dΩ) +ap_mf(u,v) = biform_mf(u,v,dΩp) +ac_mf(u,v) = biform_mf(u,v,dΩc) + +A = assemble_matrix(a_mf,Xh,Xh) +Ap = assemble_matrix(ap_mf,Zh,Zh) +Ac = assemble_matrix(ac_mf,Zh,Zh) + + + diff --git a/test/_dev/Vanka/darcy_patch.jl b/test/_dev/Vanka/darcy_patch.jl new file mode 100644 index 00000000..4c1bee29 --- /dev/null +++ b/test/_dev/Vanka/darcy_patch.jl @@ -0,0 +1,69 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +model = CartesianDiscreteModel((0,1,0,1),(2,2)) +labels = get_face_labeling(model) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) + +order = 1 +reffe_u = ReferenceFE(raviart_thomas,Float64,order-1) +reffe_p = ReferenceFE(lagrangian,Float64,order-1) + +Vh = TestFESpace(model,reffe_u) +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) + +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) + +qdegree = 2*order +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) + +Γ = BoundaryTriangulation(model,tags="boundary") +dΓ = Measure(Γ,qdegree) +n = get_normal_vector(Γ) + +α = 10.0 +f(x) = u_exact(x) - ∇(p_exact)(x) +σ(x) = p_exact(x) +function a((u,p),(v,q)) + c = ∫(u⋅v + (∇⋅v)*p + (∇⋅u)*q)dΩ + if !iszero(α) + c += ∫((∇⋅u)⋅(∇⋅v))*dΩ + end + return c +end +l((v,q)) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n) )dΓ + +op = AffineFEOperator(a,l,Xh,Yh) +A = get_matrix(op) +b = get_vector(op) +cond(Matrix(A)) +x_exact = A\b + +uh_exact, ph_exact = FEFunction(Xh,x_exact) +l2_error(uh_exact,u_exact,dΩ) +l2_error(ph_exact,p_exact,dΩ) +l2_norm(∇⋅uh_exact,dΩ) + diff --git a/test/_dev/Vanka/darcy_vanka.jl b/test/_dev/Vanka/darcy_vanka.jl new file mode 100644 index 00000000..bfd185f3 --- /dev/null +++ b/test/_dev/Vanka/darcy_vanka.jl @@ -0,0 +1,79 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +model = CartesianDiscreteModel((0,1,0,1),(8,8)) +labels = get_face_labeling(model) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) + +order = 1 +reffe_u = ReferenceFE(raviart_thomas,Float64,order-1) +reffe_p = ReferenceFE(lagrangian,Float64,order-1) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) + +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) + +qdegree = 2*order +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) + +Γ = BoundaryTriangulation(model,tags="newman") +dΓ = Measure(Γ,qdegree) +n = get_normal_vector(Γ) + +f(x) = u_exact(x) - ∇(p_exact)(x) +σ(x) = p_exact(x) +a((u,p),(v,q)) = ∫(u⋅v + (∇⋅v)*p - (∇⋅u)*q)dΩ +l((v,q)) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n) )dΓ + +op = AffineFEOperator(a,l,Xh,Yh) +A = get_matrix(op) +b = get_vector(op) +cond(Matrix(A)) +x_exact = A\b + +uh_exact, ph_exact = FEFunction(Xh,x_exact) +l2_error(uh_exact,u_exact,dΩ) +l2_error(ph_exact,p_exact,dΩ) +l2_norm(∇⋅uh_exact,dΩ) + +PD = PatchDecomposition(model) +P = VankaSolver(Xh,PD) +P_ns = numerical_setup(symbolic_setup(P,A),A) + +x = zeros(num_free_dofs(Yh)) +r = b - A*x +solve!(x,P_ns,r) +norm(x_exact - x) + + +solver = FGMRESSolver(10,P;rtol=1.e-8,verbose=true) +#solver = GMRESSolver(10,verbose=true) +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = zeros(num_free_dofs(Yh)) +solve!(x,ns,b) +norm(A*x - b) diff --git a/test/_dev/Vanka/darcy_vanka_GMG.jl b/test/_dev/Vanka/darcy_vanka_GMG.jl new file mode 100644 index 00000000..68ae5bc3 --- /dev/null +++ b/test/_dev/Vanka/darcy_vanka_GMG.jl @@ -0,0 +1,133 @@ +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays, Gridap.Adaptivity +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +cmodel = CartesianDiscreteModel((0,1,0,1),(4,4)) +labels = get_face_labeling(cmodel) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) +model = refine(cmodel,2) + +order = 1 +reffe_u = ReferenceFE(raviart_thomas,Float64,order-1) +reffe_p = ReferenceFE(lagrangian,Float64,order-1) + +VH = TestFESpace(cmodel,reffe_u;dirichlet_tags="dirichlet") +UH = TrialFESpace(VH,u_exact) +QH = TestFESpace(cmodel,reffe_p,conformity=:L2) +XH = MultiFieldFESpace([UH,QH]) +YH = MultiFieldFESpace([VH,QH]) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) + +qdegree = 2*order +Ωh = Triangulation(model) +dΩh = Measure(Ωh,qdegree) +ΩH = Triangulation(cmodel) +dΩH = Measure(ΩH,qdegree) +dΩHh = Measure(ΩH,Ωh,qdegree) + +Γh = BoundaryTriangulation(model,tags="newman") +dΓh = Measure(Γh,qdegree) +n = get_normal_vector(Γh) + +α = 10.0 +f(x) = u_exact(x) - ∇(p_exact)(x) +σ(x) = p_exact(x) +graddiv(u,v,dΩ) = ∫((∇⋅u)⋅(∇⋅v))*dΩ +function a((u,p),(v,q),dΩ) + c = ∫(u⋅v + (∇⋅v)*p + (∇⋅u)*q)dΩ + if α > 0.0 + c += graddiv(u,v,dΩ) + end + return c +end +l((v,q),dΩ,dΓ) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n))dΓ + +ah(x,y) = a(x,y,dΩh) +aH(x,y) = a(x,y,dΩH) +lh(y) = l(y,dΩh,dΓh) + +op_h = AffineFEOperator(ah,lh,Xh,Yh) +Ah = get_matrix(op_h) +bh = get_vector(op_h) +xh_exact = Ah\bh + +AH = assemble_matrix(aH,XH,YH) +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Xh,Xh) + +uh_exact, ph_exact = FEFunction(Xh,xh_exact) +l2_error(uh_exact,u_exact,dΩh) +l2_error(ph_exact,p_exact,dΩh) +l2_norm(∇⋅uh_exact,dΩh) + +# NOTE: Convergence seem quite dependent on the damping weight. For ω=0.2, we diverge... +PD = PatchDecomposition(model) +smoother = RichardsonSmoother(VankaSolver(Xh,PD),10,0.1) +smoother_ns = numerical_setup(symbolic_setup(smoother,Ah),Ah) + +function project_f2c(rh) + Qrh = Mhh\rh + uh, ph = FEFunction(Yh,Qrh) + ll((v,q)) = ∫(v⋅uh + q*ph)*dΩHh + assemble_vector(ll,YH) +end +function interp_c2f(xH) + get_free_dof_values(interpolate(FEFunction(YH,xH),Yh)) +end + +xh = randn(size(Ah,2)) +rh = bh - Ah*xh +niters = 10 + +iter = 0 +error0 = norm(rh) +error = error0 +e_rel = error/error0 +while iter < niters && e_rel > 1.0e-10 + println("Iter $iter:") + println(" > Initial: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + println(" > Pre-smoother: ", norm(rh)) + + rH = project_f2c(rh) + qH = AH\rH + qh = interp_c2f(qH) + + rh = rh - Ah*qh + xh = xh + qh + println(" > Post-correction: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + + iter += 1 + error = norm(rh) + e_rel = error/error0 + println(" > Final: ",error, " - ", e_rel) +end + +uh, ph = FEFunction(Xh,xh) +l2_error(uh,u_exact,dΩh) +l2_error(ph,p_exact,dΩh) diff --git a/test/_dev/Vanka/distributed_vanka.jl b/test/_dev/Vanka/distributed_vanka.jl new file mode 100644 index 00000000..116acc30 --- /dev/null +++ b/test/_dev/Vanka/distributed_vanka.jl @@ -0,0 +1,32 @@ + +using Gridap +using PartitionedArrays +using GridapDistributed + + +np = (2,1) +ranks = with_debug() do distribute + distribute(LinearIndices((prod(np),))) +end + +model = CartesianDiscreteModel(ranks,np,(0,1,0,1),(4,4)) + +order = 1 +reffe = ReferenceFE(lagrangian,Float64,order) +Vh = TestFESpace(model,reffe) + +Ω = Triangulation(model) +dΩ = Measure(Ω,2*order) + +a(u,v) = ∫(v⋅u)*dΩ +A = assemble_matrix(a,Vh,Vh) + + +rows, cols = axes(A) + +nbors_snd, nbors_rcv = assembly_neighbors(partition(rows)) + +map(partition(rows),partition(cols),partition(A)) do rows,cols,mat + own_rows = own_to_local(rows) +end + diff --git a/test/_dev/Vanka/kernel_analysis.jl b/test/_dev/Vanka/kernel_analysis.jl new file mode 100644 index 00000000..89f4ef7a --- /dev/null +++ b/test/_dev/Vanka/kernel_analysis.jl @@ -0,0 +1,70 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers +using GridapSolvers.MultilevelTools + +Dc = 2 +model = CartesianDiscreteModel((0,1,0,1),(4,4)) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +V = TestFESpace(model,reffe_u) +Q = TestFESpace(model,reffe_p,conformity=:L2) +X = MultiFieldFESpace([V,Q]) + +qdegree = 2*order +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) + +PD = PatchDecomposition(model) +PD.patch_cells + +ν = 1.e-8 +η = -1.0 +ε = 1.e-4 +Π = LocalProjectionMap(divergence,reffe_p) +#graddiv(u,v,dΩ) = ∫(η*(∇⋅v)⋅Π(u))*dΩ +graddiv(u,v,dΩ) = ∫(η*Π(v)⋅Π(u))*dΩ +lap(u,v,dΩ) = ∫(ν*∇(u)⊙∇(v))*dΩ +mass(u,v,dΩ) = ∫(u*v)*dΩ + +function a((u,p),(v,q),dΩ) + c = lap(u,v,dΩ) + c += ∫((∇⋅v)*p + (∇⋅u)*q)dΩ + if η > 0.0 + c += graddiv(u,v,dΩ) + end + if ε > 0.0 + c += mass(p,q,dΩ) + end + return c +end + +A = assemble_matrix((u,v) -> a(u,v,dΩ),X,X) +x = A\ones(size(A,1)) + +cells = [6,7,10,11] +dof_ids = get_cell_dof_ids(X) +patch_ids = unique(sort(vcat(map(ids -> vcat(ids.array...), dof_ids[cells])...))) + +A_vanka = Matrix(A[patch_ids,patch_ids]) +cond(A_vanka) +x_vanka = A_vanka\randn(size(A_vanka,1)) + +s = VankaSolver(X,PD) +ns = numerical_setup(symbolic_setup(s,A),A) + +for i in 1:10 + x = zeros(size(A,1)) + x_exact = randn(size(A,1)) + b = A*x_exact + solve!(x,ns,b) + println(norm(x-x_exact)) +end diff --git a/test/_dev/Vanka/mhd_uj_vanka_GMG.jl b/test/_dev/Vanka/mhd_uj_vanka_GMG.jl new file mode 100644 index 00000000..7a1569f8 --- /dev/null +++ b/test/_dev/Vanka/mhd_uj_vanka_GMG.jl @@ -0,0 +1,223 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays, Gridap.Adaptivity +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +function add_labels_2d(model) + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"newman",[5]) + add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) +end + +function add_labels_3d(model) + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"newman",[21]) + add_tag_from_tags!(labels,"dirichlet",[collect(1:20)...,22,23,24,25,26]) +end + +B = VectorValue(0.0,0.0,1.0) +u_exact(x) = VectorValue(x[1],-x[2],0.0) +p_exact(x) = sum(x) +j_exact(x) = VectorValue(-x[2],-x[1],0.0) + VectorValue(1.0,1.0,0.0) +φ_exact(x) = x[1] + x[2] + +Dc = 3 +domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) +ncells = (Dc == 2) ? (4,4) : (2,2,2) +cmodel = CartesianDiscreteModel(domain,ncells) +if Dc == 2 + add_labels_2d(cmodel) +else + add_labels_3d(cmodel) +end +model = refine(cmodel,2) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{Dc,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) +reffe_j = ReferenceFE(raviart_thomas,Float64,order-1) +reffe_φ = ReferenceFE(lagrangian,Float64,order-1) + +using Badia2024 +using SpecialPolynomials +using Polynomials +reffe_j = Badia2024.SpecialRTReffe(ChebyshevT,order,Dc) + +VH = TestFESpace(cmodel,reffe_u;dirichlet_tags="dirichlet") +UH = TrialFESpace(VH,u_exact) +DH = TestFESpace(cmodel,reffe_j;dirichlet_tags="dirichlet") +JH = TrialFESpace(DH,j_exact) +XH = MultiFieldFESpace([UH,JH]) +YH = MultiFieldFESpace([VH,DH]) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Dh = TestFESpace(model,reffe_j;dirichlet_tags="dirichlet") +Jh = TrialFESpace(Dh,j_exact) +Xh = MultiFieldFESpace([Uh,Jh]) +Yh = MultiFieldFESpace([Vh,Dh]) + +qdegree = 2*order +Ωh = Triangulation(model) +dΩh = Measure(Ωh,qdegree) +ΩH = Triangulation(cmodel) +dΩH = Measure(ΩH,qdegree) +dΩHh = Measure(ΩH,Ωh,qdegree) + +Γh = BoundaryTriangulation(model,tags="newman") +dΓh = Measure(Γh,qdegree) +n = get_normal_vector(Γh) + +ν = 1.e-8 +α = 1.0 +f(x) = -(-Δ(u_exact)(x) - cross(j_exact(x),B)) +g(x) = j_exact(x) - cross(u_exact(x),B) +σ(x) = ∇(u_exact)(x) + +crossB(u,v,dΩ) = ∫(cross(u,B)⋅v)dΩ +Π = GridapSolvers.MultilevelTools.LocalProjectionMap(divergence,reffe_p,qdegree) +Πgraddiv(u,v,dΩ) = ∫(α*Π(u)⋅(∇⋅v))*dΩ +graddiv(u,v,dΩ) = ∫(α*(∇⋅u)⋅(∇⋅v))*dΩ +function a((u,j),(v,d),dΩ) + c = ∫((-ν)*∇(u)⊙∇(v) + j⋅d)dΩ + c = c - crossB(u,d,dΩ) + crossB(j,v,dΩ) + if α > 0.0 + c += Πgraddiv(u,v,dΩ) #+ graddiv(j,d,dΩ) + end + return c +end +function a_u(u,v,dΩ) + c = ∫(ν*∇(u)⊙∇(v))*dΩ + if α > 0.0 + c += Πgraddiv(u,v,dΩ) + end + return c +end +l((v,d),dΩ,dΓ) = ∫(v⋅f + d⋅g)dΩ + ∫(v⋅(σ⋅n))dΓ + +ah(x,y) = a(x,y,dΩh) +aH(x,y) = a(x,y,dΩH) +lh(y) = l(y,dΩh,dΓh) + +op_h = AffineFEOperator(ah,lh,Xh,Yh) +Ah = get_matrix(op_h) +bh = get_vector(op_h) +xh_exact = Ah\bh + +AH = assemble_matrix(aH,XH,YH) +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Xh,Xh) + +uh_exact, jh_exact = FEFunction(Xh,xh_exact) +l2_error(uh_exact,u_exact,dΩh) +l2_error(jh_exact,j_exact,dΩh) +l2_norm(∇⋅uh_exact,dΩh) +l2_norm(∇⋅jh_exact,dΩh) + +PD = PatchDecomposition(model) +smoother = RichardsonSmoother(VankaSolver(Xh,PD),10,ν) +smoother_ns = numerical_setup(symbolic_setup(smoother,Ah),Ah) + +function project_f2c(rh) + Qrh = Mhh\rh + uh, jh = FEFunction(Yh,Qrh) + ll((v,d)) = ∫(v⋅uh + d⋅jh)*dΩHh + assemble_vector(ll,YH) +end +function interp_c2f(xH) + get_free_dof_values(interpolate(FEFunction(YH,xH),Yh)) +end + +############################################################################################ + +patches_mask = PatchBasedSmoothers.get_coarse_node_mask(model,model.glue) +Ih = PatchFESpace(Vh,PD,reffe_u;conformity=H1Conformity(),patches_mask=patches_mask) + +Ωp = Triangulation(PD) +dΩp = Measure(Ωp,qdegree) +ap(x,y) = a_u(x,y,dΩp) +Ai = assemble_matrix(ap,Ih,Ih) + +function P1(dxH) + interp_c2f(dxH) +end +function R1(rh) + project_f2c(rh) +end + +function P2(dxH) + xh = interpolate(FEFunction(YH,dxH),Yh) + dxh = get_free_dof_values(xh) + + uh, jh = xh + r̃h = assemble_vector(v -> Πgraddiv(uh,v,dΩp),Ih) + dx̃ = Ai\r̃h + Pdxh = zero_free_values(Xh) + Pdxh_u = Gridap.MultiField.restrict_to_field(Xh,Pdxh,1) + PatchBasedSmoothers.inject!(Pdxh_u,Ih,dx̃) + y = dxh - Pdxh + return y +end +function R2(rh) + rh_u = Gridap.MultiField.restrict_to_field(YH,rh,1) + r̃h = zero_free_values(Ih) + PatchBasedSmoothers.prolongate!(r̃h,Ih,rh_u) + dr̃h = Ai\r̃h + + Pdxh = zero_free_values(Vh) + PatchBasedSmoothers.inject!(Pdxh,Ih,dr̃h) + uh = FEFunction(Vh,Pdxh) + ll((v,q)) = Πgraddiv(uh,v,dΩh) + drh = assemble_vector(ll,Yh) + rH = project_f2c(rh - drh) + return rH +end + +############################################################################################ + +xh = zeros(size(Ah,2)) +rh = bh - Ah*xh +niters = 10 + +iter = 0 +error0 = norm(rh) +error = error0 +e_rel = error/error0 +while iter < niters && e_rel > 1.0e-10 + println("Iter $iter:") + println(" > Initial: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + println(" > Pre-smoother: ", norm(rh)) + + rH = R2(rh) + qH = AH\rH + qh = P2(qH) + + rh = rh - Ah*qh + xh = xh + qh + println(" > Post-correction: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + + iter += 1 + error = norm(rh) + e_rel = error/error0 + println(" > Final: ",error, " - ", e_rel) +end + +uh, jh = FEFunction(Xh,xh) +l2_error(uh,u_exact,dΩh) +l2_error(jh,j_exact,dΩh) diff --git a/test/_dev/Vanka/mhd_vanka_GMG.jl b/test/_dev/Vanka/mhd_vanka_GMG.jl new file mode 100644 index 00000000..771b729b --- /dev/null +++ b/test/_dev/Vanka/mhd_vanka_GMG.jl @@ -0,0 +1,175 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays, Gridap.Adaptivity +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +function add_labels_2d(model) + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"newman",[5]) + add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) +end + +function add_labels_3d(model) + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"newman",[21]) + add_tag_from_tags!(labels,"dirichlet",[collect(1:20)...,22,23,24,25,26]) +end + +B = VectorValue(0.0,0.0,1.0) +u_exact(x) = VectorValue(x[1],-x[2],0.0) +p_exact(x) = sum(x) +j_exact(x) = VectorValue(-x[2],-x[1],0.0) + VectorValue(1.0,1.0,0.0) +φ_exact(x) = x[1] + x[2] + +Dc = 3 +domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) +ncells = (Dc == 2) ? (4,4) : (2,2,2) +cmodel = CartesianDiscreteModel(domain,ncells) +if Dc == 2 + add_labels_2d(cmodel) +else + add_labels_3d(cmodel) +end +model = refine(cmodel,2) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{Dc,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) +reffe_j = ReferenceFE(raviart_thomas,Float64,order-1) +reffe_φ = ReferenceFE(lagrangian,Float64,order-1) + +VH = TestFESpace(cmodel,reffe_u;dirichlet_tags="dirichlet") +UH = TrialFESpace(VH,u_exact) +QH = TestFESpace(cmodel,reffe_p,conformity=:L2) +DH = TestFESpace(cmodel,reffe_j;dirichlet_tags="dirichlet") +JH = TrialFESpace(DH,j_exact) +ΦH = TestFESpace(cmodel,reffe_φ,conformity=:L2) + +XH = MultiFieldFESpace([UH,QH,JH,ΦH]) +YH = MultiFieldFESpace([VH,QH,DH,ΦH]) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) +Dh = TestFESpace(model,reffe_j;dirichlet_tags="dirichlet") +Jh = TrialFESpace(Dh,j_exact) +Φh = TestFESpace(model,reffe_φ,conformity=:L2) +Xh = MultiFieldFESpace([Uh,Qh,Jh,Φh]) +Yh = MultiFieldFESpace([Vh,Qh,Dh,Φh]) + +qdegree = 2*order +Ωh = Triangulation(model) +dΩh = Measure(Ωh,qdegree) +ΩH = Triangulation(cmodel) +dΩH = Measure(ΩH,qdegree) +dΩHh = Measure(ΩH,Ωh,qdegree) + +Γh = BoundaryTriangulation(model,tags="newman") +dΓh = Measure(Γh,qdegree) +n = get_normal_vector(Γh) + +ν = 1.e-2 +α = -1.0 +I_tensor = one(TensorValue{Dc,Dc,Float64}) +f(x) = -Δ(u_exact)(x) - ∇(p_exact)(x) - cross(j_exact(x),B) +g(x) = j_exact(x) - ∇(φ_exact)(x) - cross(u_exact(x),B) +σ(x) = ∇(u_exact)(x) + p_exact(x)*I_tensor +γ(x) = φ_exact(x)*I_tensor +crossB(u,v,dΩ) = ∫(cross(u,B)⋅v)dΩ +Π = GridapSolvers.MultilevelTools.LocalProjectionMap(divergence,reffe_p,qdegree) +Πgraddiv(u,v,dΩ) = ∫(α*Π(u)⋅(∇⋅v))*dΩ +graddiv(u,v,dΩ) = ∫(α*(∇⋅u)⋅(∇⋅v))*dΩ +function a((u,p,j,φ),(v,q,d,ψ),dΩ) + c = ∫(ν*∇(u)⊙∇(v) + (∇⋅v)*p - (∇⋅u)*q)dΩ + c = c + ∫(j⋅d + (∇⋅d)*φ - (∇⋅j)*ψ)dΩ + c = c - crossB(u,d,dΩ) - crossB(j,v,dΩ) + if α > 0.0 + c += Πgraddiv(u,v,dΩ) + graddiv(j,d,dΩ) + end + return c +end +l((v,q,d,ψ),dΩ,dΓ) = ∫(v⋅f + d⋅g)dΩ + ∫(v⋅(σ⋅n) + d⋅(γ⋅n))dΓ + +ah(x,y) = a(x,y,dΩh) +aH(x,y) = a(x,y,dΩH) +lh(y) = l(y,dΩh,dΓh) + +op_h = AffineFEOperator(ah,lh,Xh,Yh) +Ah = get_matrix(op_h) +bh = get_vector(op_h) +xh_exact = Ah\bh + +AH = assemble_matrix(aH,XH,YH) +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Xh,Xh) + +uh_exact, ph_exact, jh_exact, φh_exact = FEFunction(Xh,xh_exact) +l2_error(uh_exact,u_exact,dΩh) +l2_error(ph_exact,p_exact,dΩh) +l2_error(jh_exact,j_exact,dΩh) +l2_error(φh_exact,φ_exact,dΩh) +l2_norm(∇⋅uh_exact,dΩh) +l2_norm(∇⋅jh_exact,dΩh) + +PD = PatchDecomposition(model) +smoother = RichardsonSmoother(VankaSolver(Xh,PD),10,0.05) +smoother_ns = numerical_setup(symbolic_setup(smoother,Ah),Ah) + +function project_f2c(rh) + Qrh = Mhh\rh + uh, ph, jh, φh = FEFunction(Yh,Qrh) + ll((v,q,d,ψ)) = ∫(v⋅uh + q*ph + d⋅jh + ψ⋅φh)*dΩHh + assemble_vector(ll,YH) +end +function interp_c2f(xH) + get_free_dof_values(interpolate(FEFunction(YH,xH),Yh)) +end + +xh = randn(size(Ah,2)) +rh = bh - Ah*xh +niters = 10 + +iter = 0 +error0 = norm(rh) +error = error0 +e_rel = error/error0 +while iter < niters && e_rel > 1.0e-10 + println("Iter $iter:") + println(" > Initial: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + println(" > Pre-smoother: ", norm(rh)) + + rH = project_f2c(rh) + qH = AH\rH + qh = interp_c2f(qH) + + rh = rh - Ah*qh + xh = xh + qh + println(" > Post-correction: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + + iter += 1 + error = norm(rh) + e_rel = error/error0 + println(" > Final: ",error, " - ", e_rel) +end + +uh, ph, jh, φh = FEFunction(Xh,xh) +l2_error(uh,u_exact,dΩh) +l2_error(ph,p_exact,dΩh) +l2_error(jh,j_exact,dΩh) +l2_error(φh,φ_exact,dΩh) diff --git a/test/_dev/Vanka/stokes_vanka.jl b/test/_dev/Vanka/stokes_vanka.jl new file mode 100644 index 00000000..67f9a5bc --- /dev/null +++ b/test/_dev/Vanka/stokes_vanka.jl @@ -0,0 +1,80 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +model = CartesianDiscreteModel((0,1,0,1),(8,8)) +labels = get_face_labeling(model) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) + +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) + +qdegree = 2*order +Ω = Triangulation(model) +dΩ = Measure(Ω,qdegree) + +Γ = BoundaryTriangulation(model,tags="newman") +dΓ = Measure(Γ,qdegree) +n = get_normal_vector(Γ) + +I_tensor = one(TensorValue{Dc,Dc,Float64}) +f(x) = -Δ(u_exact)(x) - ∇(p_exact)(x) +σ(x) = ∇(u_exact)(x) + p_exact(x)*I_tensor +a((u,p),(v,q)) = ∫(∇(u)⊙∇(v) + (∇⋅v)*p - (∇⋅u)*q)dΩ +l((v,q)) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n) )dΓ + +op = AffineFEOperator(a,l,Xh,Yh) +A = get_matrix(op) +b = get_vector(op) +cond(Matrix(A)) +x_exact = A\b + +uh_exact, ph_exact = FEFunction(Xh,x_exact) +l2_error(uh_exact,u_exact,dΩ) +l2_error(ph_exact,p_exact,dΩ) +l2_norm(∇⋅uh_exact,dΩ) + +PD = PatchDecomposition(model) +P = VankaSolver(Xh,PD) +P_ns = numerical_setup(symbolic_setup(P,A),A) + +x = zeros(num_free_dofs(Yh)) +r = b - A*x +solve!(x,P_ns,r) +norm(x_exact - x) + + +solver = FGMRESSolver(10,P;rtol=1.e-8,verbose=true) +#solver = GMRESSolver(10,verbose=true) +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = zeros(num_free_dofs(Yh)) +solve!(x,ns,b) +norm(A*x - b) diff --git a/test/_dev/Vanka/stokes_vanka_GMG.jl b/test/_dev/Vanka/stokes_vanka_GMG.jl new file mode 100644 index 00000000..a91838dd --- /dev/null +++ b/test/_dev/Vanka/stokes_vanka_GMG.jl @@ -0,0 +1,189 @@ +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays, Gridap.Adaptivity +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +cmodel = CartesianDiscreteModel((0,1,0,1),(4,4)) +labels = get_face_labeling(cmodel) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) +model = refine(cmodel,2) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +VH = TestFESpace(cmodel,reffe_u;dirichlet_tags="dirichlet") +UH = TrialFESpace(VH,u_exact) +QH = TestFESpace(cmodel,reffe_p,conformity=:L2) +XH = MultiFieldFESpace([UH,QH]) +YH = MultiFieldFESpace([VH,QH]) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) + +qdegree = 2*order +Ωh = Triangulation(model) +dΩh = Measure(Ωh,qdegree) +ΩH = Triangulation(cmodel) +dΩH = Measure(ΩH,qdegree) +dΩHh = Measure(ΩH,Ωh,qdegree) + +Γh = BoundaryTriangulation(model,tags="newman") +dΓh = Measure(Γh,qdegree) +n = get_normal_vector(Γh) + +ν = 1.e-8 +α = 10.0 +I_tensor = one(TensorValue{Dc,Dc,Float64}) +f(x) = -ν*Δ(u_exact)(x) - ∇(p_exact)(x) +σ(x) = ν*∇(u_exact)(x) + p_exact(x)*I_tensor + +Π = GridapSolvers.MultilevelTools.LocalProjectionMap(divergence,reffe_p,qdegree) +#graddiv(u,v,dΩ) = ∫(α*(∇⋅u)⋅(∇⋅v))*dΩ +graddiv(u,v,dΩ) = ∫(α*Π(u)⋅(∇⋅v))*dΩ +function a_u(u,v,dΩ) + c = ∫(ν*∇(u)⊙∇(v))*dΩ + if α > 0.0 + c += graddiv(u,v,dΩ) + end + return c +end +function a((u,p),(v,q),dΩ) + c = ∫(ν*∇(u)⊙∇(v) + (∇⋅v)*p - (∇⋅u)*q)dΩ + if α > 0.0 + c += graddiv(u,v,dΩ) + end + return c +end +l((v,q),dΩ,dΓ) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n) )dΓ + +ah(x,y) = a(x,y,dΩh) +aH(x,y) = a(x,y,dΩH) +lh(y) = l(y,dΩh,dΓh) + +op_h = AffineFEOperator(ah,lh,Xh,Yh) +Ah = get_matrix(op_h) +bh = get_vector(op_h) +xh_exact = Ah\bh + +AH = assemble_matrix(aH,XH,YH) +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Xh,Xh) + +uh_exact, ph_exact = FEFunction(Xh,xh_exact) +l2_error(uh_exact,u_exact,dΩh) +l2_error(ph_exact,p_exact,dΩh) +l2_norm(∇⋅uh_exact,dΩh) + +PD = PatchDecomposition(model) +smoother = RichardsonSmoother(VankaSolver(Xh,PD),10,0.1) +smoother_ns = numerical_setup(symbolic_setup(smoother,Ah),Ah) + +function project_f2c(rh) + Qrh = Mhh\rh + uh, ph = FEFunction(Yh,Qrh) + ll((v,q)) = ∫(v⋅uh + q*ph)*dΩHh + assemble_vector(ll,YH) +end +function interp_c2f(xH) + get_free_dof_values(interpolate(FEFunction(YH,xH),Yh)) +end + +patches_mask = PatchBasedSmoothers.get_coarse_node_mask(model,model.glue) +Ih = PatchFESpace(Vh,PD,reffe_u;conformity=H1Conformity(),patches_mask=patches_mask) + +Ωp = Triangulation(PD) +dΩp = Measure(Ωp,qdegree) +ap(x,y) = a_u(x,y,dΩp) +Ai = assemble_matrix(ap,Ih,Ih) + +function P1(dxH) + interp_c2f(dxH) +end +function R1(rh) + project_f2c(rh) +end + +function P2(dxH) + xh = interpolate(FEFunction(YH,dxH),Yh) + dxh = get_free_dof_values(xh) + + uh, ph = xh + r̃h = assemble_vector(v -> graddiv(uh,v,dΩp),Ih) + dx̃ = Ai\r̃h + Pdxh = zero_free_values(Xh) + Pdxh_u = Gridap.MultiField.restrict_to_field(Xh,Pdxh,1) + PatchBasedSmoothers.inject!(Pdxh_u,Ih,dx̃) + y = dxh - Pdxh + return y +end +function R2(rh) + rh_u = Gridap.MultiField.restrict_to_field(YH,rh,1) + r̃h = zero_free_values(Ih) + PatchBasedSmoothers.prolongate!(r̃h,Ih,rh_u) + dr̃h = Ai\r̃h + + Pdxh = zero_free_values(Vh) + PatchBasedSmoothers.inject!(Pdxh,Ih,dr̃h) + uh = FEFunction(Vh,Pdxh) + ll((v,q)) = graddiv(uh,v,dΩh) + drh = assemble_vector(ll,Yh) + rH = project_f2c(rh - drh) + return rH +end + + +xh = randn(size(Ah,2)) +rh = bh - Ah*xh +niters = 20 + +iter = 0 +error0 = norm(rh) +error = error0 +e_rel = error/error0 +while iter < niters && e_rel > 1.0e-8 + println("Iter $iter:") + println(" > Initial: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + println(" > Pre-smoother: ", norm(rh)) + + rH = R2(rh) + qH = AH\rH + qh = P2(qH) + + rh = rh - Ah*qh + xh = xh + qh + println(" > Post-correction: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + + iter += 1 + error = norm(rh) + e_rel = error/error0 + println(" > Final: ",error, " - ", e_rel) +end + +uh, ph = FEFunction(Xh,xh) +l2_error(uh,u_exact,dΩh) +l2_error(ph,p_exact,dΩh) +l2_norm(∇⋅uh,dΩh) diff --git a/test/_dev/Vanka/stokes_vanka_PCont.jl b/test/_dev/Vanka/stokes_vanka_PCont.jl new file mode 100644 index 00000000..75da9b4c --- /dev/null +++ b/test/_dev/Vanka/stokes_vanka_PCont.jl @@ -0,0 +1,90 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +using GridapSolvers.PatchBasedSmoothers: PatchBoundaryExclude, PatchBoundaryInclude + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +model = CartesianDiscreteModel((0,1,0,1),(8,8)) +labels = get_face_labeling(model) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) + +PD_e = PatchDecomposition(model,patch_boundary_style=PatchBoundaryExclude()) +PD_i = PatchDecomposition(model,patch_boundary_style=PatchBoundaryInclude()) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p) + +PD = PD_i +Ph = PatchFESpace(Vh,PD_i,reffe_u;conformity=H1Conformity()) +Lh = PatchFESpace(Qh,PD_e,reffe_p;conformity=H1Conformity()) + +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) +Zh = MultiFieldFESpace([Ph,Lh]) + +qdegree = 2*order +Ω = Triangulation(model) +Ωp = Triangulation(PD) + +dΩ = Measure(Ω,qdegree) +dΩp = Measure(Ωp,qdegree) + +Γ = BoundaryTriangulation(model,tags="newman") +dΓ = Measure(Γ,qdegree) +n = get_normal_vector(Γ) + +I_tensor = one(TensorValue{Dc,Dc,Float64}) +f(x) = -Δ(u_exact)(x) - ∇(p_exact)(x) +σ(x) = ∇(u_exact)(x) + p_exact(x)*I_tensor +biform((u,p),(v,q),dΩ) = ∫(∇(u)⊙∇(v) + (∇⋅v)*p + (∇⋅u)*q)dΩ +a(x,y) = biform(x,y,dΩ) +l((v,q)) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n) )dΓ + +ap(x,y) = biform(x,y,dΩp) + +op = AffineFEOperator(a,l,Xh,Yh) +A = get_matrix(op) +b = get_vector(op) +cond(Matrix(A)) +x_exact = A\b +uh_exact, ph_exact = FEFunction(Xh,x_exact) +l2_error(uh_exact,u_exact,dΩ) +l2_error(ph_exact,p_exact,dΩ) +l2_norm(∇⋅uh_exact,dΩ) + +P = PatchBasedSmoothers.PatchBasedLinearSolver(ap,Zh,Yh) +P_ns = numerical_setup(symbolic_setup(P,A),A) +x = zeros(num_free_dofs(Yh)) +solve!(x,P_ns,b) + +solver = FGMRESSolver(10,P;rtol=1.e-8,verbose=true) +#solver = GMRESSolver(10,verbose=true) +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = zeros(num_free_dofs(Yh)) +solve!(x,ns,b) + +norm(A*x - b) diff --git a/test/_dev/Vanka/stokes_vanka_Pdisc.jl b/test/_dev/Vanka/stokes_vanka_Pdisc.jl new file mode 100644 index 00000000..5a2e7b8d --- /dev/null +++ b/test/_dev/Vanka/stokes_vanka_Pdisc.jl @@ -0,0 +1,99 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +using GridapSolvers.PatchBasedSmoothers: PatchBoundaryExclude, PatchBoundaryInclude, num_patches + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +u_exact(x) = VectorValue(-x[1],x[2]) +p_exact(x) = x[1] - 0.5 + +Dc = 2 +model = CartesianDiscreteModel((0,1,0,1),(2,2)) +labels = get_face_labeling(model) +add_tag_from_tags!(labels,"newman",[5]) +add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) + +PD_e = PatchDecomposition(model,patch_boundary_style=PatchBoundaryExclude()) +PD_i = PatchDecomposition(model,patch_boundary_style=PatchBoundaryInclude()) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +#Vh = TestFESpace(model,reffe_u) +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) + +PD = PD_i +topo = get_grid_topology(model) +#patches_mask = fill(false,num_patches(PD)) +patches_mask = get_isboundary_face(topo,0) + +Ph = PatchFESpace(Vh,PD,reffe_u;conformity=H1Conformity(),patches_mask) +Lh = PatchFESpace(Qh,PD,reffe_p;conformity=L2Conformity(),patches_mask) + +Xh = MultiFieldFESpace([Uh,Qh]) +Yh = MultiFieldFESpace([Vh,Qh]) +Zh = MultiFieldFESpace([Ph,Lh]) + +qdegree = 2*order +Ω = Triangulation(model) +Ωp = Triangulation(PD) + +dΩ = Measure(Ω,qdegree) +dΩp = Measure(Ωp,qdegree) + +Γ = BoundaryTriangulation(model,tags="newman") +dΓ = Measure(Γ,qdegree) +n = get_normal_vector(Γ) + +I_tensor = one(TensorValue{Dc,Dc,Float64}) +f(x) = -Δ(u_exact)(x) - ∇(p_exact)(x) +σ(x) = ∇(u_exact)(x) + p_exact(x)*I_tensor +biform((u,p),(v,q),dΩ) = ∫(∇(u)⊙∇(v) + (∇⋅v)*p - (∇⋅u)*q)dΩ +a(x,y) = biform(x,y,dΩ) +l((v,q)) = ∫(v⋅f)dΩ + ∫(v⋅(σ⋅n) )dΓ + +ap(x,y) = biform(x,y,dΩp) + +op = AffineFEOperator(a,l,Xh,Yh) +A = get_matrix(op) +b = get_vector(op) +cond(Matrix(A)) +x_exact = A\b + +uh_exact, ph_exact = FEFunction(Xh,x_exact) +l2_error(uh_exact,u_exact,dΩ) +l2_error(ph_exact,p_exact,dΩ) +l2_norm(∇⋅uh_exact,dΩ) + +P = PatchBasedSmoothers.PatchBasedLinearSolver(ap,Zh,Yh) +P_ns = numerical_setup(symbolic_setup(P,A),A) +x = zeros(num_free_dofs(Yh)) +r = b - A*x +solve!(x,P_ns,r) + +Ap = assemble_matrix(ap,Zh,Zh) + +solver = FGMRESSolver(10,P;rtol=1.e-8,verbose=true) +#solver = GMRESSolver(10,verbose=true) +ns = numerical_setup(symbolic_setup(solver,A),A) + +x = zeros(num_free_dofs(Yh)) +solve!(x,ns,b) + +norm(A*x - b) diff --git a/test/_dev/Vanka/stokesdarcy_vanka_GMG.jl b/test/_dev/Vanka/stokesdarcy_vanka_GMG.jl new file mode 100644 index 00000000..2a1aa183 --- /dev/null +++ b/test/_dev/Vanka/stokesdarcy_vanka_GMG.jl @@ -0,0 +1,172 @@ + +using Gridap +using GridapSolvers +using LinearAlgebra + +using Gridap.ReferenceFEs, Gridap.FESpaces, Gridap.Geometry +using Gridap.Algebra, Gridap.Arrays, Gridap.Adaptivity +using GridapSolvers.PatchBasedSmoothers, GridapSolvers.LinearSolvers + +function l2_norm(uh,dΩ) + sqrt(sum(∫(uh⋅uh)dΩ)) +end + +function l2_error(uh,u_exact,dΩ) + eh = uh - u_exact + return sqrt(sum(∫(eh⋅eh)dΩ)) +end + +function add_labels_2d(model) + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"newman",[5]) + add_tag_from_tags!(labels,"dirichlet",[collect(1:4)...,6,7,8]) +end + +function add_labels_3d(model) + labels = get_face_labeling(model) + add_tag_from_tags!(labels,"newman",[21]) + add_tag_from_tags!(labels,"dirichlet",[collect(1:20)...,22,23,24,25,26]) +end + +B = VectorValue(0.0,0.0,1.0) +u_exact(x) = VectorValue(x[1],-x[2],0.0) +p_exact(x) = sum(x) +j_exact(x) = VectorValue(-x[2],-x[1],0.0) + VectorValue(1.0,1.0,0.0) +φ_exact(x) = x[1] + x[2] + +Dc = 3 +domain = (Dc == 2) ? (0,1,0,1) : (0,1,0,1,0,1) +ncells = (Dc == 2) ? (4,4) : (2,2,2) +cmodel = CartesianDiscreteModel(domain,ncells) +if Dc == 2 + add_labels_2d(cmodel) +else + add_labels_3d(cmodel) +end +model = refine(cmodel,2) + +order = 2 +reffe_u = ReferenceFE(lagrangian,VectorValue{Dc,Float64},order) +reffe_p = ReferenceFE(lagrangian,Float64,order-1;space=:P) +reffe_j = ReferenceFE(raviart_thomas,Float64,order-1) +reffe_φ = ReferenceFE(lagrangian,Float64,order-1) + +VH = TestFESpace(cmodel,reffe_u;dirichlet_tags="dirichlet") +UH = TrialFESpace(VH,u_exact) +QH = TestFESpace(cmodel,reffe_p,conformity=:L2) +DH = TestFESpace(cmodel,reffe_j;dirichlet_tags="dirichlet") +JH = TrialFESpace(DH,j_exact) +ΦH = TestFESpace(cmodel,reffe_φ,conformity=:L2) + +XH = MultiFieldFESpace([UH,QH,JH,ΦH]) +YH = MultiFieldFESpace([VH,QH,DH,ΦH]) + +Vh = TestFESpace(model,reffe_u;dirichlet_tags="dirichlet") +Uh = TrialFESpace(Vh,u_exact) +Qh = TestFESpace(model,reffe_p,conformity=:L2) +Dh = TestFESpace(model,reffe_j;dirichlet_tags="dirichlet") +Jh = TrialFESpace(Dh,j_exact) +Φh = TestFESpace(model,reffe_φ,conformity=:L2) +Xh = MultiFieldFESpace([Uh,Qh,Jh,Φh]) +Yh = MultiFieldFESpace([Vh,Qh,Dh,Φh]) + +qdegree = 2*order +Ωh = Triangulation(model) +dΩh = Measure(Ωh,qdegree) +ΩH = Triangulation(cmodel) +dΩH = Measure(ΩH,qdegree) +dΩHh = Measure(ΩH,Ωh,qdegree) + +Γh = BoundaryTriangulation(model,tags="newman") +dΓh = Measure(Γh,qdegree) +n = get_normal_vector(Γh) + +α = -1.0 +I_tensor = one(TensorValue{Dc,Dc,Float64}) +f(x) = -Δ(u_exact)(x) - ∇(p_exact)(x) #- cross(j_exact(x),B) +g(x) = j_exact(x) - ∇(φ_exact)(x) #- cross(u_exact(x),B) +σ(x) = ∇(u_exact)(x) + p_exact(x)*I_tensor +γ(x) = φ_exact(x)*I_tensor +crossB(u,v,dΩ) = ∫(cross(u,B)⋅v)dΩ +graddiv(u,v,dΩ) = ∫((∇⋅u)⋅(∇⋅v))*dΩ +function a((u,p,j,φ),(v,q,d,ψ),dΩ) + c = ∫(∇(u)⊙∇(v) + (∇⋅v)*p - (∇⋅u)*q)dΩ + c = c + ∫(j⋅d + (∇⋅d)*φ - (∇⋅j)*ψ)dΩ +# c = c - crossB(u,d,dΩ) - crossB(j,v,dΩ) + if α > 0.0 + c += graddiv(u,v,dΩ) + graddiv(j,d,dΩ) + end + return c +end +l((v,q,d,ψ),dΩ,dΓ) = ∫(v⋅f + d⋅g)dΩ + ∫(v⋅(σ⋅n) + d⋅(γ⋅n))dΓ + +ah(x,y) = a(x,y,dΩh) +aH(x,y) = a(x,y,dΩH) +lh(y) = l(y,dΩh,dΓh) + +op_h = AffineFEOperator(ah,lh,Xh,Yh) +Ah = get_matrix(op_h) +bh = get_vector(op_h) +xh_exact = Ah\bh + +AH = assemble_matrix(aH,XH,YH) +Mhh = assemble_matrix((u,v)->∫(u⋅v)*dΩh,Xh,Xh) + +uh_exact, ph_exact, jh_exact, φh_exact = FEFunction(Xh,xh_exact) +l2_error(uh_exact,u_exact,dΩh) +l2_error(ph_exact,p_exact,dΩh) +l2_error(jh_exact,j_exact,dΩh) +l2_error(φh_exact,φ_exact,dΩh) +l2_norm(∇⋅uh_exact,dΩh) +l2_norm(∇⋅jh_exact,dΩh) + +PD = PatchDecomposition(model) +smoother = RichardsonSmoother(VankaSolver(Xh,PD),10,0.05) +smoother_ns = numerical_setup(symbolic_setup(smoother,Ah),Ah) + +function project_f2c(rh) + Qrh = Mhh\rh + uh, ph, jh, φh = FEFunction(Yh,Qrh) + ll((v,q,d,ψ)) = ∫(v⋅uh + q*ph + d⋅jh + ψ⋅φh)*dΩHh + assemble_vector(ll,YH) +end +function interp_c2f(xH) + get_free_dof_values(interpolate(FEFunction(YH,xH),Yh)) +end + +xh = randn(size(Ah,2)) +rh = bh - Ah*xh +niters = 10 + +iter = 0 +error0 = norm(rh) +error = error0 +e_rel = error/error0 +while iter < niters && e_rel > 1.0e-10 + println("Iter $iter:") + println(" > Initial: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + println(" > Pre-smoother: ", norm(rh)) + + rH = project_f2c(rh) + qH = AH\rH + qh = interp_c2f(qH) + + rh = rh - Ah*qh + xh = xh + qh + println(" > Post-correction: ", norm(rh)) + + solve!(xh,smoother_ns,rh) + + iter += 1 + error = norm(rh) + e_rel = error/error0 + println(" > Final: ",error, " - ", e_rel) +end + +uh, ph, jh, φh = FEFunction(Xh,xh) +l2_error(uh,u_exact,dΩh) +l2_error(ph,p_exact,dΩh) +l2_error(jh,j_exact,dΩh) +l2_error(φh,φ_exact,dΩh)