Pressure induced flow

Hy everybody, I am quiet new to NGSolve and so I need some help. I would like to solve the Stokes equation for a pressure induced flow (cf. flow inside a pipe).  In the tutorials there is a velocity induced flow presented. In my case, I would like to apply a pressure difference on the left and right face of my geometry, which are called inlet and outlet. In addition, I added a Neumann BC on both sides. The major problem now is, that the fluid is not moving in the centre of the geometry (clip through the y axis shows this). Furthermore, if the pressure difference is increased, the max. velocity is not rising.  I am vey thankful for every support.

from ngsolve import *
import netgen.gui
from netgen.occ import *
from netgen.geom2d import SplineGeometry
import numpy as np

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


nu = (eta/rho)  # kinematic viscosity
geo = Box((0,0,0), (3,1,1))

geo.faces.Min(X).name='inlet'
geo.faces.Max(X).name='outlet'
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))

Draw (mesh)

V = VectorH1(mesh, order=3, dirichlet="front|back|bot|top")
Q = H1(mesh, order=2,dirichlet="inlet|outlet")

u,v = V.TnT()
p,q = Q.TnT()

a = BilinearForm(V)
a += nu*InnerProduct(grad(u), grad(v))*dx

b = BilinearForm(trialspace=V, testspace=Q)
b += div(u)*p*dx+div(v)*q*dx


a.Assemble()
b.Assemble()

mp = BilinearForm(Q)
mp += SymbolicBFI(1e-10*p*q)
mp.Assemble()

f = LinearForm(V)
f+= CoefficientFunction((1,0,0))* v * ds(definedon='outlet')
f += CoefficientFunction((1,0,0))* v * ds(definedon='inlet')
f.Assemble()


g = LinearForm(Q)
g.Assemble()

gfu = GridFunction(V, name="u")
gfp = GridFunction(Q, name="p")

cf = CoefficientFunction([(2/rho) if bc=="inlet" else (0.2/rho) if bc=="outlet"  else 0 for bc in mesh.GetBoundaries()])
gfp.Set(cf, definedon=mesh.Boundaries("inlet|outlet"))

K = BlockMatrix( [ [a.mat, b.mat.T], [b.mat, None] ] )
C = BlockMatrix( [ [a.mat.Inverse(V.FreeDofs()), None], [None, mp.mat.Inverse()] ] )

rhs = BlockVector ( [f.vec, g.vec] )
sol = BlockVector( [gfu.vec, gfp.vec] )

solvers.MinRes (mat=K, pre=C, rhs=rhs, sol=sol, initialize=False)
# density free pressure --> pressure
gfp.vec.data.FV().NumPy()[:] = gfp.vec.data.FV().NumPy() * rho

pressure = gfp
velocity = gfu

Draw(gfu)

Dear Raphael,

from a mathematical point of view you can not impose a Dirichlet boundary condition on the pressure, certainly not on the pressure gradient, since $p \in L^2_0$. For Stokes equations we can only set Dirichlet boundary conditions on the velocity. The natural boundary conditions (Neumann) arising from the derivation of the weak formulation, are so called force boundary conditions, namely
$ \int_{\Gamma_N} \left( \mathbf{\tau} - p\mathbf{I} \right) n \cdot \mathbf{v}$.

What you can do, is to think of the pressure gradient as a imposed force acting on the fluid in the whole domain, and not only on the boundary. Therefore it suffices to add the imposed force to the linearform
$\int_{\Omega} \mathbf{f} \cdot \mathbf{v}$.

Best,
Jan

###Code###
from ngsolve import *
from netgen.occ 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
nu = (eta/rho) # kinematic viscosity
imposed_pressure_gradient = CF((-2, 0, 0)) # pressure decreases in x

geo = Box((0, 0, 0), (3, 1, 1))
geo.faces.Min(X).name = 'inlet'
geo.faces.Max(X).name = 'outlet'
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.25))

V = VectorH1(mesh, order=2, dirichlet="front|back|bot|top")
Q = H1(mesh, order=1)

fes = V*Q
(u, p), (v, q) = fes.TnT()

a = BilinearForm(fes)
a += nu*InnerProduct(grad(u), grad(v))*dx
a += -div(u)*p*dx-div(v)*q*dx
a += 1e-8 * p * q * dx
a.Assemble()

f = LinearForm(fes)
f += -imposed_pressure_gradient * v * dx
f.Assemble()

gfu = GridFunction(fes)
velocity = gfu.components[0]
pressure = gfu.components[1]

inv = a.mat.Inverse(freedofs=fes.FreeDofs(), inverse="umfpack")
gfu.vec.data = inv * f.vec

Draw(velocity, mesh, "velocity")
Draw(pressure, mesh, "pressure")