Hello,

I’m trying to implement Stokes equation using equal order elements. To ensure the stability of the solution, I’m using the PSPG with the following stabilization parameter:

\tau = \left[\left(\frac{2\left \| \mathbf{u} \right \|}{h}\right)^2 +\left(\frac{4\nu}{h^2}\right)^2\right]^{-\frac{1}{2}}

I’m having some problems implementing the stabilization. When I define the stabilization parameter with the same velocity used to create the variational form (`u, p = split(w)`

), the code fails at the first iteration returning `nan`

. However, if I use `w.split()[0]`

as the velocity to define the parameter, the code runs normally. I would like to know why the solver fails using the same function of the variational form and what is the correct approach to define the stabilization parameter.

```
from fenics import *
# Parameters
mu = 1e-4
rho = 1.0
nu = mu/rho
U = 1.0
H = 1.0
L = 10.0
# Mesh
mesh = RectangleMesh.create([Point(0.0, -H/2), Point(L, H/2)], [200, 20], CellType.Type.quadrilateral)
# Define Function Spaces
V = VectorElement("Lagrange", mesh.ufl_cell(), 1)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([V, Q]))
# Boundary conditions
Velocity = Constant((U, 0.0))
NoSlip = Constant((0.0, 0.0))
def Inlet(x, on_boundary):
return on_boundary and near(x[0], 0.0)
def Walls(x, on_boundary):
return on_boundary and (near(x[1], 0.5) or near(x[1], -0.5))
def Outlet(x, on_boundary):
return on_boundary and near(x[0], 10.0)
bc1 = DirichletBC(W.sub(0), Velocity, Inlet)
bc2 = DirichletBC(W.sub(0), NoSlip, Walls)
bc3 = DirichletBC(W.sub(0).sub(1), Constant(0.0), Outlet) # Tangential velocity
bcs = [bc1, bc2, bc3]
# Define test and trial functions
w_trial = TrialFunction(W)
w_test = TestFunction(W)
v, q = split(w_test)
# Define function
w = Function(W)
u, p = split(w)
# Strain rate tensor
def D(u):
return sym(grad(u))
# Identity tensor
I = Identity(u.geometric_dimension())
# Total stress
T = -p*I + 2*mu*D(u)
# Variational form
F = inner(grad(v), T)*dx + q*div(u)*dx
# Residual
res = -div(T)
# Element size
h = CellDiameter(mesh)
# Stabilization (PSPG)
#vel = u # Return nan at the first iteration!
vel = w.split()[0] #or w.sub(0) Works with both options!
norm_u = sqrt(dot(vel, vel))# + DOLFIN_EPS
tau = 1/sqrt((2.0*norm_u/h)**2 + (4.0*nu/(h**2))**2)
F += tau*inner(grad(q), res)*dx
J = derivative(F, w, w_trial)
# Solver
solve(F == 0, w, bcs, J=J)
u, p = w.split()
# Files
ufile = File("velocity.pvd")
ufile << u
pfile = File("pressure.pvd")
pfile << p
```