Solving Stokes equation with SUPG - Error in Newton Solver

General
,
1

Hello Fenics Community,

I am solving stokes equation using SUPG stabilization. I am setting this problem for lid-driven cavity flow problem. I am using mixed finite element spaces. I am facing issues in the Nonlinear solver. The solver is not converging. I have tried changing solver options, but failed. Can someone help me with it? Am I missing anything in the formulation or solver setup?

I am using FeNiCSx 0.9.0. The MWE is as follows:


import numpy as np
from dolfinx import fem, mesh, nls, io
from mpi4py import MPI
import ufl
from basix.ufl import element, mixed_element
from petsc4py import PETSc

from dolfinx.nls.petsc import NewtonSolver

# Create mesh (Unit square domain)
domain = mesh.create_unit_square(MPI.COMM_WORLD, 50, 50)

def noslip_boundary(x):
    return np.isclose(x[1], 0.0) | np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0)

def lid_boundary(x):
    return np.isclose(x[1], 1.0)

noslip_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, noslip_boundary)
lid_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, lid_boundary)

# Define equal-order elements (P1-P1 discretization)
P1v = element("Lagrange", domain.basix_cell(), 2,  shape=(domain.geometry.dim,) )  # Velocity (P1)
P1p = element("Lagrange", domain.basix_cell(), 1)  # Pressure (P1)
TH = mixed_element([P1v, P1p])

VQ= fem.functionspace(domain, TH)

(v, q) =ufl.TestFunctions(VQ)

# Define solution function
up = fem.Function(VQ)  # Nonlinear unknown (velocity-pressure)
u, p = ufl.split(up)

# Fluid properties
nu_val = 0.01  # Kinematic viscosity
f = fem.Constant(domain, PETSc.ScalarType((0.0, 0.0)))  # Body force

# SUPG stabilization parameters
h = ufl.CellDiameter(domain)
vel_norm = ufl.sqrt(ufl.inner(u, u) + 1e-6)
Pe = vel_norm * h / (2 * nu_val)
Pe = ufl.max_value(Pe, 1e-6)  # Avoid zero Pe
tau_SUPG = h / (2 * vel_norm) * (ufl.exp(Pe) - ufl.exp(-Pe)) / (ufl.exp(Pe) + ufl.exp(-Pe)) - 1/Pe

# Weak form of Stokes equations
F = (nu_val * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx
     - p * ufl.div(v) * ufl.dx
     - q * ufl.div(u) * ufl.dx
     - ufl.inner(f, v) * ufl.dx)

# SUPG stabilization term
R_m = -nu_val * ufl.div(ufl.grad(u)) + ufl.grad(p) - f
F_SUPG = tau_SUPG * ufl.inner(R_m, ufl.dot(u, ufl.grad(v))) * ufl.dx
F += F_SUPG

# Get subspace of V
V0 = VQ.sub(0)
V, _ = V0.collapse()

# Get subspace of Q
Q0 = VQ.sub(1)
Q, _ = Q0.collapse()

# Define boundary conditions
zero_velocity = fem.Function(V)
zero_velocity.interpolate(lambda x: np.zeros((2, x.shape[1])))
lid_velocity = fem.Function(V)
lid_velocity.interpolate(lambda x: np.stack((np.ones_like(x[0]), np.zeros_like(x[1]))))

noslip_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, noslip_facets)
lid_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, lid_facets)

bc_noslip = fem.dirichletbc(zero_velocity, noslip_dofs)
bc_lid = fem.dirichletbc(lid_velocity, lid_dofs)
bcs = [bc_noslip, bc_lid]

if len(bcs) == 0:
    raise ValueError("No boundary conditions found. Check domain boundaries.")

# Define nonlinear problem
uh = fem.Function(VQ)

duh = ufl.TrialFunction(VQ)
dF = ufl.derivative(F, uh, duh)

problem = fem.petsc.NonlinearProblem(F, uh, bcs=bcs, J=dF)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
#nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-3  # Looser relative tolerance
solver.atol = 1e-3  # Looser absolute tolerance
solver.max_it = 200  # Allow more iterations
solver.report = True  # Print iteration details

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"  # Use a direct solver (no Krylov iterations)
opts[f"{option_prefix}pc_type"] = "lu"        # LU factorization (direct solve)
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"  # Use MUMPS as the direct solver
ksp.setFromOptions()

# Compute the solution
solver.solve(uh)

# Extract velocity and pressure solutions
uh_velocity = fem.Function(V)
uh_pressure = fem.Function(Q)
uh_velocity.interpolate(uh.sub(0))
uh_pressure.interpolate(uh.sub(1))

# Save results
with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(uh_velocity)

with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as file:
    file.write_function(uh_pressure)
2

This code is not properly formatted. Please use

```python
#add code here
```

and ensure that people can copy-paste, then run the code.

3

Thanks for pointing out. Fixed it!

4

Hi,

Sorry I can’t run your code right now, but it might have to deal with your boundary dofs. Instead of:

zero_velocity = fem.Function(V)
zero_velocity.interpolate(lambda x: np.zeros((2, x.shape[1])))
noslip_dofs = fem.locate_dofs_topological(V, domain.topology.dim - 1, noslip_facets)

Try:

zero_velocity = fem.Function(VQ)
zero_velocity.sub(0).interpolate(lambda x: np.zeros((2, x.shape[1])))
noslip_facets = mesh.locate_entities_boundary(msh, msh.topology.dim - 1, periodic_marker)
noslip_dofs = fem.locate_dofs_topological(VQ, msh.topology.dim - 1, noslip_facets)
bc_noslip = fem.dirichletbc(zero_velocity, noslip_dofs )

Same for lid. If you want to have the dof map between VQ and V, use:

V, V_map = VQ.sub(0).collapse()
Q, Q_map = VQ.sub(1).collapse()

Let me know !

5

Hi TheauC,

Thanks for the response. The problem still persist. The error is same as previous. I am copy-pasting it here.

File "/Users/tarungangwar/Documents/0_work/Research/Projects/try/stokes.py", line 118, in <module>
    solver.solve(uh)
    ~~~~~~~~~~~~^^^^
  File "/Users/tarungangwar/anaconda3/envs/fenicsx-env/lib/python3.13/site-packages/dolfinx/nls/petsc.py", line 51, in solve
    n, converged = super().solve(u.x.petsc_vec)
                   ~~~~~~~~~~~~~^^^^^^^^^^^^^^^
RuntimeError: Newton solver did not converge because maximum number of iterations reached
6

Your dF is wrong.
It should be differentiated with respect to up, which is the variable you use within the residual formulation, not a new variable uh that is not used anywhere. ufl.derivative(F, uh, duh) is naturally zero, as F is not a function of uh.

7

As a side note, if I am not mistaken, you’ll want
F_SUPG = tau_SUPG * ufl.inner(R_m, ufl.dot(u, ufl.nabla_grad(v))) * ufl.dx
so, nabla_grad instead of grad. Acting on vectors, the resulting tensors are transposes of each other.

Also, (and I imagine this is besides the point, as you probably intend to use this code for Navier-Stokes,) SUPG for a Stokes formulation doesn’t actually stabilize. It needs the advection term of the PDE so that \int \tau R_m \cdot ((\mathbf{u}\cdot\nabla )\mathbf{v}) contributes a positive-definite term (namely \tau ((\tilde{\mathbf{u}}\cdot\nabla) \mathbf{u}) \cdot ( (\tilde{\mathbf{u}}\cdot\nabla ) \mathbf{v}) .

1 Like
8

Hi Stein,

Thanks. I tried the above suggestions; but, the solver did not converge.

I did not know that you are still active on this forum. Yes, I want to eventually implement NS solver using SUPG. I am using Stokes to setup and test FeNiCSx mixed spaces and Nonlinear solver. Perhaps, I should first implement stokes with Taylor-Hood using linear solver and then proceed for more complex problems.

9

I could well imagine that is because of this:

(It actually destabilizes)