Navier stokes lid driven cavity flow

[alex24]

Hi,

I am trying to run this older example of a lid driven cavity

Error in Solving Steady Navier Stokes Equation

I’m trying to convert my legacy FEniCS codes to FEniCSx to solve steady incompressible NSE with the Newton method for a Lid-Driven Cavity. I’m not able to set up the nonlinear solver properly. The FEniCS code is: # Define function spaces (Taylor-Hood) V = VectorElement("Lagrange", mesh.ufl_cell(), 2) Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1) # Make a mixed space TH = V * Q W = FunctionSpace(mesh, TH) # Define unknown and test function(s) (v, q) = TestFunctions(W) w = Function(W) (u, …

I was able to replace the element definitions by using basix:

# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

but I can’t find how to create a function space by combining elements in the current version of fenicsx.

This line:
TH = P2 * P1
Generates this error:
TypeError: unsupported operand type(s) for *: '_BlockedElement' and '_BasixElement'

Thank you,
Alex

The code I’m trying to update is here:

from dolfinx import mesh, fem, nls
from mpi4py import MPI
import numpy as np
import ufl
from ufl import inner, div, grad, dx
from petsc4py import PETSc

import basix

msh = mesh.create_unit_square(MPI.COMM_WORLD, 40, 40, mesh.CellType.triangle)

nu = 0.001

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity#
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

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

# Create the function space
TH = P2 * P1
W = fem.FunctionSpace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1 = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.set(0.0)
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step
F = inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v)) * \
    dx - div(v)*p*dx - q*div(u)*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = fem.petsc.NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
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()

# Compute the solution
solver.solve(w)

# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()

[alex24]

Ok I solved it using basix.ufl.mixed_element

[dokken]

See

Introduction to the Unified Form Language — FEniCS Workshop

which instructs you that you can use basix.ufl.mixed_element

[alex24]

I just found it thank you!

[alex24]

The other issue was this line:

zero.x.set(0.0)

I solved it by using:

zero.x.array[:] = 0

From:

Sudden error in boundary conditions

Hi everyone, I’ve been working on this file for quite some time now without any error, using dolfinx version 0.6, in which I set the value of a function in a function space the following way: import numpy as np import ufl from dolfinx import * import dolfinx.fem.petsc import dolfinx from mpi4py import MPI from dolfinx import fem, mesh, plot from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, dirichletbc, form, locate_dofs_geometrical, set_bc) …

[alex24]

Now the code runs but it doesn’t converge,
RuntimeError: Newton solver did not converge because maximum number of iterations reached

The viscosity and all other parameters are the same as the original code, the latest code is the following:

from dolfinx import mesh, fem, nls
from mpi4py import MPI
import numpy as np
import ufl
from ufl import inner, div, grad, dx
from petsc4py import PETSc

import basix

msh = mesh.create_unit_square(MPI.COMM_WORLD, 40, 40, mesh.CellType.triangle)

nu = 0.001

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity#
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

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

# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs   = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1    = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step
F = inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v)) * \
    dx - div(v)*p*dx - q*div(u)*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = fem.petsc.NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
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()

# Compute the solution
solver.solve(w)

# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()

[alex24]

Ok, its working now, had to make the (kinematic?) viscosity ten times higher, and took out the
PETSc options. For the xdmf files, is it possible to store both velocity and pressure in the same xdmf file using XDMFFile? How can you name them so they show the names in ParaView?

The full running code is here:

from dolfinx import mesh, fem, io
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver

from mpi4py import MPI

import numpy as np

import ufl
from ufl import inner, div, grad, dx

import basix



msh = mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, mesh.CellType.triangle)

nu = 0.01

# Function to mark x = 0, x = 1 and y = 0
def noslip_boundary(x):
    return np.logical_or(np.logical_or(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0)), np.isclose(x[1], 0.0))

# Function to mark the lid (y = 1)
def lid(x):
    return np.isclose(x[1], 1.0)

# Lid velocity#
def lid_velocity_expression(x):
    return np.stack((np.ones(x.shape[1]), np.zeros(x.shape[1])))


# the same as old ufl.VectorElement
# the same as old ufl.FiniteElement
P2 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 2, shape=(msh.geometry.dim, ))
P1 = basix.ufl.element("Lagrange", msh.topology.cell_name(), 1)

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

# Create the function space
TH = basix.ufl.mixed_element([P2, P1])
W = fem.functionspace(msh, TH)
W0, _ = W.sub(0).collapse()

# No slip boundary condition
noslip = fem.Function(V)
facets = mesh.locate_entities_boundary(msh, 1, noslip_boundary)
dofs = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc0 = fem.dirichletbc(noslip, dofs, W.sub(0))

# Driving velocity condition u = (1, 0) on top boundary (y = 1)
lid_velocity = fem.Function(W0)
lid_velocity.interpolate(lid_velocity_expression)
facets = mesh.locate_entities_boundary(msh, 1, lid)
dofs   = fem.locate_dofs_topological((W.sub(0), V), 1, facets)
bc1    = fem.dirichletbc(lid_velocity, dofs, W.sub(0))

# Since for this problem the pressure is only determined up to a constant, we pin the pressure at the point (0, 0)
zero = fem.Function(Q)
zero.x.array[:] = 0
dofs = fem.locate_dofs_geometrical(
    (W.sub(1), Q), lambda x: np.isclose(x.T, [0, 0, 0]).all(axis=1))
bc2 = fem.dirichletbc(zero, dofs, W.sub(1))

# Collect Dirichlet boundary conditions
bcs = [bc0, bc1, bc2]

# Define variational problem
w = fem.Function(W)

(u, p) = ufl.split(w)
(v, q) = ufl.TestFunctions(W)

# Tentative velocity step
F = inner(grad(u)*u, v)*dx + nu*inner(grad(u), grad(v)) * \
    dx - div(v)*p*dx - q*div(u)*dx

dw = ufl.TrialFunction(W)
dF = ufl.derivative(F, w, dw)

problem = NonlinearProblem(F, w, bcs=bcs, J=dF)

solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol   = 1e-6
solver.report = True
solver.error_on_nonconvergence = True
#solver.max_it = 10


# Compute the solution
solver.solve(w)


# Split the mixed solution and collapse
u = w.sub(0).collapse()
p = w.sub(1).collapse()


# Write the solution to file
with io.XDMFFile(MPI.COMM_WORLD, "pressure.xdmf", "w") as pfile_xdmf:
    p.x.scatter_forward()
    pfile_xdmf.write_mesh(msh)
    pfile_xdmf.write_function(p)

with io.XDMFFile(MPI.COMM_WORLD, "velocity.xdmf", "w") as ufile_xdmf:
    u.x.scatter_forward()
    P1 = basix.ufl.element("Lagrange", msh.basix_cell(), 1, shape=(msh.geometry.dim,))
    u1 = fem.Function(fem.functionspace(msh, P1))
    u1.interpolate(u)
    ufile_xdmf.write_mesh(msh)
    ufile_xdmf.write_function(u1)