From c53babbe68b57f673e8ab95aed6168c468689d16 Mon Sep 17 00:00:00 2001 From: amartin Date: Sun, 22 May 2022 18:22:59 +1000 Subject: [PATCH] * Added convergence analysis to HHO driver * Added stabilization --- test/PoissonHHOTests.jl | 148 ++++++++++++++++++--------- test/StressAssistedHenckyHDGTests.jl | 40 ++------ 2 files changed, 107 insertions(+), 81 deletions(-) diff --git a/test/PoissonHHOTests.jl b/test/PoissonHHOTests.jl index f5ed56d..0f33423 100644 --- a/test/PoissonHHOTests.jl +++ b/test/PoissonHHOTests.jl @@ -11,6 +11,8 @@ module PoissonHHOTests reffe_c = ReferenceFE(monomial_basis , Float64, order+1; subspace=:OnlyConstant) reffe_nc = ReferenceFE(monomial_basis , Float64, order+1; subspace=:ExcludeConstant) + Ω = dΩ.quad.trian + VKR = TestFESpace(Ω , refferecᵤ; conformity=:L2) UKR = TrialFESpace(VKR) UKR_NZM = TrialFESpace(TestFESpace(Ω, reffe_nzm; conformity=:L2)) @@ -42,57 +44,107 @@ module PoissonHHOTests LocalFEOperator((m,n),UK_U∂K,VK_V∂K; field_type_at_common_faces=MultiValued()) end - u(x)=x[1]+x[2] - f(x)=-Δ(u)(x) - - model=CartesianDiscreteModel((0,1,0,1),(2,2)) - D = num_cell_dims(model) - Ω = Triangulation(ReferenceFE{D},model) - Γ = Triangulation(ReferenceFE{D-1},model) - ∂K = GridapHybrid.Skeleton(model) - - order = 1 - reffeᵤ = ReferenceFE(lagrangian,Float64,order ;space=:P) - - VK = TestFESpace(Ω , reffeᵤ; conformity=:L2) - V∂K = TestFESpace(Γ , reffeᵤ; conformity=:L2,dirichlet_tags=collect(5:8)) - UK = TrialFESpace(VK) - U∂K = TrialFESpace(V∂K,u) - VK_V∂K = MultiFieldFESpace([VK,V∂K]) - UK_U∂K = MultiFieldFESpace([UK,U∂K]) - - degree = 2*order+1 - dΩ = Measure(Ω,degree) - dΓ = Measure(Γ,degree) - d∂K = Measure(∂K,degree) - - R=setup_reconstruction_operator(model, order, dΩ, d∂K) - diff_op=setup_difference_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K) - - function r(u,v) - uK,u∂K=R(u) - vK,v∂K=R(v) - ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ + p = 3 + u(x) = x[1]^p+x[2]^p + f(x)=-Δ(u)(x) + + #u(x)=x[1]+x[2] + #f(x)=-Δ(u)(x) + + function solve_hho(cells,order) + partition = (0,1,0,1) + model = CartesianDiscreteModel(partition, cells) + D = num_cell_dims(model) + Ω = Triangulation(ReferenceFE{D},model) + Γ = Triangulation(ReferenceFE{D-1},model) + ∂K = GridapHybrid.Skeleton(model) + + reffeᵤ = ReferenceFE(lagrangian,Float64,order ;space=:P) + + VK = TestFESpace(Ω , reffeᵤ; conformity=:L2) + V∂K = TestFESpace(Γ , reffeᵤ; conformity=:L2,dirichlet_tags=collect(5:8)) + UK = TrialFESpace(VK) + U∂K = TrialFESpace(V∂K,u) + VK_V∂K = MultiFieldFESpace([VK,V∂K]) + UK_U∂K = MultiFieldFESpace([UK,U∂K]) + + degree = 2*order+1 + dΩ = Measure(Ω,degree) + dΓ = Measure(Γ,degree) + d∂K = Measure(∂K,degree) + + R=setup_reconstruction_operator(model, order, dΩ, d∂K) + diff_op=setup_difference_operator(UK_U∂K,VK_V∂K,R,dΩ,d∂K) + + function r(u,v) + uK,u∂K=R(u) + vK,v∂K=R(v) + ∫(∇(vK)⋅∇(uK))dΩ + ∫(∇(vK)⋅∇(u∂K))dΩ + + ∫(∇(v∂K)⋅∇(uK))dΩ + ∫(∇(v∂K)⋅∇(u∂K))dΩ + end + + function s(u,v) + h_T=CellField(get_array(∫(1)dΩ),Ω) + h_T_1=1.0/h_T + h_T_2=1.0/(h_T*h_T) + + δuK,δu∂K=diff_op(u) + δvK,δv∂K=diff_op(v) + δvK_K,δvK_∂K=δvK + δuK_K,δuK_∂K=δuK + δv∂K_K,δv∂K_∂K=δv∂K + δu∂K_K,δu∂K_∂K=δu∂K + ∫(h_T_2*δvK_K*δuK_K)dΩ+ + ∫(h_T_2*δvK_K*δuK_∂K)dΩ+ + ∫(h_T_2*δvK_∂K*δuK_K)dΩ+ + ∫(h_T_2*δvK_∂K*δuK_∂K)dΩ+ + ∫(h_T_1*δv∂K_∂K*δu∂K_∂K)d∂K + end + + a(u,v)=r(u,v)+s(u,v) + l((vK,))=∫(vK*f)dΩ + + + @time op=HybridAffineFEOperator((u,v)->(a(u,v),l(v)), UK_U∂K, VK_V∂K, [1], [2]) + @time xh=solve(op) + + uhK,uh∂K=xh + e = u -uhK + # @test sqrt(sum(∫(e⋅e)dΩ)) < 1.0e-12 + return sqrt(sum(∫(e⋅e)dΩ)) end - function s(u,v) - δuK,δu∂K=diff_op(u) - δvK,δv∂K=diff_op(v) - δvK_K,δvK_∂K=δvK - δuK_K,δuK_∂K=δuK - δv∂K_K,δv∂K_∂K=δv∂K - δu∂K_K,δu∂K_∂K=δu∂K - ∫(δvK_K*δuK_K)dΩ+∫(δvK_K*δuK_∂K)dΩ+∫(δvK_∂K*δuK_K)dΩ+∫(δvK_∂K*δuK_∂K)dΩ+ - ∫(δv∂K_∂K*δu∂K_∂K)d∂K + function conv_test(ns,order) + el2 = Float64[] + hs = Float64[] + for n in ns + l2 = solve_hho((n,n),order) + println(l2) + h = 1.0/n + push!(el2,l2) + push!(hs,h) + end + println(el2) + el2 end - a(u,v)=r(u,v)+s(u,v) - l((vK,))=∫(vK*f)dΩ - - op=HybridAffineFEOperator((u,v)->(a(u,v),l(v)), UK_U∂K, VK_V∂K, [1], [2]) - xh=solve(op) + function slope(hs,errors) + x = log10.(hs) + y = log10.(errors) + linreg = hcat(fill!(similar(x), 1), x) \ y + linreg[2] + end - uhK,uh∂K=xh - e = u -uhK - @test sqrt(sum(∫(e⋅e)dΩ)) < 1.0e-12 + ns=[8,16,32,64,128] + order=1 + el, hs = conv_test(ns,order) + println("Slope L2-norm u: $(slope(hs,el))") + slopek =[Float64(ni)^(-(order)) for ni in ns] + slopekp1=[Float64(ni)^(-(order+1)) for ni in ns] + slopekp2=[Float64(ni)^(-(order+2)) for ni in ns] + plot(hs,[el slopek slopekp1 slopekp2], + xaxis=:log, yaxis=:log, + label=["L2u (measured)" "slope k" "slope k+1" "slope k+2"], + shape=:auto, + xlabel="h",ylabel="L2 error",legend=:bottomright) end diff --git a/test/StressAssistedHenckyHDGTests.jl b/test/StressAssistedHenckyHDGTests.jl index a2ad7e1..83a28b4 100644 --- a/test/StressAssistedHenckyHDGTests.jl +++ b/test/StressAssistedHenckyHDGTests.jl @@ -165,31 +165,16 @@ function solve_stress_assisted_diffusion_hencky_hdg(cells,order;write_results=fa MuΓ_TR = TrialFESpace(MuΓ,u₀) Y_TR = MultiFieldFESpace([W_TR, W_TR, Ψ_TR, S_TR, S_TR, V_TR, MϕΓ_TR, MuΓ_TR]) - # FE formulation params - degree = 2 * (order + 1) # TO-DO: To think which is the minimum degree required - dΩ = Measure(Ω, degree) - n = get_cell_normal_vector(∂K) + # Integration + degree = 2 * (order + 1) + dΩ = Measure(Ω, degree) d∂K = Measure(∂K, degree) - # Stabilization operator (solute concentration) - τρΓ = 1.0 # To be defined??? - # Stabilization operator (displacements) - τuΓ = Float64(cells[1]) # To be defined??? - - println("h=$(1.0/cells[1]) τuΓ=$(τuΓ) τρΓ=$(τρΓ)") - - Y_basis = get_fe_basis(Y) - Y_TR_basis = get_trial_fe_basis(Y_TR) - (_,_,_,_,_,_,_,uhΓ_basis) = Y_TR_basis + # Outer normal, cell boundary + n = get_cell_normal_vector(∂K) - # # Testing for correctness of Pₘ projection - # xh = FEFunction(Y_TR,rand(num_free_dofs(Y_TR))) - # _,_,_,_,_,_,_,μ = Y_basis - # _, uh, _, _ = xh - # dc=∫(μ⋅(uh-Pₘ(uh, uhΓ_basis, μ, d∂K)))d∂K - # for k in dc.dict[d∂K.quad.trian] - # @test all(k[8] .< 1.0e-14) - # end + # Stabilization operator (displacements) + τuΓ = Float64(cells[1]) # Residual # current_iterate => (ωh,ρh,ϕh,th,σh,uh,ϕhΓ,uhΓ) @@ -207,16 +192,7 @@ function solve_stress_assisted_diffusion_hencky_hdg(cells,order;write_results=fa ∫( μh⋅(σh⋅n) )d∂K - ∫( τuΓ*μh⋅Pm_uh)d∂K + ∫( τuΓ*μh⋅uhΓ )d∂K end - # free_values = zeros(num_free_dofs(Y_TR)) - # xh = FEFunction(Y_TR,free_values) - # (ωh,ρh,ϕh,th,σh,uh,ϕhΓ,uhΓ) = xh - # (rh,nh,ψh,sh,τh,vh,ηh,μh) = Y_basis - # dc=residual((ωh,ρh,ϕh,th,σh,uh,ϕhΓ,uhΓ),(rh,nh,ψh,sh,τh,vh,ηh,μh)) - # res=assemble_vector(dc,Y) - # dc=jacobian(x->residual(x,Y_basis),xh) - # A=assemble_matrix(dc,Y_TR,Y) @time op=HybridFEOperator(residual,Y_TR,Y,collect(1:6),[7,8]) - # op=FEOperator(residual,Y_TR,Y) nls = NLSolver(show_trace=true, method=:newton, ftol=1.0e-14) solver = FESolver(nls) @@ -259,8 +235,6 @@ function solve_stress_assisted_diffusion_hencky_hdg(cells,order;write_results=fa end norm2_ω,norm2_ρ,norm2_t,norm2_ϕ,norm2_σ,norm2_u - #@test sqrt(sum(∫(eu ⋅ eu)dΩ)) < 1.0e-12 - #@test sqrt(sum(∫(eσ ⊙ eσ)dΩ)) < 1.0e-12 end