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,
Martin
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()
geo.AddRectangle((0,0),(1,1),bc="rectangle")
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)))
mesh.Curve(order=order)
print("ElementBoundary=",mesh.GetBoundaries())
print("NumEles=",mesh.ne)
print("NumEdges=",mesh.nedge)
print("NumVertex=",mesh.nv)
# 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)
a.Assemble()
# assembling the rhs
f = LinearForm(X)
source = CoefficientFunction(src)
f += SymbolicLFI(InnerProduct(source,v))
f.Assemble()
#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()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += 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)