Problems with defining pressure BC for Stokes equation

Hy everybody,I try to solve the Stokes equation. If I define a velocity BC (U on surface ‘bot’) everything works fine and the expected result and the calculated result correspond. But, if I define a pressure BC, the results do not correspond to my expectations. On the surface ‘left’ I defined a pressure of 2 MPa and on the right surface I defined a pressure of 0.2 MPa. The velocity on the other surfaces is 0. My problem is, that after solving the Stokes equation, the pressure field on the left and right surface, does not have the value of the BC. So also, the velocity field is quite strange.So, maybe somebody can help me to define the pressure BC correctly.Many thanks to all of you.

Followed you can find my code: 
# V = Periodic(VectorH1(mesh,order=3, dirichlet='bot|top'))  # velocity
Q = H1(mesh,order=2, dirichlet='left|right') # pressure

X = V * Q 

# obtain test and trial functions
u, p = X.TrialFunction()
v, q = X.TestFunction()
# space for solution data
gfu = GridFunction(X)
# boundary condition
#uin = CoefficientFunction((U,0,0))
#gfu.components[0].Set(uin, definedon=mesh.Boundaries("bot"))
cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right"  else 0 for bc in mesh.GetBoundaries()])
    # Set BC
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))



nu = (eta/rho)  # kinematic viscosity

SetNumThreads(24) #Anzahl der Kerne zum Parallelisieren
with TaskManager():
    # right hand side of the equation
    f = LinearForm(X)
    f.Assemble()

    # left hand side of the equation
    stokes = nu*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
    a = BilinearForm(X, symmetric=True)
    a += stokes*dx
    a.Assemble()

    # calculate solution data
    inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
    res = f.vec.CreateVector()
    res.data = f.vec - a.mat*gfu.vec
    gfu.vec.data += inv_stokes * res

    # calculate drag force
    res.data = a.mat*gfu.vec
    F_drag_j = InnerProduct(res,gfu_help.vec) * rho

# density free pressure --> pressure
gfu.components[1].vec.data.FV().NumPy()[:] = gfu.components[1].vec.data.FV().NumPy() * rho


# pressure solution field
pressure = gfu.components[1]
# Define a coefficient function for val

# Draw velocity
#Draw(Norm(gfu.components[0]), mesh, "velocity",sd=3)   
#
Draw((gfu.components[0]), mesh, "velocity",sd=3)          # Draw velocity
# Draw the pressure
#Draw(pressure,mesh)                                         # Draw pressure

Hi Raphael_S,

could you please provide a minimal  (non)-working example? I ran the code below and this set the pressure correctly (although I'm not sure the equations make much sense in this form).

Best wishes, Henry

from netgen.geom2d import unit_square
from ngsolve import *
from ngsolve.webgui import Draw

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))
V = VectorH1(mesh, order=3, dirichlet='bottom|top')  # velocity
Q = H1(mesh,order=2, dirichlet='left|right') # pressure
X = V * Q 
(u, p), (v, q) = X.TnT()

gfu = GridFunction(X)
cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right" else 0 for bc in mesh.GetBoundaries()])
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))

f = LinearForm(X)
f.Assemble()

stokes = 0.1*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
a = BilinearForm(X, symmetric=True)
a += stokes*dx
a.Assemble()

inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
res = f.vec.CreateVector()
res.data = f.vec - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res

print(gfu.components[1](mesh(1,0.5)))
print(gfu.components[1](mesh(0,0.5)))
 

Hy Henry,  Many thanks for your support. As I see, I did not add the whole code. Sorry for this. Many thanks that you support me. Best wishes, Raphael

from netgen.occ import *
#from netgen.webgui import Draw as DrawGeo
from ngsolve import *
from netgen import*

# Velocity
U = 5*10**3 # velcoity in mm/s
rho = 857.22*10**-3/10**9   # density of the oil // t/mm^3
eta = 0.022474*10**-6       # dynamic viscosity of the oil // N/mm^2*s

box = Box((0,0,0), (3,1,1))
geo = box

geo.faces.Min(X).name='left'
geo.faces.Max(X).name='right'

geo.faces.Min(Y).name='front'
geo.faces.Max(Y).name='back'

geo.faces.Min(Z).name='bot'
geo.faces.Max(Z).name='top'

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.1)).Curve(3)
#DrawGeo(geo)
Draw(mesh)

# Define a finite element space
V = VectorH1(mesh,order=3,dirichlet='front|back|bot|top') # velocity
# V = Periodic(VectorH1(mesh,order=3, dirichlet='bot|top'))  # velocity
Q = H1(mesh,order=2, dirichlet='left|right') # pressure

X = V * Q

# obtain test and trial functions
u, p = X.TrialFunction()
v, q = X.TestFunction()
# space for solution data
gfu = GridFunction(X)
# boundary condition
#uin = CoefficientFunction((U,0,0))
#gfu.components[0].Set(uin, definedon=mesh.Boundaries("bot"))
cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right"  else 0 for bc in mesh.GetBoundaries()])
    # Set BC
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))

nu = (eta/rho)  # kinematic viscosity

SetNumThreads(24) #Anzahl der Kerne zum Parallelisieren
with TaskManager():
    # right hand side of the equation
    f = LinearForm(X)
    f.Assemble()

    # left hand side of the equation
    stokes = nu*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q
    a = BilinearForm(X, symmetric=True)
    a += stokes*dx
    a.Assemble()

    # calculate solution data
    inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
    res = f.vec.CreateVector()
    res.data = f.vec - a.mat*gfu.vec
    gfu.vec.data += inv_stokes * res

    # calculate drag force


# density free pressure --> pressure
gfu.components[1].vec.data.FV().NumPy()[:] = gfu.components[1].vec.data.FV().NumPy() * rho

# pressure solution field
pressure = gfu.components[1]
# Define a coefficient function for val

# Draw velocity
#Draw(Norm(gfu.components[0]), mesh, "velocity",sd=3)   
#
#Draw((gfu.components[0]), mesh, "velocity",sd=3)          # Draw velocity
# Draw the pressure
Draw(pressure)                                         # Draw pressure

pInitLeft = gfu.components[1](mesh(0,0.5,0.5))
pInitRight = gfu.components[1](mesh(3,0.5,0.5))
## data visualisation / export

velocity = gfu.components[0]
nameOfFile = 'test'
    #nameOfFile = 'Initial_Conditions'
    # VTKOutput object
vtk = VTKOutput(ma=mesh,
                    coefs=[pressure,velocity],
                    names = ["pressure","velocity"],
                    filename=nameOfFile,
                    subdivision=0,
                    legacy = True)
    # Exporting the results:
vtk.Do()

 

Dear Raphael,

if you post code that isn't working as expected, please make a minimal working example, i.e., remove all code that isn't relevant to the issue and that runs on a laptop. Something like this:

from netgen.occ import *
from ngsolve import *

U = 5*10**3 # velcoity in mm/s
rho = 857.22*10**-3/10**9 # density of the oil // t/mm^3
eta = 0.022474*10**-6 # dynamic viscosity of the oil // N/mm^2*s

geo = Box((0,0,0), (3,1,1))

geo.faces.Min(X).name='left'
geo.faces.Max(X).name='right'
geo.faces.Min(Y).name='front'
geo.faces.Max(Y).name='back'
geo.faces.Min(Z).name='bot'
geo.faces.Max(Z).name='top'

mesh = Mesh(OCCGeometry(geo).GenerateMesh(maxh=0.5))

V = VectorH1(mesh,order=2,dirichlet='front|back|bot|top') # velocity
Q = H1(mesh,order=1, dirichlet='left|right') # pressure
X = V * Q
u, p = X.TrialFunction()
v, q = X.TestFunction()

gfu = GridFunction(X)

cf = CoefficientFunction([(2) if bc=="left" else (0.2) if bc=="right" else 0 for bc in mesh.GetBoundaries()])
gfu.components[1].Set(cf, definedon=mesh.Boundaries("left|right"))

a = BilinearForm(X, symmetric=True)
a += ((eta/rho)*InnerProduct(grad(u), grad(v))+div(u)*q+div(v)*p - 1e-10*p*q) * dx
a.Assemble()

inv_stokes = a.mat.Inverse(X.FreeDofs(), inverse="sparsecholesky")
res = gfu.vec.CreateVector()
res.data = - a.mat*gfu.vec
gfu.vec.data += inv_stokes * res
print(gfu.components[1](mesh(0, 0.5, 0.5)))
print(gfu.components[1](mesh(3, 0.5, 0.5)))

which again results in the expected output for me.

Best wishes,
Henry