Error in Solving Steady Navier Stokes Equation

Errors
1

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
F =   inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v))*dx \
    - 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
F = inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v))*dx - div(v)*p*dx - q*div(u)*dx

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):
  File "/Users/varunkumar/Downloads/FEniCSx v0.6/SteadyNSE-Lid.py", line 67, in <module>
    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.

2

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

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

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
F = inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v))*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
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.

4

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

5

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
when compiling your pybind11 module.

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}")

but this too gave the same error. Please help me.

6

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
F = inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v)) * \
    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()