diff --git a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl index d00bf98..d836dcf 100644 --- a/src/PatchBasedSmoothers/PatchBasedSmoothers.jl +++ b/src/PatchBasedSmoothers/PatchBasedSmoothers.jl @@ -13,7 +13,7 @@ using GridapSolvers.MultilevelTools export PatchDecomposition export PatchFESpace -export PatchBasedLinearSolver +export PatchBasedLinearSolver, VankaSolver export PatchProlongationOperator, PatchRestrictionOperator export setup_patch_prolongation_operators, setup_patch_restriction_operators @@ -26,10 +26,12 @@ include("seq/PatchTriangulations.jl") # FESpaces include("seq/PatchFESpaces.jl") include("mpi/PatchFESpaces.jl") +include("seq/ZeroMeanPatchFESpaces.jl") include("seq/PatchMultiFieldFESpaces.jl") # Solvers include("seq/PatchBasedLinearSolvers.jl") include("seq/PatchTransferOperators.jl") +include("seq/VankaSolvers.jl") end \ No newline at end of file diff --git a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl b/src/PatchBasedSmoothers/seq/PatchFESpaces.jl index d019307..930dd77 100644 --- a/src/PatchBasedSmoothers/seq/PatchFESpaces.jl +++ b/src/PatchBasedSmoothers/seq/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 """ @@ -103,13 +115,13 @@ function PatchFESpace( 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 diff --git a/src/PatchBasedSmoothers/seq/VankaSolvers.jl b/src/PatchBasedSmoothers/seq/VankaSolvers.jl new file mode 100644 index 0000000..c819cdf --- /dev/null +++ b/src/PatchBasedSmoothers/seq/VankaSolvers.jl @@ -0,0 +1,79 @@ + +struct VankaSolver{A} <: Algebra.LinearSolver + patch_ids::A + function VankaSolver( + patch_ids::AbstractVector{<:AbstractVector{<:Integer}} + ) + A = typeof(patch_ids) + new{A}(patch_ids) + end +end + +function VankaSolver(space::MultiFieldFESpace) + 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::MultiFieldFESpace,patch_decomposition::PatchDecomposition) + patch_cells = patch_decomposition.patch_cells + return VankaSolver(space,patch_cells) +end + +function VankaSolver(space::MultiFieldFESpace,patch_cells::Table{<:Integer}) + cell_ids = get_cell_dof_ids(space) + patch_ids = map(patch_cells) do cells + ids = vcat([vcat(cell_ids[cell].array...) for cell in cells]...) + filter!(x->x>0,ids) + sort!(ids) + unique!(ids) + return ids + end + return VankaSolver(patch_ids) +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/seq/ZeroMeanPatchFESpaces.jl b/src/PatchBasedSmoothers/seq/ZeroMeanPatchFESpaces.jl new file mode 100644 index 0000000..0785dd9 --- /dev/null +++ b/src/PatchBasedSmoothers/seq/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/test/_dev/MacroFEs/block_stokes.jl b/test/_dev/MacroFEs/block_stokes.jl new file mode 100644 index 0000000..10dae06 --- /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/stokes.jl b/test/_dev/MacroFEs/stokes.jl index ea5887b..d575fc5 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/PatchBased/ZeroMeanPatchBasedFESpaces.jl b/test/_dev/PatchBased/ZeroMeanPatchBasedFESpaces.jl new file mode 100644 index 0000000..5b86a5d --- /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 = PD.patch_cells_overlapped +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/Vanka/darcy_vanka.jl b/test/_dev/Vanka/darcy_vanka.jl new file mode 100644 index 0000000..bfd185f --- /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 0000000..cbd7eb3 --- /dev/null +++ b/test/_dev/Vanka/darcy_vanka_GMG.jl @@ -0,0 +1,132 @@ +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 !iszero(α) + 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),5,0.2) +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/stokes_vanka.jl b/test/_dev/Vanka/stokes_vanka.jl new file mode 100644 index 0000000..67f9a5b --- /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 0000000..23cbc97 --- /dev/null +++ b/test/_dev/Vanka/stokes_vanka_GMG.jl @@ -0,0 +1,125 @@ +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) + +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),dΩ) = ∫(∇(u)⊙∇(v) + (∇⋅v)*p - (∇⋅u)*q)dΩ +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),5,0.2) +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/stokes_vanka_PCont.jl b/test/_dev/Vanka/stokes_vanka_PCont.jl new file mode 100644 index 0000000..75da9b4 --- /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 0000000..5a2e7b8 --- /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)