# Error in Solving Steady Navier Stokes Equation

I’m trying to convert my legacy FEniCS codes to FEniCSx to solve steady incompressible NSE with the Newton method for a Lid-Driven Cavity. I’m not able to set up the nonlinear solver properly.

The FEniCS code is:

``````# Define function spaces (Taylor-Hood)
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

# Make a mixed space
TH = V * Q
W = FunctionSpace(mesh, TH)

# Define unknown and test function(s)
(v, q) = TestFunctions(W)
w = Function(W)
(u, p) = (as_vector((w[0], w[1])), w[2])

# Tentative velocity step
- div(v)*p*dx \
- q*div(u)*dx

dw = TrialFunction(W)
dF = derivative(F, w, dw)

nsproblem = NonlinearVariationalProblem(F, w, bc, dF)
solver = NonlinearVariationalSolver(nsproblem)
solver.solve()

# Extract solutions:
(u, p) = w.split()
``````

The FEniCSx code is:

``````msh = mesh.create_unit_square(MPI.COMM_WORLD, 40, 40, mesh.CellType.triangle)

nu = 0.001

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
return np.isclose(x[1], 1.0)

# Lid velocity
def lid_velocity_expression(x):
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = fem.FunctionSpace(msh, P2), fem.FunctionSpace(msh, P1)

# Create the function space
TH = P2 * P1
W = fem.FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.set(0.0)
dofs = fem.locate_dofs_geometrical((W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
(u, p) = ufl.TrialFunctions(W)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step

w = fem.Function(W)
dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = fem.petsc.NonlinearProblem(F, w, bcs=[bcs], J=dF)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Compute the solution
ksp.solve(w.vector)

# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()
``````

This produces an error:

``````Traceback (most recent call last):
problem = fem.petsc.NonlinearProblem(F, w, bcs=[bcs], J=dF)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 683, in __init__
self._L = _create_form(F, form_compiler_options=form_compiler_options,
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 176, in form
return _create_form(form)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 171, in _create_form
return _form(form)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/forms.py", line 145, in _form
ufcx_form, module, code = jit.ffcx_jit(mesh.comm, form,
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/jit.py", line 204, in ffcx_jit
r = ffcx.codegeneration.jit.compile_forms([ufl_object], options=p_ffcx, **p_jit)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 186, in compile_forms
impl = _compile_objects(decl, forms, form_names, module_name, p, cache_dir,
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ffcx/codegeneration/jit.py", line 250, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ffcx/compiler.py", line 97, in compile_ufl_objects
analysis = analyze_ufl_objects(ufl_objects, options)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ffcx/analysis.py", line 86, in analyze_ufl_objects
form_data = tuple(_analyze_form(form, options) for form in forms)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ffcx/analysis.py", line 86, in <genexpr>
form_data = tuple(_analyze_form(form, options) for form in forms)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ffcx/analysis.py", line 164, in _analyze_form
form_data = ufl.algorithms.compute_form_data(
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/algorithms/compute_form_data.py", line 407, in compute_form_data
check_form_arity(preprocessed_form, self.original_form.arguments(), complex_mode)  # Currently testing how fast this is
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 177, in check_form_arity
check_integrand_arity(itg.integrand(), arguments, complex_mode)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 159, in check_integrand_arity
arg_tuples = map_expr_dag(rules, expr, compress=False)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 36, in map_expr_dag
result, = map_expr_dags(function, [expression], compress=compress,
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/corealg/map_dag.py", line 99, in map_expr_dags
r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 63, in product
raise ArityMismatch("Multiplying expressions with overlapping form argument number {0}, argument is {1}.".format(x[0].number(), _afmt(x)))
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 19, in _afmt
return tuple("conj({0})".format(arg) if conj else str(arg)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/ufl/algorithms/check_arities.py", line 20, in <genexpr>
for arg, conj in atuple)
ValueError: too many values to unpack (expected 2)
``````

I think there are mistakes in the nonlinear solver setup and in the trial functions u & p. Any help will be much appreciated.

These should not be trial-functions. They should be a function,
i.e.

``````w = Function(W)
u, p = ufl.split(w)
``````
1 Like

I changed u & p to functions just as you instructed. Now I’m getting an error in the compute statement:

``````# Define variational problem
w = fem.Function(W)
u, p = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = fem.petsc.NonlinearProblem(F, w, bcs=[bcs], J=dF)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Compute the solution
ksp.solve(w)
``````

This gives an error:

``````    ksp.solve(w)
File "PETSc/KSP.pyx", line 398, in petsc4py.PETSc.KSP.solve
TypeError: solve() takes exactly 2 positional arguments (1 given)
``````

In the tutorial solve() had one argument. Please let me know how to correct this.

You shouldn’t call the ksp solver directly, you should call the Newton solver, ie, solver.solve(w)

I did precisely that earlier but that didn’t work either. When I ran this:

``````# Compute the solution
n = solver.solve(w)
``````

It gave out an error:

``````    n = solver.solve(w)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/nls/petsc.py", line 41, in solve
n, converged = super().solve(u.vector)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 730, in F
apply_lifting(b, [self._a], bcs=[self.bcs], x0=[x], scale=-1.0)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/petsc.py", line 494, in apply_lifting
assemble.apply_lifting(b_local.array_w, a, bcs, x0_r, scale, constants, coeffs)
File "/Users/varunkumar/opt/anaconda3/envs/fenicsx/lib/python3.10/site-packages/dolfinx/fem/assemble.py", line 307, in apply_lifting
_cpp.fem.apply_lifting(b, a, constants, coeffs, bcs, x0, scale)
TypeError: apply_lifting(): incompatible function arguments. The following argument types are supported:
1. (b: numpy.ndarray[numpy.float32], a: List[dolfinx::fem::Form<float>], constants: List[numpy.ndarray[numpy.float32]], coeffs: List[Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float32]]], bcs1: List[List[dolfinx::fem::DirichletBC<float>]], x0: List[numpy.ndarray[numpy.float32]], scale: float) -> None
2. (b: numpy.ndarray[numpy.float64], a: List[dolfinx::fem::Form<double>], constants: List[numpy.ndarray[numpy.float64]], coeffs: List[Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.float64]]], bcs1: List[List[dolfinx::fem::DirichletBC<double>]], x0: List[numpy.ndarray[numpy.float64]], scale: float) -> None
3. (b: numpy.ndarray[numpy.complex64], a: List[dolfinx::fem::Form<std::__1::complex<float> >], constants: List[numpy.ndarray[numpy.complex64]], coeffs: List[Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex64]]], bcs1: List[List[dolfinx::fem::DirichletBC<std::__1::complex<float> >]], x0: List[numpy.ndarray[numpy.complex64]], scale: float) -> None
4. (b: numpy.ndarray[numpy.complex128], a: List[dolfinx::fem::Form<std::__1::complex<double> >], constants: List[numpy.ndarray[numpy.complex128]], coeffs: List[Dict[Tuple[dolfinx::fem::IntegralType, int], numpy.ndarray[numpy.complex128]]], bcs1: List[List[dolfinx::fem::DirichletBC<std::__1::complex<double> >]], x0: List[numpy.ndarray[numpy.complex128]], scale: float) -> None

Invoked with: array([0., 0., 0., ..., 0., 0., 0.]), [<dolfinx.fem.forms.Form object at 0x1343608b0>], [array([], dtype=float64)], [{(<IntegralType.cell: 0>, -1): array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])}], [[[<dolfinx.fem.bcs.DirichletBC object at 0x100b6afc0>, <dolfinx.fem.bcs.DirichletBC object at 0x1342cf240>, <dolfinx.fem.bcs.DirichletBC object at 0x1342cf420>]]], [array([0., 0., 0., ..., 0., 0., 0.])], -1.0

Did you forget to `#include <pybind11/stl.h>`? Or <pybind11/complex.h>,
<pybind11/functional.h>, <pybind11/chrono.h>, etc. Some automatic
conversions are optional and require extra headers to be included
``````

I also tried (u, p) = (ufl.as_vector((w[0], w[1])), w[2]) but that didn’t work either. I tried this from the tutorial:

``````log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(uh)
assert(converged)
print(f"Number of interations: {n:d}")
``````

Replace `bcs=[bcs]` with `bcs=bcs`.

Full working code:

``````from dolfinx import mesh, fem, nls
from mpi4py import MPI
import numpy as np
import ufl
from ufl import inner, div, grad, dx
from petsc4py import PETSc

msh = mesh.create_unit_square(MPI.COMM_WORLD, 40, 40, mesh.CellType.triangle)

nu = 0.001

# Function to mark x = 0, x = 1 and y = 0

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

# Function to mark the lid (y = 1)

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

# Lid velocity#

def lid_velocity_expression(x):
return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))

P2 = ufl.VectorElement("Lagrange", msh.ufl_cell(), 2)
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(), 1)
V, Q = fem.FunctionSpace(msh, P2), fem.FunctionSpace(msh, P1)

# Create the function space
TH = P2 * P1
W = fem.FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.set(0.0)
dofs = fem.locate_dofs_geometrical(
(W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step
dx - div(v)*p*dx - q*div(u)*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = fem.petsc.NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()

# Compute the solution
solver.solve(w)

# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()

``````