Test problem 1: Channel flow (Poiseuille flow)

Authors: Anders Logg and Hans Petter Langtangen

Adapted to DOLFINx by: Jørgen S. Dokken

In this section, you will learn how to:

  • Solve the Navier-Stokes problem using a splitting scheme

  • Visualize functions from higher order Lagrangian spaces

In this section, we will compute the flow bewteen two infinite plates, so-called channel or Poiseuille flow. As we shall see, this problem has an analytical solution. Let \(H\) be the distance between the plates and \(L\) the length of the channel. There are no body forces.

We may scale the problem first to get rid of seemingly independent physical parameters. The physics of this problem is governed by viscous effects only, in the direction perpendicular to the flow, so a time scale should be based on diffusion accross the channel: \(t_v=H^2/\nu\). We let \(U\), some characteristic inflow velocity, be the velocity scale and \(H\) the spatial scale. The pressure scale is taken as the characteristic shear stress, \(\mu U/H\), since this is a primary example of shear flow. Inserting \(\bar{x}=x/H, \bar{y}=y/H, \bar{z}=z/H, \bar{u}=u/U, \bar{p}=Hp/{\mu U}\), and \(\bar{t}=H^2/\nu\) in the equations results in the scaled Navier-Stokes equations (dropping the bars after scaling)

(42)\[\begin{split}\frac{\partial u}{\partial t}+ \mathrm{Re} u \cdot \nabla u &= -\nabla p + \nabla^2 u,\\ \nabla \cdot u &=0.\end{split}\]

A detailed derivation for scaling of the Navier-Stokes equation for a large variety of physical situations can be found in [LP16] (Chapter 4.2) by Hans Petter Langtangen and Geir K. Pedersen.

Here, \(\mathrm{Re}=\rho UH/\mu\) is the Reynolds number. Because of the time and pressure scales, which are different from convection dominated fluid flow, the Reynolds number is associated with the convective term and not the viscity term.

The exact solution is derived by assuming \(u=(u_x(x,y,z),0,0)\) with the \(x\)-axis pointing along the channel. Since \(\nabla \cdot u = 0\), \(u\) cannot be dependent on \(x\).

The physics of channel flow is also two-dimensional so we can omit the \(z\)-coordinate (more precisely: \(\partial/\partial z = 0\)). Inserting \(u=(u_x, 0, 0)\) in the (scaled) governing equations gives \(u_x''(y)=\frac{\partial p}{\partial x)}\). Differentiating this equation with respect to \(x\) shows that \(\frac{\partial^2p}{\partial^2x}=0\) so \(\partial p/\partial x\) s a constant here called \(-\beta\). This is the driving force of the flow and can be specified as a known parameter in the problem. Integrating \(u_x''(x)=-\beta\) over the width of the channel, \([0,1]\), and requiring \(u=(0,0,0)\) at the channel walls, results in \(u_x=\frac{1}{2}\beta y(1-y)\). The characteristic inlet velocity \(U\) can be taken as the maximum inflow at \(y=0.5\), implying \(\beta=8\). The length of the channel, \(L/H\) in the scaled model, has no impact on the result, so for simplicity we just compute on the unit square. Mathematically, the pressure must be prescribed at a point, but since \(p\) does not depend on \(y\), we can set \(p\) to a known value, e.g. zero, along the outlet boundary \(x=1\). The result is \(p(x)=8(1-x)\) and \(u_x=4y(1-y)\).

The boundary conditions can be set as \(p=8\) at \(x=0\), \(p=0\) at \(x=1\) and \(u=(0,0,0)\) on the walls \(y=0.1\). This defines the pressure drop and should result in unit maximum velocity at the inlet and outlet and a parabolic velocity profile without no further specifications. Note that it is only meaningful to solve the Navier-Stokes equations i n 2D or 3D geometries, although the underlying mathematical problem collapses to two \(1D\) problems, one for \(u_x(y)\) and one for \(p(x)\).

The scaled model is not so easy to simulate using a standard Navier-Stokes solver with dimensions. However, one can argue that the convection term is zero, so the Re coefficient in front of this term in the scaled PDEs is not important and can be set to unity. In that case, setting \(\rho=\mu=1\) in the original Navier-Stokes equations resembles the scaled model.

For a specific engineering problem one wants to simulate a specific fluid and set corresponding parameters. A general solver is therefore most naturally implemented with dimensions and using the original physical parameters. However, scaling may greatly simplify numerical simulations. First of all, it shows that all fluids behave in the same way; it does not matter whether we have oil, gas, or water flowing betwen two plates, and it does not matter how fast the flow is (up to some critical value of the Reynolds number where the flow becomes unstable and transitions to a complicated turbulent flow of totally differnent nature.) This means that one simulation is enough to cover all types of channel flow! In other applications, scaling shows that it might be necessary to just set the fraction of some parameters (dimensionless numbers) rather than the parameters themselves. This simplifies exploring the input parameter space which is often the purpose of simulation. Frequently, the scaled problem is run by setting some of the input parameters with dimension to fixed values (often unity).

Implementation

Author: Jørgen S. Dokken

As in the previous example, we load the DOLFINx module, along with the mpi4py module, and create the unit square mesh and define the run-time and temporal discretization

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np

from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_matrix, 
                         assemble_vector, apply_lifting, create_vector, form, locate_dofs_geometrical,
                        set_bc)
from dolfinx.mesh import create_unit_square

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
t = 0
T = 10
num_steps = 500
dt = T/num_steps

As oposed to the previous demos, we will create our two function spaces using the ufl element definitions as input

import ufl
v_cg2 = ufl.VectorElement("CG", mesh.ufl_cell(), 2)
s_cg1 = ufl.FiniteElement("CG", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, v_cg2)
Q = FunctionSpace(mesh, s_cg1)

The first space V is a vector valued function space for the velocity, while Q is a scalar valued function space for pressure. We use piecewise quadratic elements for the velocity and piecewise linear elements for the pressure. When createing the vector finite element, the dimension of the vector element will be set to the geometric dimension of the mesh. One can easily create vector-valued function spaces with other dimensions by adding the keyword argument dim, i.e. v_cg = ufl.VectorElement("CG", mesh.ufl_cell(), 2, dim=10).

Stable finite element spaces for the Navier-Stokes equation

It is well-known that certain finite element spacces are not stable for the Navier-Stokes equations, or even for the simpler Stokes equation. The prime example of an unstable pair of finite element spacecs is to use first order degree continuous piecewise polynomials for both the velocity and the pressure. Using an unstable pair of spaces typically results in a solution with spurious (unwanted, non-physical) oscillations in the pressure solution. The simple remedy is to use continuous piecewise quadratic elements for the velocity and continuous piecewise linear elements for the pressure. Together, these elements form the so-called Taylor-Hood element. Spurious oscillations may occur also for splitting methids if an unstable element pair is used. } Since we have two different function spaces, we need to create two sets of trial and test functions:

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

As we have seen in Linear elasticty problem we can use Python-functions to create the different Dirichlet conditions. For this problem, we have three Dirichlet condition: First, we will set \(u=0\) at the walls of the channel, that is at \(y=0\) and \(y=1\). In this case, we will use dolfinx.fem.locate_dofs_geometrical

def walls(x):
    return np.logical_or(np.isclose(x[1],0), np.isclose(x[1],1))
wall_dofs = locate_dofs_geometrical(V, walls)
u_noslip = np.array((0,) * mesh.geometry.dim, dtype=PETSc.ScalarType)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

Second, we will set \(p=8\) at the inflow (\(x=0\))

def inflow(x):
    return np.isclose(x[0], 0)
inflow_dofs = locate_dofs_geometrical(Q, inflow)
bc_inflow = dirichletbc(PETSc.ScalarType(8), inflow_dofs, Q)

And finally, \(p=0\) at the outflow (\(x=1\)). This will result in a pressure gradient that will accelerate the flow from the initial state with zero velocity. At the end, we collect the boundary conditions for the velocity and pressure in Python lists so we can easily access them in the following computation.

def outflow(x):
    return np.isclose(x[0], 1)
outflow_dofs = locate_dofs_geometrical(Q, outflow)
bc_outflow = dirichletbc(PETSc.ScalarType(0), outflow_dofs, Q)
bcu = [bc_noslip]
bcp = [bc_inflow, bc_outflow]

We now move on to the definition of the three variational forms, one for each step in the IPCS scheme. Let us look at the definition of the first variational problem and the relevant parameters.

u_n = Function(V)
u_n.name = "u_n"
U = 0.5 * (u_n + u)
n = ufl.FacetNormal(mesh)
f = Constant(mesh, PETSc.ScalarType((0,0)))
k = Constant(mesh, PETSc.ScalarType(dt))
mu = Constant(mesh, PETSc.ScalarType(1))
rho = Constant(mesh, PETSc.ScalarType(1))

Usage of “dolfinx.Constant”

Note that we have wrapped several parameters as constants. This is to reduce the compilation-time of the variational formulations. By wrapping them as a constant, we can change the variable

The next step is to set up the variational form of the first step. As the variational problem contains a mix of known and unknown quantities, we will use the following naming convention: u is the known (mathematically \(u^{n+1}\)) as a trial function in the variational form. u_ is the most recently computed approximation (\(u^{n+1}\) available as a Function object), u_n is \(u^n\), and the same convetion goes for p,p_ (\(p^{n+1}\)) and p_n (p^n).

from ufl import Identity, div, dot, ds, dx, inner, lhs, nabla_grad, rhs, sym
# Define strain-rate tensor
def epsilon(u):
    return sym(nabla_grad(u))

# Define stress tensor
def sigma(u, p):
    return 2*mu*epsilon(u) - p*Identity(u.geometric_dimension())

# Define the variational problem for the first step
p_n = Function(Q)
p_n.name = "p_n"
F1 = rho*dot((u - u_n) / k, v)*dx
F1 += rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx
F1 += inner(sigma(U, p_n), epsilon(v))*dx
F1 += dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds
F1 -= dot(f, v)*dx
a1 = form(lhs(F1))
L1 = form(rhs(F1))

Note that we have used the ufl-functions lhs and rhs to sort out the bilinear form \(a(u,v)\) and linear form \(L(v)\). This is particulary convenient in longer and more complicated variational forms. With our particular discretization \(a(u,v)\) a1 is not time dependent, and only has to be assembled once, while the right hand side is dependent on the solution from the previous time step (u_n). Thus, we do as for the A known analytical solution, and create the matrix outside the time-loop.

A1 = assemble_matrix(a1, bcs=bcu)
A1.assemble()
b1 = create_vector(L1)

We now set up similar variational formulations and structures for the second and third step

# Define variational problem for step 2
u_ = Function(V)
a2 = form(dot(nabla_grad(p), nabla_grad(q))*dx)
L2 = form(dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx)
A2 = assemble_matrix(a2, bcs=bcp)
A2.assemble()
b2 = create_vector(L2)

# Define variational problem for step 3
p_ = Function(Q)
a3 = form(dot(u, v)*dx)
L3 = form(dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx)
A3 = assemble_matrix(a3)
A3.assemble()
b3 = create_vector(L3)

As we have create all the linear structures for the problem, we can now create a solver for each of them using PETSc. We can therefore customize the solution strategy for each step. For the tentative velocity step and pressure correction step, we will use the Stabilized version of BiConjugate Gradient to solve the linear system, and using algebraic myltigrid for preconditioning. For the last step, the velocity update, we use a conjugate gradient method with successive over relaxation, Gauss Seidel (SOR) preconditioning.

# Solver for step 1
solver1 = PETSc.KSP().create(mesh.comm)
solver1.setOperators(A1)
solver1.setType(PETSc.KSP.Type.BCGS)
pc1 = solver1.getPC()
pc1.setType(PETSc.PC.Type.HYPRE)
pc1.setHYPREType("boomeramg")

# Solver for step 2
solver2 = PETSc.KSP().create(mesh.comm)
solver2.setOperators(A2)
solver2.setType(PETSc.KSP.Type.BCGS)
pc2 = solver2.getPC()
pc2.setType(PETSc.PC.Type.HYPRE)
pc2.setHYPREType("boomeramg")

# Solver for step 3
solver3 = PETSc.KSP().create(mesh.comm)
solver3.setOperators(A3)
solver3.setType(PETSc.KSP.Type.CG)
pc3 = solver3.getPC()
pc3.setType(PETSc.PC.Type.SOR)

We prepare output files for the velocity and pressure data, and write the mesh and initial conditions to file

import dolfinx.io
xdmf = dolfinx.io.XDMFFile(MPI.COMM_WORLD, "poiseuille.xdmf", "w")
xdmf.write_mesh(mesh)
xdmf.write_function(u_n, t)
xdmf.write_function(p_n, t)

We also interpolate the analytical solution into our function-space and create a variational formulation for the \(L^2\)-error.

import numpy as np
def u_exact(x):
    values = np.zeros((2, x.shape[1]), dtype=PETSc.ScalarType)
    values[0] = 4*x[1]*(1.0 - x[1])
    return values
u_ex = Function(V)
u_ex.interpolate(u_exact)

L2_error = form(dot(u_ - u_ex, u_ - u_ex)*dx)

The next step is to create the loop over time. Note that we for all three steps only have to assemble the right hand side and apply the boundary condition using lifting.

for i in range(num_steps):
    # Update current time step
    t += dt

    # Step 1: Tentative veolcity step
    with b1.localForm() as loc_1:
        loc_1.set(0)
    assemble_vector(b1, L1)
    apply_lifting(b1, [a1], [bcu])
    b1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b1, bcu)
    solver1.solve(b1, u_.vector)
    u_.x.scatter_forward()
    
    # Step 2: Pressure corrrection step
    with b2.localForm() as loc_2:
        loc_2.set(0)
    assemble_vector(b2, L2)
    apply_lifting(b2, [a2], [bcp])
    b2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b2, bcp)
    solver2.solve(b2, p_.vector)
    p_.x.scatter_forward()

    # Step 3: Velocity correction step
    with b3.localForm() as loc_3:
        loc_3.set(0)
    assemble_vector(b3, L3)
    b3.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    solver3.solve(b3, u_.vector)
    u_.x.scatter_forward()
    # Update variable with solution form this time step
    u_n.x.array[:] = u_.x.array[:]
    p_n.x.array[:] = p_.x.array[:]

    # Write solutions to file
    xdmf.write_function(u_n, t)
    xdmf.write_function(p_n, t)

    # Compute error at current time-step
    error_L2 = np.sqrt(MPI.COMM_WORLD.allreduce(dolfinx.fem.assemble_scalar(L2_error), op=MPI.SUM))
    error_max = MPI.COMM_WORLD.allreduce(np.max(u_.vector.array - u_ex.vector.array), op=MPI.MAX)
    # Print error only every 20th step and at the last step
    if (i % 20 == 0) or (i == num_steps - 1):
        print(f"Time {t:.2f}, L2-error {error_L2:.2e}, Max error {error_max:.2e}")
# Close xmdf file
xdmf.close()
Time 0.02, L2-error 5.88e-01, Max error 1.60e-01
Time 0.42, L2-error 1.09e-02, Max error 1.28e-04
Time 0.82, L2-error 2.11e-04, Max error 2.65e-04
Time 1.22, L2-error 1.98e-05, Max error 1.56e-04
Time 1.62, L2-error 9.53e-06, Max error 8.46e-05
Time 2.02, L2-error 5.89e-06, Max error 5.24e-05
Time 2.42, L2-error 4.55e-06, Max error 3.62e-05
Time 2.82, L2-error 3.99e-06, Max error 2.70e-05
Time 3.22, L2-error 3.71e-06, Max error 2.13e-05
Time 3.62, L2-error 3.55e-06, Max error 1.75e-05
Time 4.02, L2-error 3.46e-06, Max error 1.49e-05
Time 4.42, L2-error 3.41e-06, Max error 1.30e-05
Time 4.82, L2-error 3.37e-06, Max error 1.16e-05
Time 5.22, L2-error 3.35e-06, Max error 1.06e-05
Time 5.62, L2-error 3.34e-06, Max error 9.75e-06
Time 6.02, L2-error 3.33e-06, Max error 9.11e-06
Time 6.42, L2-error 3.32e-06, Max error 8.92e-06
Time 6.82, L2-error 3.32e-06, Max error 8.92e-06
Time 7.22, L2-error 3.32e-06, Max error 9.05e-06
Time 7.62, L2-error 3.31e-06, Max error 9.26e-06
Time 8.02, L2-error 3.31e-06, Max error 9.44e-06
Time 8.42, L2-error 3.31e-06, Max error 9.58e-06
Time 8.82, L2-error 3.31e-06, Max error 9.69e-06
Time 9.22, L2-error 3.31e-06, Max error 9.79e-06
Time 9.62, L2-error 3.31e-06, Max error 9.86e-06
Time 10.00, L2-error 3.31e-06, Max error 1.05e-05

Verification

As for the previous problems we compute the error at each degree of freedom and the \(L^2(\Omega)\)-error. We start with the initial condition \(u=(0,0)\). We have not specified the initial condition explicitly, and FEniCSx will initialize all Functions including u_n and u_ to zero. Since the exact solution is quadratic, we expect to reach machine precision within finite time. For our implementation, we observe that the error quickly approaches zero, and is of order \(10^{-6}\) at $T=10

Visualization of higher order functions

As the velocity field \(u\) is a second order Lagrangian space, it’s degrees of freedom does not align with the mesh vertices. To be able to visualize this, we use the convenience function dolfinx.plot.create_vtk_topology. However, we will send in the function-space of the velocity field.

import dolfinx.plot
import pyvista
pyvista.set_jupyter_backend("pythreejs")

topology, cell_types = dolfinx.plot.create_vtk_topology(V)
# To make this function work in parallel, we only consider ghosts owned by the current process
num_dofs_local = V.dofmap.index_map.size_local
# We create a geometry for our modified mesh using the dof coordinates
geometry = V.tabulate_dof_coordinates()
# As we are dealing with a vector field, we reshape the underlying dof array to accommedate for the three dimensional space
num_dofs = V.dofmap.index_map.size_local + V.dofmap.index_map.num_ghosts
values = np.zeros((num_dofs, 3), dtype=np.float64)
values[:, :mesh.geometry.dim] = u_n.x.array.real.reshape(num_dofs, V.dofmap.index_map_bs)

# Create a point cloud of glyphs
function_grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
function_grid["u"] = values
glyphs = function_grid.glyph(orient="u", factor=0.2)

# Create a pyvista-grid for the mesh
topology_m, cell_types_m = dolfinx.plot.create_vtk_topology(mesh, mesh.topology.dim)
grid = pyvista.UnstructuredGrid(topology_m, cell_types_m, mesh.geometry.x)

# Create plotter
plotter = pyvista.Plotter()
plotter.add_mesh(grid, style="wireframe", color="k")
plotter.add_mesh(glyphs)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    pyvista.start_xvfb()
    fig_as_array = plotter.screenshot("glyphs.png")
WARNING:py.warnings:/usr/local/dolfinx-real/lib/python3.8/dist-packages/dolfinx/plot.py:120: UserWarning: Plotting of higher order functions is experimental.
  warnings.warn("Plotting of higher order functions is experimental.")

References

LP16

Hans Petter Langtangen and Geir K. Pedersen. Advanced partial differential equation models, pages 99–134. Springer International Publishing, Cham, 2016. doi:10.1007/978-3-319-32726-6_4.