# Incompressible Navier-Stokes with particular scheme & mixed elements FEniCSx

I am having trouble solving incompressible Navier-Stokes with a scheme I will detail below. The mesh generation is following Dokken’s FEniCSx tutorial for the flow past a cylinder (Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial).

Here are the equations governing my problem, first one is the momentum equation, second is the incompressibility constraint:

I am employing a BDF (Backward Differentiation formula) of order 1 or 2 to approximate the time derivative (equation (3)), relaxing the incompressibility constraint (equation (4)) and using a Newton approximation for the convective term (equation (5)):

\mathrm{BDF}_1(\mathbf u)=\frac{\mathbf u^n-\mathbf u^{n-1}} {\Delta t},\qquad\mathrm{BDF}_2(\mathbf u)=\frac{3\mathbf u^n-4\mathbf u^{n-1}+\mathbf u^{n-2}}{2\Delta t},\qquad\textcolor{blue}{(3)}\\ \mathbf\nabla\cdot\mathbf u+\epsilon p=0\qquad\textcolor{blue}{(4)}\\ \mathbf u^n\overline\otimes\mathbf{\nabla u}^n\simeq\mathbf u^n\overline\otimes\mathbf{\nabla u}^{n-1}+\mathbf u^{n-1}\overline\otimes\mathbf{\nabla u}^{n}-\mathbf u^{n-1}\overline\otimes\mathbf{\nabla u}^{n-1}\qquad\textcolor{blue}{(5)}

With that in mind, we can derive the weak formulation (after some tensor calculations) and obtain equation (6):

\int_\Omega\mathrm{BDF}_k(\mathbf u)\overline\otimes\mathbf v\mathrm dx+\int_\Omega\left(\mathbf u^n\overline\otimes\mathbf{\nabla u}^{n-1}+\mathbf u^{n-1}\overline\otimes\mathbf{\nabla u}^{n}-\mathbf u^{n-1}\overline\otimes\mathbf{\nabla u}^{n-1}\right)\overline\otimes\mathbf v\mathrm d x\\ +\frac{1}{Re}\int_\Omega\nabla\mathbf u^n\overline{\overline\otimes}\mathbf{\nabla v}\mathrm d x-\int_\Omega p^n\mathbf \nabla\cdot\mathbf v\mathrm d x-\int_\Omega\mathbf f^n\overline\otimes\mathbf v\mathrm d x\qquad\textcolor{blue}{(6)}\\ -\int_\Omega q\mathbf\nabla\cdot\mathbf u^n\mathrm d x-\epsilon\int_\Omega p^nq\mathrm d x=0

where
\mathbf\nabla,\;\mathbf\nabla\cdot,\;\overline\otimes,\;\overline{\overline\otimes},\;\Omega,\;\mathbf u,\;p,\;\mathbf v,\;q,\;\mathbf f,\;n,\;k,\;\epsilon,Re
are respectively gradient, divergence operators, simple and double dot products, the space we integrate on, the velocity vector (2D), pressure, test functions for velocity and pressure, volume forces, step, order of the BDF (Backward Differential Formula, defining the time derivative of velocity), a relaxation constant (used to relax the incompressibility constraint) and the Reynolds number (equations are not dimensional).

In my script, I use first order BDF for the first time step then the second order BDF, a Reynolds number of 500, a small relaxation constant and no volumic forces such that:
Re=500, k\in\lbrace1,2\rbrace, \epsilon=10^{-6}, \mathbf f=\mathbf0

Let’s now describe the boundary conditions, they are almost the same as in the tutorial, I just modified the inlet profile so that it does not depend on time (equation (7)):

1. null velocity on walls and cylinder
2. null pressure on the outlet
3. inlet velocity profile following
\forall (x,y)\in\lbrace0\rbrace\times[0,H],\;\mathbf u(x,y)=\frac{4U_0y(H-y)}{H^2},\\ \forall (x,y)\in([0,L]\times\lbrace0,H\rbrace)\cup\partial S,\;\mathbf u(x,y)=\mathbf 0\qquad\textcolor{blue}{(7)}\\ \forall (x,y)\in\lbrace L\rbrace\times[0,H],\;p(x,y)=0

where H,\;L,\;U_0,\;\partial S are respectivly the heigth and length of the rectangle, the lengthwise/longitudinal velocity at the middle of the inlet and the circle delimiting the cylinder (the obstacle).

Now I have to explain what I do in the code. The error occurs at first step, so we only consider first order BDF as it is the one used to calculate the solution at step 1. Let’s go:

First I import everything and create the mesh like in the tutorial (I modified the sizes of the rectangle and cylinder, and modified some other stuff, but I am very confident the issue does not lie here)

############################################
#--------------- Imports ------------------#
############################################
import gmsh
import os
import numpy as np
import matplotlib.pyplot as plt
import tqdm.autonotebook
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx.cpp.mesh import to_type, cell_entity_type
from dolfinx.fem import (Constant, Function, FunctionSpace, VectorFunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_topological, set_bc, Expression)
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix, assemble_vector, create_vector, create_matrix, set_bc, LinearProblem, NonlinearProblem)
from dolfinx.geometry import BoundingBoxTree, compute_collisions, compute_colliding_cells
from dolfinx.io import (VTXWriter, distribute_entity_data, gmshio)
from dolfinx.mesh import create_mesh, meshtags_from_entities
from ufl import (FacetNormal, FiniteElement, Identity, Measure, TestFunction, TestFunctions, TrialFunction, VectorElement, MixedElement, as_vector, div, dot, ds, dx, inner, lhs, grad, nabla_grad, rhs, sym, curl, derivative, split)
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from os import path, mkdir
from dolfinx import log

mesh_comm = MPI.COMM_WORLD

############################################
#--------------- Mesh creation ------------#
############################################

model_rank = 0
gmsh.initialize()

# Dimensions
L = 1.0 # Length of the rectangle
H = 0.6 # Height of the rectangle
c_x = c_y = 0.3 # coordinates of the center of the cylinder (also refered a the obstacle)
r = 0.05 # radius of the cylinder (obstacle)
gdim = 2
model_rank = 0

## Creation of the geometry entities (a cylinder within a rectangle)
if mesh_comm.rank == model_rank:
rectangle = gmsh.model.occ.addRectangle(0, 0, 0, L, H, tag=1)
obstacle = gmsh.model.occ.addDisk(c_x, c_y, 0, r, r)

## The next step is to substract the obstacle from the channel, such that we do not mesh the interior of the circle
if mesh_comm.rank == model_rank:
fluid = gmsh.model.occ.cut([(gdim, rectangle)], [(gdim, obstacle)])
gmsh.model.occ.synchronize()

## To get GMSH to mesh the fluid, we add a physical volume marker
fluid_marker = 1
if mesh_comm.rank == model_rank:
volumes = gmsh.model.getEntities(dim=gdim)
assert(len(volumes) == 1)
gmsh.model.setPhysicalName(volumes, fluid_marker, "Fluid")

## To tag the different surfaces of the mesh, we tag the inflow (left hand side) with marker 2, the outflow (right hand side)
## with marker 3 and the fluid walls with 4 and obstacle with 5. We will do this by computing the center of mass for each geometrical
## entity.
inlet_marker, outlet_marker, wall_marker, obstacle_marker = 2, 3, 4, 5
inflow, outflow, walls, obstacle = [], [], [], []
if mesh_comm.rank == model_rank:
boundaries = gmsh.model.getBoundary(volumes, oriented=False)
for boundary in boundaries:
center_of_mass = gmsh.model.occ.getCenterOfMass(boundary, boundary)
if np.allclose(center_of_mass, [0, H/2, 0]):
inflow.append(boundary)
elif np.allclose(center_of_mass, [L, H/2, 0]):
outflow.append(boundary)
elif np.allclose(center_of_mass, [L/2, H, 0]) or np.allclose(center_of_mass, [L/2, 0, 0]):
walls.append(boundary)
else:
obstacle.append(boundary)
gmsh.model.setPhysicalName(1, wall_marker, "Walls")
gmsh.model.setPhysicalName(1, inlet_marker, "Inlet")
gmsh.model.setPhysicalName(1, outlet_marker, "Outlet")
gmsh.model.setPhysicalName(1, obstacle_marker, "Obstacle")

## Here we will have variable mesh sizes to resolve the flow solution in the area of interest; close to the circular obstacle. To
## do this, we use GMSH Fields.
res_factor = 2
res_min = r / res_factor
if mesh_comm.rank == model_rank:
gmsh.model.mesh.field.setNumbers(distance_field, "EdgesList", obstacle)
gmsh.model.mesh.field.setNumber(threshold_field, "IField", distance_field)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMin", res_min)
gmsh.model.mesh.field.setNumber(threshold_field, "LcMax", 0.5 * H)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMin", r)
gmsh.model.mesh.field.setNumber(threshold_field, "DistMax", 2 * H)
gmsh.model.mesh.field.setNumbers(min_field, "FieldsList", [threshold_field])
gmsh.model.mesh.field.setAsBackgroundMesh(min_field)

### GENERATING THE MESH
## We are now ready to generate the mesh. However, we have to decide if our mesh should
## consist of triangles or quadrilaterals. Here we use second order quadrilateral elements.
if mesh_comm.rank == model_rank:
gmsh.option.setNumber("Mesh.Algorithm", 8)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
gmsh.option.setNumber("Mesh.RecombineAll", 1)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
gmsh.model.mesh.generate(gdim)
gmsh.model.mesh.setOrder(2)
gmsh.model.mesh.optimize("Netgen")

## As we have generated the mesh, we now need to load the mesh and corresponding
## facet markers into DOLFINx. See https://jsdokken.com/src/tutorial_gmsh.html for explanations
mesh, c, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
mesh.name = "Cylinder"
c.name = "Cell markers"
ft.name = "Facet markers"


Then I define some problem specific parameters and the boundary conditions. Those are exactly the same as in the tutorial. I use a mixed function space.

############################################
#--- Parameters & boudnary conditions -----#
############################################

## We define our problem specific parameters
dt = Constant(mesh, PETSc.ScalarType(1/1600))
nu = Constant(mesh, PETSc.ScalarType(1/500)) # nu = 1/Re
relax = Constant(mesh, PETSc.ScalarType(1e-6)) # relaxation constant
U_0 = 2 # longitunal/lenghtwise velocity at coordinate (x,y)=(0,H/2)

# Inlet velocity profile
class Inlet_Velocity():
def __init__(self, U_0, H):
self.U_0 = U_0
self.H = H
def __call__(self, x):
values = np.zeros((gdim, x.shape), dtype=PETSc.ScalarType)
values = 4 * self.U_0 * x * (self.H - x)/(self.H**2)
return values

## As we have created the mesh and relevant mesh tags, we can now specify the function
## spaces V aand Q along with the boundary conditions. As the ft contains markers for facets,
## we use this class to find the facets for the inlet and walls.

Velocity_element = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Pressure_element = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, MixedElement([Velocity_element, Pressure_element]))
V, V_to_W = W.sub(0).collapse()
Q, Q_to_W = W.sub(1).collapse()

fdim = mesh.topology.dim - 1
# Inlet
u_inlet = Function(V)
inlet_velocity = Inlet_Velocity(U_0, H)
u_inlet.interpolate(inlet_velocity)
bcu_inflow = dirichletbc(u_inlet, locate_dofs_topological(V, fdim, ft.find(inlet_marker)))

# Walls
u_nonslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bcu_walls = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(wall_marker)), V)

# Obstacle
bcu_obstacle = dirichletbc(u_nonslip, locate_dofs_topological(V, fdim, ft.find(obstacle_marker)), V)

# grouping all boundary conditions involving the velocity u
bcu = [bcu_inflow, bcu_walls, bcu_obstacle]

# Outlet
p_outlet = PETSc.ScalarType(0)
bcp_outlet = dirichletbc(p_outlet, locate_dofs_topological(Q, fdim, ft.find(outlet_marker)), Q)

# grouping all boundary conditions involving the pressure p
bcp = [bcp_outlet]


bcs = bcu + bcp

Now I define the scheme as in equation (6) with BDF of order 1 (BDF 1). w_prev1 is the solution at step 0 (known) and w is the unknown of the problem (solution at step 1).

############################################
#--------- Discretization scheme ----------#
############################################

# defining the temporal scheme for the first time step (BDF1 is used, so BDF is set to 1)
v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)

w_prev1 = Function(W)
u_prev1, _ = split(w_prev1)

# discretization scheme
F = dot(u - u_prev1, v) * dx\
- p * div(v) * dt * dx\
- q * div(u) * dt * dx\
- relax * p * q * dt * dx


Finally I solve the system for only one time step.

############################################
#----------------- Solving ----------------#
############################################

# creating the problem and solving it
J = derivative(F, w)
problem = NonlinearProblem(F, w, bcs=bcu+bcp, J=J)
solver = NewtonSolver(mesh_comm, problem)
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()
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(w)
assert(converged)
print(f"Number of interations: {n:d}")


I would like to have a result like in Implementation — FEniCSx tutorial. Instead I have a core dump error, or this kind of result (no iteration at all):

2023-01-25 08:59:10.937 (   0.165s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = 0 (tol = 1e-10) r (rel) = -nan(tol = 1e-06)
2023-01-25 08:59:10.937 (   0.165s) [main            ]       NewtonSolver.cpp:255   INFO| Newton solver finished in 0 iterations and 0 linear solver iterations.
Number of interations: 0
2023-01-25 08:59:10.992 (   0.220s) [main            ]             loguru.cpp:526   INFO| atexit



I am to some extent a beginner in FEniCSx and solvers in general. Maybe I am not using the right solver because, to me, the weak formulation I use is linear as the convective term was linearized (am I right ?). In fact this way of solving incompressible Navier-Stokes equations has already been used in old FEniCS using NonLinearVariationalProblem and NonLinearVariationalSolver as follow:

J = derivative(F_BDF, w)
problem = NonlinearVariationalProblem(F_BDF, w, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.parameters["newton_solver"]["linear_solver"] = "mumps"
solver.solve()


And it worked, but now with FEniCSx it throws error. Can anyone help me findout why is the solver not converging in my case ? I hope I am clear in my explanations.

Thanks and regards.

There is a lot of code here, and I do not have time to parse all of it. The thing that stands out for me is

which I believe should be
u,p = ufl.split(w)

I would also set the linear solver to ‘preonly’ and ksp_type ‘lu’, see Implementation — FEniCSx tutorial

This means that your initial velocity guess 0, solves the problem, which likely indicates that you have an issue with your boundary conditions.
I would strongly recommend making sure that the correct parts of the domain is marked, and that the boundary condition is correct by first:

1. Write the facet markers to file, and check that they are sensible
2. Use dolfinx.fem.petsc.set_bc on a Function and write it to file, to check that the correct boundary condition is applied. (see for instance: Setting initial condition on irregular boundary from marked input mesh - #2 by dokken)

Thanks for your response. You are right, the problem seems to come from my boundary conditions. I have updated my code as you said:

First I checked if the facet markers were correctly assigned thanks to the adding of those lines:

with XDMFFile(mesh.comm, f"out_gmsh/mesh_rank_{model_rank}.xdmf", "w") as file:
file.write_mesh(mesh)
file.write_meshtags(ft)


As we can see in the figure Figure_1, the facet markers seem to be ok.

Next I imported set_bc from dolfinx.fem.petsc and used it like that:

u, p = w.split()
set_bc(u.vector, bcs)
set_bc(p.vector, bcs)

u.name = "Velocity"
p.name = "Pressure"
xdmf = XDMFFile(mesh_comm, "bc_test.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u)
xdmf.write_function(p)


This shows a problem we can see in Figure_2, the initial condition is chaotic yet not really random since the values range from 0 to 2 as expected by the inlet velocity profile (I set U_0=2).

When I change u, p = w.split() by u, p = split(w) as you suggested, I get this error (I guess because ufl.split() and fem.Function.split() are not the same):

AttributeError: 'ListTensor' object has no attribute 'vector'


Now if I code that:

u = Function(V)
p = Function(Q)
set_bc(u.vector, bcs)
set_bc(p.vector, bcs)

u.name = "Velocity"
p.name = "Pressure"
xdmf = XDMFFile(mesh_comm, "bc_test.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u)
xdmf.write_function(p)


I get the result seen in Figure_3. I check if all boundary conditions were sensible, and it is the case. So it works that way.

As I understand, the problem seems to come from the splitting method. I also wanted to add that I now get this error message when solving:

2023-01-25 19:21:21.123 (   0.154s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 0: r (abs) = 412.453 (tol = 1e-06) r (rel) = inf(tol = 1e-06)
2023-01-25 19:21:21.140 (   0.171s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-01-25 19:21:21.216 (   0.248s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 1: r (abs) = 5.33965e+06 (tol = 1e-06) r (rel) = 0.93322(tol = 1e-06)
2023-01-25 19:21:21.241 (   0.272s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-01-25 19:21:21.292 (   0.323s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 2: r (abs) = 1.86976e+11 (tol = 1e-06) r (rel) = 32678.1(tol = 1e-06)
2023-01-25 19:21:21.317 (   0.348s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-01-25 19:21:21.364 (   0.395s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 3: r (abs) = 7.38665e+15 (tol = 1e-06) r (rel) = 1.29098e+09(tol = 1e-06)
2023-01-25 19:21:21.388 (   0.419s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-01-25 19:21:21.435 (   0.467s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 4: r (abs) = 3.64202e+20 (tol = 1e-06) r (rel) = 6.36522e+13(tol = 1e-06)
2023-01-25 19:21:21.460 (   0.491s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-01-25 19:21:21.507 (   0.538s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 5: r (abs) = 3.1705e+25 (tol = 1e-06) r (rel) = 5.54114e+18(tol = 1e-06)
2023-01-25 19:21:21.532 (   0.563s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
...
2023-01-25 19:21:24.677 (   3.708s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 49: r (abs) = inf (tol = 1e-06) r (rel) = inf(tol = 1e-06)
2023-01-25 19:21:24.701 (   3.733s) [main            ]              petsc.cpp:677   INFO| PETSc Krylov solver starting to solve system.
2023-01-25 19:21:24.750 (   3.782s) [main            ]       NewtonSolver.cpp:36    INFO| Newton iteration 50: r (abs) = inf (tol = 1e-06) r (rel) = inf(tol = 1e-06)
Traceback (most recent call last):
File "/home/remi/Bureau/Scripts/Python/FEniCSx/Cases/Tests/first_step_clean.py", line 240, in <module>
n, converged = solver.solve(w)
File "/home/remi/anaconda3/envs/FEniCSx_env/lib/python3.10/site-packages/dolfinx/nls/petsc.py", line 41, in solve
n, converged = super().solve(u.vector)
RuntimeError: Newton solver did not converge because maximum number of iterations reached
2023-01-25 19:21:24.811 (   3.843s) [main            ]             loguru.cpp:526   INFO| atexit



Maybe the solver is not converging because of the chaotic boundary condition. Now the r (abs) and r (rel) are diverging toward inf. I believe it is better than before (before they were nan or 0…).

So, for the moment, it seems that the problem first lies in the definition and application of boundary conditions in a mixed element space.

I did not see any example of the usage of Newton Solver, for a problem with mixed elements, and Dirichlet BC, (maybe there are, I just didn’t find them) in FEniCSx. This kind of example could really help me to fix my problem.

You should set the bc on w, and then split w (with dolfinx.split() when saving data to xdmf, or using w.sub(0), w.sub(1) in outputting).

These boundary conditions are not correct for mixed spaces, see for instance
https://jsdokken.com/dolfinx-tutorial/chapter3/component_bc.html#boundary-conditions