Stationary Navier-Stokes equations

[marhatom]

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_uu = ufl.inner(ufl.grad(u), ufl.grad(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(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():
     b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
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")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

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.

[francesco-ballarin]

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.

[marhatom]

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.

[francesco-ballarin]

Newton iteration against time-stepping for steady solution of the incompressible Navier-Stokes equations

Hi! I’m interested in using FEniCSx to compute the steady solution of incompressible Navier-Stokes equations. I’ve found several examples (in the tutorials or in this forum) but they are all based on time-stepping. But I’m coming from the world of FreeFem++ where we usually rather use a Newton iteration (as explained in section 2.1 of this article and shown in the code here). So here comes my question: is there such an example available for FEniCSx that I missed in my research? If not then w…

[dokken]

@marhatom

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

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

@marhatom

u, p = ufl.TrialFunction(V), ufl.TrialFunction(Q)

[marhatom]

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")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

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():
             b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        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?

[dokken]

@marhatom

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():
             b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
        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

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.

[marhatom]

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

[dokken]

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

[marhatom]

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")
ksp.getPC().setFieldSplitType(PETSc.PC.CompositeType.ADDITIVE)

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():
        b_sub.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE)
    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.

[Yubin_Huang]

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

[francesco-ballarin]

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.