Lid Driven Cavity Flow: PETSc error 73, Object not found

Errors
1

I am trying to use FEniCSx to simulate lid driven cavity flow. When I try and solve the problem using a nonlinear solver, I am given the following error.

RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 73, Object is in wrong state

I am using a mixed function space (Taylor-Hood) elements to solve the problem, and I am given the above error, and I do not know what is causing it. I am using FEniCSx version 0.8, installed on linux using conda.

Here is my MWE

# The required modules are first imported:

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
import ufl
from basix.ufl import element, mixed_element
from dolfinx import fem
from dolfinx.fem import (Function, dirichletbc, functionspace, locate_dofs_topological)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import dx, grad, inner

# Create mesh
msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)

# Function to mark outer edges
def boundary(x):
    return np.isclose(x[0], 0.0) | np.isclose(x[0], 1.0) | np.isclose(x[1], 0.0) | np.isclose(x[1], 1.0)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,))
P1 = element("Lagrange", msh.basix_cell(), 1)
V, Q = functionspace(msh, P2), functionspace(msh, P1)

# No-slip condition on boundaries
noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType)  # type: ignore
facets = locate_entities_boundary(msh, 1, boundary)
bc0 = dirichletbc(noslip, locate_dofs_topological(V, 1, facets), V)

# Create the Taylot-Hood function space
TH = mixed_element([P2, P1])
W = functionspace(msh, TH)

# Boundary condition
W0 = W.sub(0)
Q, _ = W0.collapse()
noslip = Function(Q)
facets = locate_entities_boundary(msh, 1, boundary)
dofs = locate_dofs_topological((W0, Q), 1, facets)
bc0 = dirichletbc(noslip, dofs, W0)

uh = fem.Function(W)
u, p = uh.split()
(v, q) = ufl.TestFunctions(W)

# Define the residual, just the viscous term is needed to produce the error
R = inner(grad(u),grad(v))*dx

Problem = NonlinearProblem(R, uh, [bc0]) # Formulate nonlinear problem

solver = NewtonSolver(MPI.COMM_WORLD, Problem) # Uses a newton solver
ksp = solver.krylov_solver
opts = PETSc.Options()
solver.convergence_criterion = "incremental"
solver.solve(uh) # Error occurs on this line

Note, for my MWE I am using only one part of the Navier-Stokes equations, and I have simplified the boundary conditions to be zero along all 4 edges

2

Hi, welcome to the FEniCS community!.

Please note that when using mixed element, you need to split the function with ufl.split(fun) to properly write the variational formulation.

u, p = ufl.split(uh)

This have been discussed already in, for example, Picard Iteration Not Updating Properly - #2 by kamensky .

Cheers.

1 Like
3

Thank you for the help, I got my code working.