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

Errors
1

i am new to fenicsx and Attempting to formulate the problem and solve the nonlinear system results in the following Error code 73 in PETSc. i have tried few provided solutions in discourse but still gets the same error.
below i am attaching my code. can you please help me with error.
additional info:
fenicsx version 0.8.0
conda

from mpi4py import MPI
from petsc4py import PETSc
import numpy as np
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, io, mesh, default_scalar_type,log
import ufl
from dolfinx.io import XDMFFile
import h5py
from math import pi
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.fem import functionspace, Function, locate_dofs_geometrical, dirichletbc, Constant

#Initialize MPI
comm = MPI.COMM_WORLD

#Load the mesh from the XDMF file
with io.XDMFFile(MPI.COMM_WORLD, "Circle_hole_small.xdmf", "r") as xdmf:
    msh = xdmf.read_mesh(name ='mesh')
    #write function to file
with io.XDMFFile(comm,'Circle_hole_small.xdmf','w') as ofile:
    ofile.write_mesh(msh)
    encoding=io.XDMFFile.Encoding.HDF5

    #Open the HDF5 file for reading
with h5py.File("Circle_hole_small.h5", "r") as ifile:
    # Read the mesh dataset from HDF5
    mesh_dataset = ifile["/Mesh"]

#Define measures
dx = ufl.Measure("dx", domain=msh)

p_u = element("Lagrange", msh.topology.cell_name(), 1, shape=(msh.topology.dim, ),dtype=default_real_type )
p_p = element("Lagrange", msh.topology.cell_name(), 1, dtype=default_real_type)
W_elem = mixed_element([p_u, p_p])

#Create mixed function space

W = fem.functionspace(msh, W_elem)

#Define Dirichlet boundary conditions
def outer_boundary(x):
    r = np.sqrt(x[0]**2 + x[1]**2)
    return np.isclose(r, 1.0)

def inner_boundary(x):
    r = np.sqrt(x[0]**2 + x[1]**2)
    return np.isclose(r, 0.001)
#Locate boundary facets
fdim = msh.topology.dim - 1
outer_facets = mesh.locate_entities_boundary(msh, fdim, outer_boundary)
inner_facets = mesh.locate_entities_boundary(msh, fdim, inner_boundary)

marked_facets = np.hstack([outer_facets, inner_facets])
marked_values = np.hstack([np.full_like(outer_facets, 2), np.full_like(inner_facets, 3)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(msh, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])

ds = ufl.Measure("ds", domain=msh, subdomain_data=facet_tag)
 
W1 = W.sub(1)
Q1,_ = W1.collapse()
p_bci = Function(Q1)
p_bci.x.array[:] = 0.2
p_bco = Function(Q1)
p_bco.x.array[:] = 0.0

dofs_p_out = fem.locate_dofs_topological((W1,Q1), fdim, outer_facets)
dofs_p_in = fem.locate_dofs_topological((W1,Q1), fdim, inner_facets)
outer_boundary_pressure = (fem.dirichletbc(p_bco, dofs_p_out, (W1)))
inner_boundary_pressure = (fem.dirichletbc(p_bci, dofs_p_in, (W1)))
bcs = [inner_boundary_pressure,outer_boundary_pressure]

E =  0.770
nu =  0.4 / 1.4
G = (E / (2 * (1 + nu)))

#Define strain and stress
def epsilon(u):
    SS = ufl.sym(ufl.grad(u))
    return SS

def sigma(u, p):
    
    """Return an expression for the stress σ given a displacement field"""
    #S = 2.0 * G * ufl.sym(ufl.grad(u)) + λ * ufl.tr(ufl.sym(ufl.grad(u))) * ufl.Identity(len(u)) - p * ufl.Identity(len(u))
    S = 2*G*ufl.sym(ufl.grad(u))+(2*G*nu)/(1-2*nu)*ufl.tr(ufl.sym(ufl.grad(u)))* ufl.Identity(len(u))-p* ufl.Identity(len(u))
    return S

w_coupled = fem.Function(W)  
(u, p)= ufl.split(w_coupled)
(q, v) = ufl.TestFunctions(W)
# Compute stress tensor
Stress = sigma(u, p)

def penalty1(u):
    return ufl.sqrt((u[0])**2 + (u[1])**2)

#Momentum balance
B = Constant(msh,default_scalar_type((0, 0)))
T = Constant(msh,default_scalar_type(0))
n = ufl.FacetNormal(msh)
Tr = T *n
Contact_contribution = pen*ufl.dot(q, (u))*ds(2)
L1 = (ufl.inner(Stress, ufl.grad(q)) * dx) - (ufl.inner(B, q) * dx) + Contact_contribution
#Add L of fluid problem
f = Constant(msh,0.0)
k = Constant(msh,1.0)
pressure_drop = 0.33
r1 =0.001
r2 = 1.0
g1 = Constant(msh,(-pressure_drop / (r1 * (np.log(r1) - np.log(r2)))))
A1 = 2.0 * pi * 0.001
A2 = 2.0 * pi * 1.0
g2 = -g1 * (A1 / A2)  # g2 --> outer boundary (2)
L2 = f * v * dx + g1 * v * ds(3) - k * ufl.inner(ufl.grad(p), ufl.grad(v)) * dx 

#State full problem

F = L1 + L2

#Solve variational problem
w_coupled = fem.Function(W)
my_problem = NonlinearProblem(F,w_coupled,bcs)
solver =  NewtonSolver(MPI.COMM_WORLD, my_problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.report = True
log.set_log_level(log.LogLevel.INFO)
n, converged = solver.solve(w_coupled)
assert (converged)
print(f"Number of interations: {n:d}")

when i am trying to run this code i am getting error of

Traceback (most recent call last):
  File "/home/user/workspace/code/try /trial.py", line 198, in <module>
    n, converged = solver.solve(w_coupled)
                   ^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/user/anaconda3/envs/fenicsx/lib/python3.12/site-packages/dolfinx/nls/petsc.py", line 47, in solve
    n, converged = super().solve(u.x.petsc_vec)
                   ^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: Failed to successfully call PETSc function 'KSPSolve'. PETSc error code is: 73, Object is in wrong state
2

Please read Read before posting: How do I get my question answered?, and please post the code within “```”. I have edited your post because the code as you posted it was unreadable. Please edit it back posting the properly formatted code, and a code which is reproducible (i.e., anyone can copy and paste it, and run it).

3

You are redefining w_coupled here, meaning that the automatically computed Jacobian is singular. Please remove this definition.

Secondly, please note that we cannot run your code as you haven’t provided a mesh.

1 Like