Forum Message



We have moved the forum to . This is an archived version of the topics until 05/05/23. All the topics were moved to the new forum and conversations can be continued there. This forum is just kept as legacy to not invalidate old links. If you want to continue a conversation just look for the topic in the new forum.


The forum is in read only mode.

Stokes problem (pressure order convergence)

1 year 3 months ago - 1 year 3 months ago #4425 by martinebert
Dear developer,
I am new to NGSolve and trying to solve the Stokes problem to check the velocity and pressure convergence order for the known solution. The velocity error converges to the right order but the pressure doesn't converge. Anyone can help me to get the right convergence orders please? Thanks
Best regards,

The code:

from ngsolve import *
# import netgen.gui
from netgen.geom2d import unit_square
from ngsolve.webgui import Draw
from math import pi
from netgen.geom2d import SplineGeometry
import numpy as np

geo = SplineGeometry()

nu = 1.0
# exact velocity
ex_vel = (sin(2*pi*x) * sin(2*pi*y),
          cos(2*pi*x) * cos(2*pi*y))
# exact pressure
ex_pre = 0.25*(cos(4*pi*x) - cos(4*pi*y))
# source term
src = ( 8*pi*pi*sin(2*pi*x)*sin(2*pi*y) + pi*sin(4*pi*x),
       8*pi*pi*cos(2*pi*y)*cos(2*pi*x) + pi*sin(4*pi*y))

mesh = Mesh(geo.GenerateMesh(maxh=1/(2*N)))
# finite element formulation
V = VectorH1(mesh,order=order, dirichlet="rectangle")
Q = H1(mesh,order=order-1)
X = FESpace([V,Q])

u,p = X.TrialFunction()
v,q = X.TestFunction()

print("velocity dofs ", V.ndof, X.Range(0))
print("pressure dofs ", Q.ndof, X.Range(1))
print("Total dofs ", X.ndof)
# assembling the system matrix
stokes = nu*InnerProduct(grad(u), grad(v))-div(u)*q-div(v)*p - 1e-10*p*q
a = BilinearForm(X)
a += SymbolicBFI(stokes)
# assembling the rhs
f = LinearForm(X)
source = CoefficientFunction(src)
f += SymbolicLFI(InnerProduct(source,v))
#Exact Dirichelt BC for velocity
gfu = GridFunction(X)

u_exact = CoefficientFunction(ex_vel)
gfu.components[0].Set(u_exact, definedon=mesh.Boundaries("rectangle"))    
# solve the linear system
inv_stokes = a.mat.Inverse(X.FreeDofs())

res = f.vec.CreateVector() = f.vec - a.mat*gfu.vec += inv_stokes * res

velocity = gfu.components[0]
pressure = gfu.components[1]
# Draw( velocity, mesh, "velocity" )
# Draw( Norm(velocity), mesh, "absvel(hdiv)")
# Draw( pressure, mesh, "pressure" )
# computing the error estimates
p_exact = CoefficientFunction(ex_pre)
pressure = CoefficientFunction(gfu.components[1])
errp = sqrt(Integrate((p_exact-pressure)*(p_exact-pressure), mesh))

velocity = CoefficientFunction(gfu.components[0])    
erru = sqrt (Integrate ( (u_exact-velocity )*(u_exact-velocity), mesh))

print("L2(u):", erru)
print("L2(p):", errp) 
Last edit: 1 year 3 months ago by martinebert. Reason: formating code
1 year 3 months ago #4426 by THaubold

your source term has a small typo. The sign in the pressure part should be negative, i.e.

src = ( 8*pi*pi*sin(2*pi*x)*sin(2*pi*y) - pi*sin(4*pi*x),
       8*pi*pi*cos(2*pi*y)*cos(2*pi*x) + pi*sin(4*pi*y))

1 year 3 months ago #4427 by martinebert
Dear Tim,
thank you very much. It was a stupid mistake but I appreciate your help.
Best regards,
Time to create page: 0.124 seconds