Problems with defining pressure BC for Stokes equation

NGSolve
1

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

2

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.1InnerProduct(grad(u), grad(v))+div(u)q+div(v)p - 1e-10pq
a = BilinearForm(X, symmetric=True)
a += stokesdx
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.components1)
print(gfu.components1)

3

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 = 510**3 # velcoity in mm/s
rho = 857.2210**-3/109 # 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.components1
pInitRight = gfu.components1

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

4

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 = 510**3 # velcoity in mm/s
rho = 857.2210**-3/109 # 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-10pq) * 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