Hello everyone!
I have struggled with nonlinear stokes problem for 2 weeks.
I use
mixed_element to define the FEspace, and I use NonlinearProblem solver to solve this problem. Really need your help!! 
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 
. 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!