Solution blow-up for Stokes problem

[ConnorMallon]

Refining the mesh in test/GridapTests/StokesTaylorHoodTests.jl ( to e.g. partition = (200,200) ) results in some unexpected behaviour. The solution blows up at some locations on the boundary:

The driver reproducing the results:

module StokesTaylorHoodTests

using Test
using Gridap
import Gridap: ∇

u(x) = VectorValue( x[1]^2 + 2*x[2]^2, -x[1]^2 )
p(x) = x[1] + 3*x[2]
f(x) = -Δ(u)(x) + ∇(p)(x)
g(x) = (∇⋅u)(x)
∇u(x) = ∇(u)(x)

domain = (0,2,0,2)
partition = (200,200)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"dirichlet",[1,2,5])
add_tag_from_tags!(labels,"neumann",[6,7,8])

order = 2

reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffe_p = ReferenceFE(lagrangian,Float64,order-1)

V = TestFESpace(model,reffe_u,labels=labels,dirichlet_tags="dirichlet",conformity=:H1)
Q = TestFESpace(model,reffe_p,labels=labels,conformity=:H1)

U = TrialFESpace(V,u)
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model,labels,tags="neumann")
n_Γ = get_normal_vector(Γ)

degree = order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)

a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )*dΩ

l((v,q)) = ∫( v⋅f + q*g )*dΩ + ∫( v⋅(n_Γ⋅∇u) - (n_Γ⋅v)*p )*dΓ

op = AffineFEOperator(a,l,X,Y)

uh, ph = solve(op)
writevtk(Ω,"uh",cellfields=["uh"=>uh])

eu = u - uh
ep = p - ph

l2(u) = sqrt(sum( ∫( u⊙u )*dΩ ))
h1(u) = sqrt(sum( ∫( u⊙u + ∇(u)⊙∇(u) )*dΩ ))

eu_l2 = l2(eu)
eu_h1 = h1(eu)
ep_l2 = l2(ep)

tol = 1.0e-9
@test eu_l2 < tol
@test eu_h1 < tol
@test ep_l2 < tol

end # module

[fverdugo]

Hi @ConnorMallon can you try with larger values of the integration degree degree ?

[ConnorMallon]

Hey @fverdugo! The problem still remains after increasing degree.

[ConnorMallon]

Hey @oriolcg , have you ever had a problem similar to this when solving for flow problems in a channel ? If I solve the Stokes problem in a channel I get the same solution blowup on the Neumann boundary when using a fine mesh:

module StokesChannel

using Test
using Gridap
import Gridap: ∇

u_in((x,y)) = VectorValue(y*(1-y),0) 
u_s(x) = VectorValue(0,0)
hN(x) = VectorValue(0,0)

domain = (0,4,0,1)
partition = (400,100)
model = CartesianDiscreteModel(domain,partition)

labels = get_face_labeling(model)
add_tag_from_tags!(labels,"inlet",[1,3,7])
add_tag_from_tags!(labels,"no_slip",[2,4,5,6,])
add_tag_from_tags!(labels,"outlet",[8])

order = 2

reffe_u = ReferenceFE(lagrangian,VectorValue{2,Float64},order)
reffe_p = ReferenceFE(lagrangian,Float64,order-1)

V = TestFESpace(model,reffe_u,labels=labels,dirichlet_tags=["inlet","no_slip"],conformity=:H1)
Q = TestFESpace(model,reffe_p,labels=labels,conformity=:H1)

U = TrialFESpace(V,[u_in,u_s])
P = TrialFESpace(Q)

Y = MultiFieldFESpace([V,Q])
X = MultiFieldFESpace([U,P])

Ω = Triangulation(model)
Γ = BoundaryTriangulation(model,labels,tags="outlet")
n_Γ = get_normal_vector(Γ)

degree = order
dΩ = Measure(Ω,degree)
dΓ = Measure(Γ,degree)

a((u,p),(v,q)) = ∫( ∇(v)⊙∇(u) - (∇⋅v)*p + q*(∇⋅u) )*dΩ

l((v,q)) = ∫( v⋅hN )*dΓ

op = AffineFEOperator(a,l,X,Y)

uh, ph = solve(op)

writevtk(Ω,"stokes_sol",cellfields=["uh"=>uh,"ph"=>ph])

end # module

Not that the blowup doesn't occur for courser meshes. The problem goes away with P0 elements for the pressure.

[ConnorMallon]

problem is related to the linear solver. change the ls or add terms in the bottom right block to avoid the issue (e.g. h^2 (∇p ⋅ ∇q) )