KSPSolve error while solving steady Navier Stokes

General
1

I’m trying to solve the NS equation and get the following error:
" Failed to successfully call PETSc function ‘KSPSolve’. PETSc error code is: 77, Petsc has generated inconsistent data"

this is the MWC:

import dolfinx.fem.petsc

import numpy as np
import pyvista

from dolfinx import fem

import dolfinx.plot as dplt

from dolfinx.fem import (Constant,  VectorFunctionSpace, Function, FunctionSpace, assemble_scalar, 
                         dirichletbc, form, locate_dofs_topological,extract_function_spaces)

from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.mesh import create_unit_square, locate_entities, meshtags, create_rectangle, CellType, Mesh
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc

from ufl import (VectorElement,FiniteElement,FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, TrialFunctions, TestFunctions,
                 div, dot, dx, grad, inner, lhs, rhs, MixedElement, nabla_grad, split)
from dolfinx.io import XDMFFile
from dolfinx.plot import create_vtk_mesh

import matplotlib.pyplot as plt

from dolfinx import cpp as _cpp

from dolfinx.nls.petsc import NewtonSolver

def vel_in(x):
    val = -4*x[0]*x[0] + 4*x[0]
    return np.stack((np.zeros(x.shape[1]), -val*np.ones(x.shape[1])))

mesh = create_rectangle(comm=MPI.COMM_WORLD,
                            points=((0.0, 0.0), (1.0, 3.0)), n=(25, 75),
                            cell_type=CellType.triangle)

ds = Measure("ds", domain=mesh)

fdim = mesh.topology.dim - 1

P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

TH = MixedElement([P2, P1])

#TH = P2*P1
W = FunctionSpace(mesh, TH)

V, V_to_W = W.sub(0).collapse()

Q, Q_to_W = W.sub(1).collapse()

# parabolic inlet
u_in = Function(V)
u_in.interpolate(vel_in)
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[1], 3.0))
dofs = locate_dofs_topological(V, fdim, where)
bc0 = dirichletbc(u_in, dofs)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[0], 0.0))
dofs = locate_dofs_topological(V, fdim, where)
bc1 = dirichletbc(noslip, dofs, V)

# noslip
noslip = Constant(mesh,(0.0,0.0))
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[0], 1.0))
dofs = locate_dofs_topological(V, fdim, where)
bc2 = dirichletbc(noslip, dofs, V)

#outflow
outf = Constant(mesh,(0.0))
where = locate_entities(mesh, fdim,marker=lambda t: np.isclose(t[1], 0.0))
dofs = locate_dofs_topological(Q, fdim, where)
bc3 = dirichletbc(outf, dofs, Q)

bcs = [bc0,bc1,bc2,bc3]

up = Function(W)

u, p = split(up)

(v, q) = TestFunctions(W)

f = Constant(mesh, (ScalarType(0), ScalarType(0)))
g = Constant(mesh, (ScalarType(0)))

Re = Constant(mesh,ScalarType(10.0))

F= (1/Re)*inner(grad(u), grad(v)) * dx + inner(dot(u, nabla_grad(u)),v) * dx \
  - inner(div(v),p) * dx - inner(div(u),q) * dx + inner(f, v) * dx + inner(g,q) * ds

problem = NonlinearProblem(F, up, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()

opts[f"{option_prefix}ksp_type"] = "fgmres"
opts[f"{option_prefix}pc_type"] = "fieldsplit"
opts[f"{option_prefix}pc_fieldsplit_type"] = "schur"

# Set the direct solver types for each block
opts[f"{option_prefix}pc_fieldsplit_0_ksp_type"] = "preonly"
opts[f"{option_prefix}pc_fieldsplit_0_pc_type"] = "lu"  # or "mumps", "superlu"

opts[f"{option_prefix}pc_fieldsplit_1_ksp_type"] = "preonly"
opts[f"{option_prefix}pc_fieldsplit_1_pc_type"] = "jacobi"  # or "ilu"

ksp.setFromOptions()

r = solver.solve(up)

print(f"Num iterations: {r[0]}")
2

Start by making sure that normal solver works:

import dolfinx.fem.petsc
import dolfinx.plot as dplt
import matplotlib.pyplot as plt
import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace,
                         VectorFunctionSpace, assemble_scalar, dirichletbc,
                         extract_function_spaces, form,
                         locate_dofs_topological)
from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.io import VTXWriter, XDMFFile
from dolfinx.mesh import (CellType, Mesh, create_rectangle, create_unit_square,
                          locate_entities, meshtags)
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.plot import create_vtk_mesh
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FacetNormal, FiniteElement, Measure, MixedElement,
                 SpatialCoordinate, TestFunction, TestFunctions, TrialFunction,
                 TrialFunctions, VectorElement, div, dot, dx, grad, inner, lhs,
                 nabla_grad, rhs, split)

from dolfinx import cpp as _cpp
from dolfinx import fem


def vel_in(x):
    val = -4*x[0]*x[0] + 4*x[0]
    return np.stack((np.zeros(x.shape[1]), -val*np.ones(x.shape[1])))


mesh = create_rectangle(comm=MPI.COMM_WORLD,
                        points=((0.0, 0.0), (1.0, 3.0)), n=(25, 75),
                        cell_type=CellType.triangle)

ds = Measure("ds", domain=mesh)

fdim = mesh.topology.dim - 1

P2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

TH = MixedElement([P2, P1])

# TH = P2*P1
W = FunctionSpace(mesh, TH)

V, V_to_W = W.sub(0).collapse()

Q, Q_to_W = W.sub(1).collapse()

# parabolic inlet
u_in = Function(V)
u_in.interpolate(vel_in)
where = locate_entities(mesh, fdim, marker=lambda t: np.isclose(t[1], 3.0))
dofs = locate_dofs_topological((W.sub(0), V), fdim, where)
bc0 = dirichletbc(u_in, dofs, W.sub(0))

# noslip
noslip = dolfinx.fem.Function(V)
noslip.x.array[:] = 0
where = locate_entities(mesh, fdim, marker=lambda t: np.isclose(t[0], 0.0))
dofs = locate_dofs_topological((W.sub(0), V), fdim, where)
bc1 = dirichletbc(noslip, dofs, W.sub(0))

# noslip
where = locate_entities(mesh, fdim, marker=lambda t: np.isclose(t[0], 1.0))
dofs = locate_dofs_topological((W.sub(0), V), fdim, where)
bc2 = dirichletbc(noslip, dofs, W.sub(0))


bcs = [bc0, bc1, bc2]

up = Function(W)

u, p = split(up)

(v, q) = TestFunctions(W)

f = Constant(mesh, (ScalarType(0), ScalarType(0)))
g = Constant(mesh, (ScalarType(0)))

Re = Constant(mesh, ScalarType(10.0))

F = (1/Re)*inner(grad(u), grad(v)) * dx + inner(dot(u, nabla_grad(u)), v) * dx \
    - inner(div(v), p) * dx - inner(div(u), q) * \
    dx + inner(f, v) * dx + inner(g, q) * ds

problem = NonlinearProblem(F, up, bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)

solver.rtol = 1e-6

ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "gmres"
opts[f"{option_prefix}pc_type"] = "lu"

ksp.setFromOptions()

r = solver.solve(up)

print(f"Num iterations: {r[0]}")

u = up.sub(0).collapse()
with VTXWriter(mesh.comm, "u.bp", [u]) as vtx:
    vtx.write(0.0)

Here I’ve modified your boundary conditions to be applied correctly for a mixed space.