Implementing the Stokes equation with PSPG stabilization

variational formulation
1

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
2

@fforteneto Can I ask the question: what is W.sub(0).sub(1)? why is not use the code: W.sub(1)

3

Hi, I’m using the command W.sub(0).sub(1) to impose the boundary condition in the y-componente of the vector field, since I want to avoid some exit effects. Using W.sub(1), the boundary condition will be imposed on the pressure, however, the pressure is already prescribed since I’m considering the traction equal to zero at the outlet.

5

I still can’t make the solver converges, but I come up with the following hypotheses:

  • I think that to use u instead of w.split()[0], is necessary to change the solver configuration. I have tried to use GMRES with different preconditioners, but I still get the same error. Reading some topics in the group I noticed that a similar problem occurred with the SUPG stabilization link.

  • Is there a command in FEniCS similar to COMSOL’s nojac()?. I’m asking that because I think that use w.split()[0] might have a similar effect as nojac(), since w.split()[0] it is a deepcopy of the w and is probably not contributing to the Jacobian which is calculated with respect to w.

I’m insisting on this topic because I’m going to work with a more complex stabilization parameter that involves other variable and their gradients. That way, I would be very grateful if someone could help me with this topic.