Hello everyone,

I was trying to solve a steady pressure driven Stokes channel flow by using the Chorin scheme with neglecting the innertial term. But surprisingly I found the velocity field converge to a false time step length dependent value. I list my MWE below, with which I found when dt=0.1 , 0.2, 0.4, u_max=1.1, 1.2, 1.4 respectively. Did i miss some fundamental modification that required for solving a steady problem with false time approaching? Thanks!

```
from dolfin import *
# Print log messages only from the root process in parallel
parameters["std_out_all_processes"] = False;
# Load mesh from file
N = 32
mesh = UnitSquareMesh(N, N)
# Inflow boundary
class InflowBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[0] < DOLFIN_EPS
# Outflow boundary
class OutflowBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1 - DOLFIN_EPS
# No-slip boundary
class NoslipBoundary(SubDomain):
def inside(self, x, on_boundary):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
# Define function spaces (P2-P1), but not mixed space
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
# Define trial and test functions
u = TrialFunction(V)
p = TrialFunction(Q)
v = TestFunction(V)
q = TestFunction(Q)
# Create functions
u0 = Function(V)
u1 = Function(V)
p1 = Function(Q)
# ------------------------------------------------------------------------
# Define boundary conditions
# ------------------------------------------------------------------------
# Define pressure boundary condition
def pressure_bc(Q):
element = FiniteElement("CG", triangle, 1)
return Expression("1 - x[0]", element=element)
# no-slip
bv = DirichletBC(V, Constant((0.0, 0.0)), NoslipBoundary())
# Create boundary conditions for pressure
bp0 = DirichletBC(Q, pressure_bc(Q), InflowBoundary())
bp1 = DirichletBC(Q, pressure_bc(Q), OutflowBoundary())
bcu = [bv]
bcp = [bp0, bp1]
# Set parameter values
dt = 0.4
T = 30
#nu = 0.01
nu = Constant(1.0/8.0)
# Define coefficients
k = Constant(dt)
f = Constant((0, 0))
# Tentative velocity step
#F1 = (1/k)*inner(u - u0, v)*dx + \
# nu*inner(grad(u), grad(v))*dx - inner(f, v)*dx
#a1 = lhs(F1)
#L1 = rhs(F1)
a1 = (1/k)*inner(u, v)*dx + nu*inner(grad(u), grad(v))*dx
L1 = (1/k)*inner(u0, v)*dx + inner(f, v)*dx
# Pressure update
a2 = inner(grad(p), grad(q))*dx
L2 = -(1/k)*div(u1)*q*dx
# Velocity update
a3 = inner(u, v)*dx
L3 = inner(u1, v)*dx - k*inner(grad(p1), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
# Create files for storing solution
ufile = File("results/velocity.pvd")
pfile = File("results/pressure.pvd")
# Time-stepping
#t = dt
t = 0.0
while t < T + DOLFIN_EPS:
# Update pressure boundary condition
#p_in.t = t
# Compute tentative velocity step
print("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(A1, b1) for bc in bcu]
solve(A1, u1.vector(), b1, "gmres", "default")
# Pressure correction
print("Computing pressure correction")
b2 = assemble(L2)
[bc.apply(A2, b2) for bc in bcp]
solve(A2, p1.vector(), b2, "cg", prec)
# Velocity correction
print("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(A3, b3) for bc in bcu]
solve(A3, u1.vector(), b3, "gmres", "default")
# calculate the norm of difference between two time step
u_delta = u1.vector() - u0.vector()
print("error nomr =", norm(u_delta, 'linf'))
# Move to next time step
u0.assign(u1)
t += dt
print("t =", t)
print("u_max", u1.vector().max())
```