Legacy FEniCS works, DOLFINx doesn't in Stokes channel flow

dolfinx
1

I’m really at a loss here. I have two versions of what is essentially the same Stokes channel flow script, one written in the new API, one in the legacy API. Both run, but the legacy produces correct results whereas the other one does not. I just cannot find the error and I would really appreciate some help.

This is the DOLFINx script:

import numpy as np
from dolfinx import io
from dolfinx.fem import (
    Constant,
    Function,
    FunctionSpace,
    dirichletbc,
    locate_dofs_geometrical,
    petsc,
)
from dolfinx.mesh import create_rectangle
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (
    FiniteElement,
    TestFunctions,
    TrialFunctions,
    VectorElement,
    div,
    dx,
    grad,
    inner,
)

# Numerical parameters
nx, ny = 90, 30

# Geometry definition
L = 3
H = 1
domain = create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([L, H])],
    [nx, ny],
)

# Physical parameters
mu = Constant(domain, 1.0)
f = Constant(domain, (0.0, 0.0))

# Define the function spaces
V_el = VectorElement("P", domain.ufl_cell(), 2)
Q_el = FiniteElement("DG", domain.ufl_cell(), 0)
W = FunctionSpace(domain, V_el * Q_el)
V, _ = W.sub(0).collapse()


# Define the boundary conditions
def inlet(x):
    return np.isclose(x[0], 0)


def outlet(x):
    return np.isclose(x[0], L)


def walls(x):
    # Top wall
    loc_top = np.isclose(x[1], H)
    # Bottom wall
    loc_bot = np.isclose(x[1], 0)
    return np.any([loc_top, loc_bot], axis=0)


def inlet_v(x):
    return np.vstack([-4 * x[1] * (x[1] - 1), np.zeros(x.shape[1])])


u_inlet = Function(V)
u_inlet.interpolate(inlet_v)
inlet_dofs = locate_dofs_geometrical(V, inlet)
bc_inlet = dirichletbc(u_inlet, inlet_dofs)

u_noslip = np.array([0.0, 0.0], dtype=PETSc.ScalarType)
wall_dofs = locate_dofs_geometrical(V, walls)
bc_noslip = dirichletbc(u_noslip, wall_dofs, V)

bcs = [bc_inlet, bc_noslip]

# Forming and solving the linear system
u, p = TrialFunctions(W)
v, q = TestFunctions(W)

a = mu * inner(grad(u), grad(v)) * dx - p * div(v) * dx + q * div(u) * dx
L = inner(v, f) * dx

problem = petsc.LinearProblem(a, L, bcs=bcs)
wh = problem.solve()
uh, ph = wh.split()
uh.name = "velocity"
ph.name = "pressure"

# External post-processing
with io.XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)
    xdmf.write_function(ph)

This is the legacy script (taken from a course on finite elements).

from __future__ import print_function
from dolfin import *
from fenics import *
from mshr import *
import numpy as np

# Problem parameters
    # Fluid parameters
mu = Constant(1)            # dynamic viscosity
    # External force
f  = Constant((0,0))

# Geometry setup and meshing
    # Create domain
domain = Rectangle(Point(0, 0), Point(3, 1))
    # Specify boundaries
def wall(x, on_boundary):
    return on_boundary and (near(x[1], 0) or near(x[1],1))
def inlet(x, on_boundary):
    return on_boundary and near(x[0], 0)
def outlet(x, on_boundary):
    return on_boundary and near(x[0], 3)
    # Create mesh-triangulation
mesh = generate_mesh(domain, 64)
plot(mesh, title='Mesh')

# Set up discrete setting
    # Finite Element degree
deg_velocity = 2
deg_pressure = 0
    # Define Elements
V_el = VectorElement("P", mesh.ufl_cell(), deg_velocity)
Q_el = FiniteElement("DG", mesh.ufl_cell(), deg_pressure)
W_el = MixedElement([V_el, Q_el])
    # Define function spaces
V = VectorFunctionSpace(mesh, "P", deg_velocity)
Q = FunctionSpace(mesh, "DG", deg_pressure)
W = FunctionSpace(mesh, W_el)
    # Define trial and test functions
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)
    # Define solution (functions)
w_sol = Function(W)
u_sol = Function(V)
p_sol = Function(Q)

# Set up boundary data
    # Define boundary values
bv_inlet  = Expression(('-4*(x[1]-0)*(x[1]-1)','0'), degree=2)
bv_wall   = Constant((0,0))
    # Define boundary conditions
bc_inlet  = DirichletBC(W.sub(0), bv_inlet,  inlet)
bc_wall   = DirichletBC(W.sub(0), bv_wall,   wall)

bcs = [bc_inlet, bc_wall]

# Define variational problem
F = ( mu*inner(grad(u), grad(v)) - div(v)*p + q*div(u) - inner(v, f))*dx
a = lhs(F)
L = rhs(F)

# Solve linear system
solve(a==L, w_sol, bcs)

# Split the mixed solution
(u_sol, p_sol) = w_sol.split()

# Plot solution
plot(u_sol, title='Velocity')
plot(p_sol, title='Pressure')
u_sol.rename('velocity', 'u')
p_sol.rename('pressure', 'p')
vtkfile = File('Stokes_velocity.pvd')
vtkfile << u_sol
vtkfile = File('Stokes_pressure.pvd')
vtkfile << p_sol

fenicsx_solution
fenicsx_solution1512×804 169 KB

I can only post one media file, but the legacy version produces the correct results. I’m running both programs on their respective Docker containers.

All help is appreciated!

2

The solution on the legacy FEniCS is attached here. Even though it uses a different grid, I’ve used a similar grid on the DOLFINx version, but the solution stayed the same.

legacy_solution
legacy_solution1512×804 189 KB

3

You have not set up your boundary conditions correctly, see:

import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace, dirichletbc,
                         locate_dofs_geometrical, locate_dofs_topological, petsc)
from dolfinx.mesh import create_rectangle, locate_entities_boundary
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (FiniteElement, MixedElement, TestFunctions, TrialFunctions,
                 VectorElement, div, dx, grad, inner)

from dolfinx import io

# Numerical parameters
nx, ny = 90, 30

# Geometry definition
L = 3
H = 1
domain = create_rectangle(
    MPI.COMM_WORLD,
    [np.array([0, 0]), np.array([L, H])],
    [nx, ny],
)

# Physical parameters
mu = Constant(domain, 1.0)
f = Constant(domain, (0.0, 0.0))

# Define the function spaces
V_el = VectorElement("Lagrange", domain.ufl_cell(), 2)
Q_el = FiniteElement("DG", domain.ufl_cell(), 0)
W = FunctionSpace(domain, MixedElement([V_el, Q_el]))
V, _ = W.sub(0).collapse()


# Define the boundary conditions
def inlet(x):
    return np.isclose(x[0], 0)


def outlet(x):
    return np.isclose(x[0], L)


def walls(x):
    # Top wall
    loc_top = np.isclose(x[1], H)
    # Bottom wall
    loc_bot = np.isclose(x[1], 0)
    return np.any([loc_top, loc_bot], axis=0)


def inlet_v(x):
    return (-4 * x[1] * (x[1] - 1), np.zeros(x.shape[1]))


u_inlet = Function(V)
u_inlet.interpolate(inlet_v)
inlet_dofs = locate_dofs_geometrical((W.sub(0), V), inlet)
bc_inlet = dirichletbc(u_inlet, inlet_dofs, W.sub(0))


def zero(x):
    return np.zeros((domain.topology.dim, x.shape[1]), dtype=PETSc.ScalarType)


u_walls = Function(V)
u_walls.interpolate(zero)
wall_facets = locate_entities_boundary(domain, domain.topology.dim - 1, walls)
wall_dofs = locate_dofs_topological((W.sub(0), V), domain.topology.dim - 1, wall_facets)
bc_noslip = dirichletbc(u_walls, wall_dofs, W.sub(0))

bcs = [bc_inlet, bc_noslip]

# Forming and solving the linear system
u, p = TrialFunctions(W)
v, q = TestFunctions(W)

a = mu * inner(grad(u), grad(v)) * dx - p * div(v) * dx + q * div(u) * dx
L = inner(v, f) * dx

# Use iterative solvers as pressure is not pinned at a single point or set the nullspace of the solver
problem = petsc.LinearProblem(a, L, bcs=bcs)
wh = problem.solve()
uh, ph = wh.split()
uh.name = "velocity"
ph.name = "pressure"

# External post-processing
with io.XDMFFile(domain.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)
    xdmf.write_function(ph)
4

Thank you so much @dokken! I will have a good look at the changes you suggested to better understand how to correctly call the API.