Confused about how to solve nonlinear stokes equation


Hello everyone!
I have struggled with nonlinear stokes problem for 2 weeks. :innocent:
I use mixed_element to define the FEspace, and I use NonlinearProblem solver to solve this problem. Really need your help!! :heart:
The error is:

yls@hty:~$  /usr/bin/env /usr/bin/python3 /home/yls/.vscode-server/extensions/ms-python.debugpy-2025.6.0-linux-x64/bundled/libs/debugpy/adapter/../../debugpy/launcher 45545 -- /home/yls/ns-dolfinx/nonlinear-stokes.py 
At column        0, pivotL() encounters zero diagonal at line 723 in file ./SRC/symbfact.c

I tried to fix the pressure at the point: (0.0, -0.25), and use dirichletbc.
But I think it not works as I think :neutral_face:. Perhaps there are other bugs I haven’t noticed…

def marker(x): return np.isclose(x[0], 0.0) & np.isclose(x[1], -0.25)
p_dofs = fem.locate_dofs_geometrical((ME.sub(1), S), marker) 
bc_pressure = dirichletbc(p_bc, p_dofs, ME.sub(1))

This is the full code:

from mpi4py import MPI
import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, la, log, mesh
from dolfinx.fem import (
    Constant,
    Function,
    dirichletbc,
    form,
    functionspace,
)
from petsc4py import PETSc
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block, NonlinearProblem, LinearProblem
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, inner, nabla_grad, dot, TestFunction, TrialFunction, TestFunctions, TrialFunctions, SpatialCoordinate
import ufl.functionspace
from dolfinx.nls.petsc import NewtonSolver


# Solution
def u_x(x):
    return x[0]**2 * x[1]**2 + np.exp(-x[1])


def u_y(x):
    return -2/3 * x[0] * x[1]**3 + 2 - np.pi * np.sin(np.pi * x[0])


def pressure(x):
    return -(2 - np.pi * np.sin(np.pi * x[0])) * np.cos(2 * np.pi * x[1])


def true_velocity(x):
    return [u_x(x), u_y(x)]


def f_x(x_list):
    x, y = x_list  
    term1 = -2 * nu * x ** 2 - 2 * nu * y ** 2 - nu * ufl.exp(-y)
    term2 = ufl.pi ** 2 * ufl.cos(ufl.pi * x) * ufl.cos(2 * ufl.pi * y)
    term3 = 2 * x * y ** 2 * (x ** 2 * y ** 2 + ufl.exp(-y))
    term4 = (-2 * x * y ** 3 / 3 + 2 - ufl.pi * ufl.sin(ufl.pi * x)) * \
        (2 * x ** 2 * y - ufl.exp(-y))
    return term1 + term2 + term3 + term4


def f_y(x_list):
    x, y = x_list  
    term1 = 4 * nu * x * y - nu * ufl.pi ** 3 * ufl.sin(ufl.pi * x)
    term2 = 2 * ufl.pi * (2 - ufl.pi * ufl.sin(ufl.pi * x)) * ufl.sin(2 * ufl.pi * y)
    term3 = (x ** 2 * y ** 2 + ufl.exp(-y)) * \
        (-2 * y ** 3 / 3 - ufl.pi ** 2 * ufl.cos(ufl.pi * x))
    term4 = (-2 * x * y ** 3 / 3 + 2 - ufl.pi *
             ufl.sin(ufl.pi * x)) * (-2 * x * y ** 2)
    return term1 + term2 + term3 + term4

def f(x):
    return [f_x(x), f_y(x)]


def D(u):
    return 0.5 * (ufl.nabla_grad(u) + ufl.nabla_grad(u).T)

# Params
n_list = [
    [4, 1],
    [8, 2],
    [16, 4],
    [32, 8],
    [64, 16],
]
idx = 3
nu = 1

msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, -0.25]), np.array(
        [1, 0])], n=n_list[idx], cell_type=CellType.triangle
)

# Function Space
P2 = element(
    "Lagrange", msh.basix_cell(), degree=2, shape=(msh.geometry.dim,), dtype=default_real_type
)
P1 = element("Lagrange", msh.basix_cell(), degree=1, dtype=default_real_type)
mel = mixed_element([P2, P1])  
ME = functionspace(msh, mel)  

x = SpatialCoordinate(msh)  


Q, _ = ME.sub(0).collapse()
S, _ = ME.sub(1).collapse()
u_bc = Function(Q)  
p_bc = Function(S)  
u_bc.interpolate(true_velocity)
p_bc.interpolate(pressure)

sol = Function(ME)  
sol.x.array[:] = 0.0  

u, p = sol.split()  # trial: u, p; test: v, q
v, q = TestFunctions(ME)  # Test Functions

# Weak Formulation
F = 2 * nu * inner(D(u), D(v)) * dx  
F += inner(div(u), q) * dx
F -= inner(div(v), p) * dx
F += inner(dot(u, nabla_grad(u)), v) * dx  
f_expr = ufl.as_vector(f(x)) 
F -= inner(f_expr, v) * dx  
 

# BCs
fdim = msh.topology.dim - 1
facets = mesh.locate_entities_boundary(
    msh, fdim, lambda x: np.isclose(x[0], 1.0) | np.isclose(
        x[0], 0.0) | np.isclose(x[1], -0.25) | np.isclose(x[1], 0.0)
)

v_dofs = fem.locate_dofs_topological((ME.sub(0), Q), fdim, facets)
bc_vel = dirichletbc(u_bc, v_dofs, ME.sub(0))  

# Not Sure: How to set pressure?
def marker(x): return np.isclose(x[0], 0.0) & np.isclose(x[1], -0.25)
p_dofs = fem.locate_dofs_geometrical((ME.sub(1), S), marker) 
bc_pressure = dirichletbc(p_bc, p_dofs, ME.sub(1))


# Create nonlinear problem and Newton solver
problem = NonlinearProblem(F, sol, bcs=[bc_vel, bc_pressure])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options()  # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
sys = PETSc.Sys()  # type: ignore
# For factorisation prefer superlu_dist, then MUMPS, then default
if sys.hasExternalPackage("superlu_dist"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "superlu_dist"
elif sys.hasExternalPackage("mumps"):
    opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()


n, converged = solver.solve(sol)

Looking forward to your help! :heart:

Use ufl.split rather than sol.split

Thank you so much. It works!
But newton solver doesn’t converge. I need to check my script again:)