Skip to content

Commit

Permalink
Restore Stokes tests after devs
Browse files Browse the repository at this point in the history
  • Loading branch information
ericneiva committed Oct 7, 2024
1 parent 8d0476b commit 76230e5
Show file tree
Hide file tree
Showing 2 changed files with 64 additions and 84 deletions.
137 changes: 62 additions & 75 deletions test/DistributedTests/StokesTests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ p(x) = x[1] - x[2]
f(x) = - Δ(u)(x) + (p)(x)
g(x) = (∇u)(x)

function main(distribute,parts;n=4,cells=(n,n),order=2,name="a")
function main(distribute,parts;n=4,cells=(n,n),order=2)

@show parts
ranks = distribute(LinearIndices((prod(parts),)))
Expand All @@ -36,92 +36,79 @@ function main(distribute,parts;n=4,cells=(n,n),order=2,name="a")
degree = order == 1 ? 3 : 2*order
vquad = Quadrature(algoim,phi,degree,phase=IN)
v_cell_quad,is_a = CellQuadratureAndActiveMask(model₀,vquad)
# squad = Quadrature(algoim,phi,degree)
# s_cell_quad,is_c = CellQuadratureAndActiveMask(model₀,squad)
squad = Quadrature(algoim,phi,degree)
s_cell_quad,is_c = CellQuadratureAndActiveMask(model₀,squad)

# model,aggregates = aggregate(model₀,is_a,is_c,IN)
model = model₀
# map(aggregates) do agg
# @show agg
# end
model,aggregates = aggregate(model₀,is_a,is_c,IN)

Ω = Triangulation(model)
Ωᵃ,dΩᵃ = TriangulationAndMeasure(Ω,v_cell_quad,is_a)
# map(local_views(dΩᵃ)) do dΩᵃ
# @show dΩᵃ.quad.cell_point.ptrs
# end
# map(local_views(dΩᵃ)) do dΩᵃ
# @show dΩᵃ.quad.cell_weight.ptrs
# end

# _,dΓ = TriangulationAndMeasure(Ω,s_cell_quad,is_c)
# n_Γ = normal(phi,Ω)
_,dΓ = TriangulationAndMeasure(Ω,s_cell_quad,is_c)
n_Γ = normal(phi,Ω)

dbg = parts[1]
writevtk(dΓ,"cut_quad_$dbg")

@show (( ( 1.0 )dΩᵃ ))
# @show √(∑( ∫( 1.0 )dΩᵃ ))
# @show √(∑( ∫( 1.0 )dΓ ))

# return

# # reffeᵘ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
# # reffeᵖ = ReferenceFE(lagrangian,Float64,order-1,space=:P)
reffeᵘ = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffeᵖ = ReferenceFE(lagrangian,Float64,order-1,space=:P)

# # Vstdᵘ = TestFESpace(Ωᵃ,reffeᵘ)
# # Vstdᵖ = TestFESpace(Ωᵃ,reffeᵖ)
Vstdᵘ = TestFESpace(Ωᵃ,reffeᵘ)
Vstdᵖ = TestFESpace(Ωᵃ,reffeᵖ)

# # Vᵘ = AgFEMSpace(model,Vstdᵘ,aggregates)
# # Vᵖ = AgFEMSpace(model,Vstdᵖ,aggregates)
# # Vˡ = ConstantFESpace(model)
Vᵘ = AgFEMSpace(model,Vstdᵘ,aggregates)
Vᵖ = AgFEMSpace(model,Vstdᵖ,aggregates)
= ConstantFESpace(model)

# # Uᵘ = TrialFESpace(Vᵘ)
# # Uᵖ = TrialFESpace(Vᵖ)
# # Uˡ = TrialFESpace(Vˡ)
Uᵘ = TrialFESpace(Vᵘ)
Uᵖ = TrialFESpace(Vᵖ)
= TrialFESpace(Vˡ)

# # Y = MultiFieldFESpace([Vᵘ,Vᵖ,Vˡ])
# # X = MultiFieldFESpace([Uᵘ,Uᵖ,Uˡ])

# # Nitsche method
# γᵈ = 2.0*order^2
# h = (domain[2]-domain[1])/cells[1]

# a((u,p,l),(v,q,ℓ)) =
# ∫( ∇(u)⊙∇(v) - q*(∇⋅u) - p*(∇⋅v) )dΩᵃ +
# ∫( p*ℓ )dΩᵃ + ∫( q*l )dΩᵃ +
# ∫( (γᵈ/h)*(u⋅v) -
# v⋅(n_Γ⋅∇(u)) - u⋅(n_Γ⋅∇(v)) +
# p*(n_Γ⋅v) + q*(n_Γ⋅u) )dΓ

# l((v,q,ℓ)) =
# ∫( v⋅f - q*g )dΩᵃ +
# ∫( (γᵈ/h)*(u⋅v) - u⋅(n_Γ⋅∇(v)) + (q*n_Γ)⋅u )dΓ

# # FE problem
# # op = AffineFEOperator(a,l,X,Y)
# # uₕ,pₕ,lₕ = solve(op)
# # writevtk(Ωᵃ,"stokes",cellfields=["u"=>uₕ,"p"=>pₕ])

# # uₕ = zero(Uᵘ)
# # pₕ = zero(Uᵖ)

# euₕ = CellField(u,Ωᵃ) # - uₕ
# epₕ = CellField(p,Ωᵃ) # - pₕ

# writevtk(Ωᵃ,"stokes_$name",cellfields=["u"=>euₕ,"p"=>epₕ])

# l2(e) = √(∑( ∫( e⋅e )dΩᵃ ))
# h1(e) = √(∑( ∫( e⋅e + ∇(e)⊙∇(e) )dΩᵃ ))

# el2ᵘ = l2(euₕ)
# eh1ᵘ = h1(euₕ)
# el2ᵖ = l2(epₕ)
# # el2ˡ = l2(lₕ)

# # @test el2ᵘ < 1.0e-014
# # @test eh1ᵘ < 1.0e-013
# # @test el2ᵖ < 1.0e-014
# @show el2ᵘ
# @show eh1ᵘ
# @show el2ᵖ
# # @show el2ˡ
Y = MultiFieldFESpace([Vᵘ,Vᵖ,Vˡ])
X = MultiFieldFESpace([Uᵘ,Uᵖ,Uˡ])

# Nitsche method
γᵈ = 2.0*order^2
h = (domain[2]-domain[1])/cells[1]

a((u,p,l),(v,q,ℓ)) =
( (u)(v) - q*(∇u) - p*(∇v) )dΩᵃ +
( p*ℓ )dΩᵃ + ( q*l )dΩᵃ +
( (γᵈ/h)*(uv) -
v(n_Γ(u)) - u(n_Γ(v)) +
p*(n_Γv) + q*(n_Γu) )dΓ

l((v,q,ℓ)) =
( vf - q*g )dΩᵃ +
( (γᵈ/h)*(uv) - u(n_Γ(v)) + (q*n_Γ)u )dΓ

# FE problem
op = AffineFEOperator(a,l,X,Y)
uₕ,pₕ,lₕ = solve(op)

writevtk(Ωᵃ,"stokes",cellfields=["u"=>uₕ,"p"=>pₕ])

euₕ = u - uₕ
epₕ = p - pₕ

l2(e) = (( ( ee )dΩᵃ ))
h1(e) = (( ( ee + (e)(e) )dΩᵃ ))

el2ᵘ = l2(euₕ)
eh1ᵘ = h1(euₕ)
el2ᵖ = l2(epₕ)
# el2ˡ = l2(lₕ)

# @test el2ᵘ < 1.0e-014
# @test eh1ᵘ < 1.0e-013
# @test el2ᵖ < 1.0e-014
@show el2ᵘ
@show eh1ᵘ
@show el2ᵖ
# @show el2ˡ

end

Expand Down
11 changes: 2 additions & 9 deletions test/DistributedTests/sequential/StokesTests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,7 @@ module StokesTestsSeq
using PartitionedArrays
include("../StokesTests.jl")
with_debug() do distribute
StokesTests.main(distribute,(1,1),cells=(2,2),name="a")
# StokesTests.main(distribute,(1,2),cells=(2,2),name="b")
# StokesTests.main(distribute,(2,1),cells=(2,2),name="c")
StokesTests.main(distribute,(2,2),cells=(2,2),name="d")
# StokesTests.main_zeromean(distribute,(1,1),cells=(4,4))
# StokesTests.main_zeromean(distribute,(1,2),cells=(4,4))
# StokesTests.main_zeromean(distribute,(2,1),cells=(4,4))
# StokesTests.main_zeromean(distribute,(2,2),cells=(4,4))
# StokesTests.main(distribute,(4,4),cells=(16,16))
StokesTests.main(distribute,(1,1),cells=(16,16))
StokesTests.main(distribute,(2,2),cells=(16,16))
end
end

0 comments on commit 76230e5

Please sign in to comment.