Forum Message

 

 

We have moved the forum to https://forum.ngsolve.org . 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.

Notice

The forum is in read only mode.

Stokes problem (pressure order convergence)

More
1 year 9 months ago - 1 year 9 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,
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) 
Last edit: 1 year 9 months ago by martinebert. Reason: formating code
More
1 year 9 months ago #4426 by THaubold
Hi,

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))

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