# Stationary Navier-Stokes equations

Hello,
I am recently new to using fenics. I would like to solve the stationary Navier-Stokes equations. I have found some posts from the legacy fenics, but it seems that in the new fenics it is not as straight-forward.

My starting point is a working example of the Stokes equation with an inlet and outlet like so

``````from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, io

import dolfinx.fem.petsc
import ufl

msh, cell_markers, facet_markers = io.gmshio.read_from_msh("mesh/DFG.msh", MPI.COMM_WORLD, gdim=2)

P2 = ufl.VectorElement("Lagrange", msh.topology.cell_types[0].name, degree=2)
P1 = ufl.FiniteElement("Lagrange", msh.topology.cell_types[0].name, degree=1)

V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)

u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)
v, q = ufl.TestFunction(V), ufl.TestFunction(Q)

a_up = ufl.inner(p, ufl.div(v)) * ufl.dx
a_pu = ufl.inner(ufl.div(u), q) * ufl.dx
a = fem.form([[a_uu, a_up], [a_pu, fem.Constant(msh, PETSc.ScalarType(0))*p*q*ufl.dx]])

f = fem.Constant(domain, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
L = fem.form([ufl.inner(f, v) * ufl.dx,
ufl.inner(fem.Constant(msh, PETSc.ScalarType(0)), q) * ufl.dx])

bc = []
tdim = msh.topology.dim

noslip = fem.Constant(domain, [0., 0.])
noslip_facet_id = 110  # Gmsh facet marker
noslip_facets = fem.locate_dofs_topological(V, tdim - 1, facet_markers.find(noslip_facet_id))
bc.append(fem.dirichletbc(noslip, noslip_facets, V))

inlet = fem.Constant(domain, [5., 0.])
inlet_facet_id = 111  # Gmsh facet marker
inlet_facets = fem.locate_dofs_topological(V, tdim - 1, facet_markers.find(inlet_facet_id))
bc.append(fem.dirichletbc(inlet, inlet_facets, V))

outlet = PETSc.ScalarType(0)
outlet_facet_id = 100  # Gmsh facet marker
outlet_facets = fem.locate_dofs_topological(Q, tdim - 1, facet_markers.find(outlet_facet_id))
bc.append(fem.dirichletbc(outlet, outlet_facets, Q))

A = fem.petsc.assemble_matrix_nest(a, bcs=bc)
A.assemble()

b = fem.petsc.assemble_vector_nest(L)
fem.petsc.apply_lifting_nest(b, a, bcs=bc)
for b_sub in b.getNestSubVecs():
spaces = fem.extract_function_spaces(L)
bcs0 = fem.bcs_by_block(spaces, bc)
fem.petsc.set_bc_nest(b, bcs0)

ksp = PETSc.KSP().create(comm)
ksp.setOperators(A)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")

pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

u, p = fem.Function(V), fem.Function(Q)
w = PETSc.Vec().createNest([u.vector, p.vector])
ksp.solve(b, w)
``````

Now for the stationary Navier-Stokes equation, if I simply add the convective nonlinear term to the bilinear form like this

``````a_uu = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(ufl.dot(u, ufl.grad(u)), v)*ufl.dx
``````

I get this error regarding to the creation of the bilinear form

``````Exception has occurred: TypeError
cannot unpack non-iterable bool object
``````

So I assume that the nonlinear term has to somehow be linearized, so maybe there has to be a linerization loop, but I am not quite sure how exactly.

please use the search button, there were posts about stationary NS equations in dolfinx in the last month or so, to which I (and others) replied.

I spent some time searching through the forum, but the only similar topic I could find is Convergence problem solving steady Navier-Stokes with multiphenicsx but even there the multiphenicsx library is used, so i was hoping somebody could provide some minimalistic example, just to get the basic idea of solving stationary NS.

Please note that you need to replace you trial function definition by a function definition in the case of a non-linear problem

Thank you guys for help, but I was thinking something more like this

``````u_old = fem.Function(V)

u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)
v, q = ufl.TestFunction(V), ufl.TestFunction(Q)

a_uu = ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(ufl.dot(u_old, ufl.grad(u)), v)*ufl.dx
a_up = ufl.inner(p, ufl.div(v)) * ufl.dx
a_pu = ufl.inner(ufl.div(u), q) * ufl.dx
a = fem.form([[a_uu, a_up], [a_pu, fem.Constant(msh, PETSc.ScalarType(0))*p*q*ufl.dx]])

A = dolfinx.fem.petsc.assemble_matrix_nest(a, bcs=bcs)
A.assemble()

ksp = PETSc.KSP().create(comm)
ksp.setOperators(A)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")

pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

niter = 20
for i in range(niter):
A = fem.petsc.assemble_matrix_nest(a, bcs=bc)
A.assemble()

b = fem.petsc.assemble_vector_nest(L)
fem.petsc.apply_lifting_nest(b, a, bcs=bc)
for b_sub in b.getNestSubVecs():
spaces = fem.extract_function_spaces(L)
bcs0 = fem.bcs_by_block(spaces, bc)
fem.petsc.set_bc_nest(b, bcs0)

w = PETSc.Vec().createNest([u.vector, p.vector])
ksp.solve(b, w)

rez = abs(u_old.x.norm() - u.x.norm())

u_old.x.array[:] = u.x.array[:]

if abs(norm_u_old - norm_u) < 1e-10:
break
``````

Is this a feasible method?

I wouldnt re-create A, w and b at each time step.
I would intialize them outside the loop, then Zero out their entries inside the opp and pass them to the assemblers.

I am sorry, but I dont think I follow, how can I assemble the matrix in the loop, when I zero out all entries.

See for instance the loop in: Test problem 2: Flow past a cylinder (DFG 2D-3 benchmark) — FEniCSx tutorial

Ok, I see. So my aim is to solve dfg_benchmark1_re20 . My current code looks like this

``````from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, io

import dolfinx.fem.petsc
import ufl

import numpy as np

msh, cell_markers, facet_markers = io.gmshio.read_from_msh("mesh/DFG.msh", MPI.COMM_WORLD, gdim=2)

# Parameters
Umax = 0.3
Re = 20
Lchar = 0.2
nu = 2/3.*Umax * Lchar / Re

P2 = ufl.VectorElement("Lagrange", msh.topology.cell_types[0].name, degree=2)
P1 = ufl.FiniteElement("Lagrange", msh.topology.cell_types[0].name, degree=1)

V, Q = fem.functionspace(msh, P2), fem.functionspace(msh, P1)

u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)
v, q = ufl.TestFunction(V), ufl.TestFunction(Q)
u_old = fem.Function(V)

a_uu = nu * ufl.inner(ufl.grad(u), ufl.grad(v)) * ufl.dx + ufl.dot(ufl.dot(u_old, ufl.grad(u)), v) * ufl.dx
a_up = ufl.inner(p, ufl.div(v)) * ufl.dx
a_pu = ufl.inner(ufl.div(u), q) * ufl.dx
a = fem.form([[a_uu, a_up], [a_pu, fem.Constant(msh, PETSc.ScalarType(0))*p*q*ufl.dx]])

f = fem.Constant(msh, (PETSc.ScalarType(0), PETSc.ScalarType(0)))
L = fem.form([ufl.inner(f, v) * ufl.dx,
ufl.inner(fem.Constant(msh, PETSc.ScalarType(0)), q) * ufl.dx])

bcu = []
tdim = msh.topology.dim

noslip = fem.Constant(msh, [0., 0.])
noslip_facet_id = 110  # Gmsh facet marker
noslip_facets = fem.locate_dofs_topological(V, tdim - 1, facet_markers.find(noslip_facet_id))
bcu.append(fem.dirichletbc(noslip, noslip_facets, V))

def inlet_parabolic_velocity_expression(x):
values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
values[0] = 4 * Umax * x[1] * (0.41 - x[1]) / (0.41**2)
return values

inlet = fem.Function(V)
inlet_facet_id = 111  # Gmsh facet marker
inlet_facets = fem.locate_dofs_topological(V, tdim - 1, facet_markers.find(inlet_facet_id))
inlet.interpolate(inlet_parabolic_velocity_expression)
bcu.append(fem.dirichletbc(inlet, inlet_facets))

outlet = PETSc.ScalarType(0)
outlet_facet_id = 100  # Gmsh facet marker
outlet_facets = fem.locate_dofs_topological(Q, tdim - 1, facet_markers.find(outlet_facet_id))
bcu.append(fem.dirichletbc(outlet, outlet_facets, Q))

A = fem.petsc.assemble_matrix_nest(a, bcs=bcu)
A.assemble()
b = fem.petsc.assemble_vector_nest(L)

ksp = PETSc.KSP().create(msh.comm)
ksp.setOperators(A)
ksp.setType("minres")
ksp.setTolerances(rtol=1e-9)
ksp.getPC().setType("fieldsplit")

pc = ksp.getPC()
pc.setType("lu")
pc.setFactorSolverType("mumps")
pc.setFactorSetUpSolverType()
pc.getFactorMatrix().setMumpsIcntl(icntl=24, ival=1)
pc.getFactorMatrix().setMumpsIcntl(icntl=25, ival=0)

u, p = fem.Function(V), fem.Function(Q)
w = PETSc.Vec().createNest([u.vector, p.vector])

niter = 30
for i in range(niter):
A.zeroEntries()
b.zeroEntries()
w.zeroEntries()

fem.petsc.assemble_matrix_nest(A, a, bcs=bcu)
A.assemble()

fem.petsc.assemble_vector_nest(b, L)
fem.petsc.apply_lifting_nest(b, a, bcs=bcu)
for b_sub in b.getNestSubVecs():
spaces = fem.extract_function_spaces(L)
bc0 = fem.bcs_by_block(spaces, bcu)
fem.petsc.set_bc_nest(b, bc0)

ksp.solve(b, w)

norm_u = u.x.norm()
norm_u_old = u_old.x.norm()
rez = abs(u_old.x.norm() - u.x.norm())

u_old.x.array[:] = u.x.array[:]

if rez < 1e-13:
break

from pathlib import Path
name = "DFG_navier_stokes"
folder = Path("results")
vtx_u = io.VTXWriter(msh.comm, folder / f"{name}_velocity_final.bp", u, engine="BP4")
vtx_p = io.VTXWriter(msh.comm, folder / f"{name}_pressure_final.bp", p, engine="BP4")
vtx_u.write(0)
vtx_p.write(0)
``````

But the results I get are not similar to the presented results in the article, but they remind me the solution for Stokes equation.

Hi, I referenced your code. But when the kinematic viscosity is very small (such as 1e-6), the result will not converge.

The Reynolds number is very large in that case, you’ll need to add some stabilization to your FE scheme. A possible acronym to search for is SUPG, but there is a whole research field of techniques working on this.