"PETSc KSPSolve error 73: Object is in wrong state - how to troubleshoot solver configuration issues for a nonlinear PDE?"

General
1

Hello everyone,

I am very new to numerical methods and programming with dolfinx. I am currently trying to solve a 2D nonlinear reaction-diffusion equation on a given mesh. The goal is to describe a tumor growth in a specific region of the brain (defined as “unhealthy” in my code). If I solve the equation on the whole mesh, the simulation works fine. Now I want to define the region so that the simulation is only performed in the “unhealthy” region. Therefore, I have created an ufl.Measure which I replace with dx in the variation formula F. When I do this, I get the following error message:

Traceback (most recent call last):
  File "/home/jz97x/Mesh_25_06/Reaction_Diffusion_2D_v1_mesh_var.py", line 162, in <module>
    r = solver.solve(u)
  File "/usr/lib/petsc/lib/python3/dist-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

I tried to play with different solver settings and boundary conditions and check similar problems in the forum. Unfortunately with no positive result. In short, my goal is to run the simulation within the “unhealthy” domain, which is a subdomain of the “healthy_unhealthy” mesh. Please find attached the actual code.

Any help would be greatly appreciated, thanks in advance and have a nice day.
Julian

from mpi4py import MPI
from petsc4py.PETSc import ScalarType
import numpy as np
import ufl
from dolfinx import fem, io, mesh, plot, default_real_type, default_scalar_type, log
from dolfinx.fem.petsc import LinearProblem, NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from dolfinx.io import XDMFFile, gmshio
from dolfinx.mesh import CellType, create_unit_square
from ufl import dx, grad, inner
import pyvista as pv
import pyvistaqt as pvqt
from dolfinx import *

#       ∂ u 
#       ─── - div(d ∇ u) = u * (K - u)
#       ∂ t

dt = 0.02
t = 0.0
T = 0.1

def custom_function(x):
    #print('Test')
    #print(x[0],x[1])
    dummy = np.zeros(x[0].shape)
    for index, (X, Y) in enumerate(zip(x[0],x[1])):
        if X <= 0.1 and X >= 0.08 and Y >= -0.07 and Y <= -0.05:  
            dummy[index] = 0.01 + X**10 + Y**5
        else:
            dummy[index] = 0 #1 + (X**5)/10 + (X**4)/5 + (X**3)
    return dummy

def initial_cond(x):
    # print('Test')
    # print(x[0],x[1])
    dummy = np.zeros(x[0].shape)
    for index, (X, Y) in enumerate(zip(x[0],x[1])): #replace with tumor region
        #if X <= 0.106 and X >= 0.098 and Y >= -0.102 and Y <= -0.095:
        if X <= 0.1 and X >= 0.08 and Y >= -0.07 and Y <= -0.05:    
            dummy[index] = 0.01 + X**2 + Y**2
        else:
            dummy[index] = 0
    return dummy

with XDMFFile(MPI.COMM_WORLD, "Domain.xdmf", "r") as xdmf:
    msh = xdmf.read_mesh(name="Grid")
    healthy_unhealthy = xdmf.read_meshtags(msh, name="Grid") 

with XDMFFile(MPI.COMM_WORLD, "Domain_boundary.xdmf", "r") as xdmf:
    boundary_msh = xdmf.read_mesh(name="Grid")
    boundary = xdmf.read_meshtags(boundary_msh, name="Grid") 

msh.topology.create_connectivity(msh.topology.dim, msh.topology.dim - 1)    

V = fem.functionspace(msh, ("Lagrange", 1)) #check this?? problem with dimensions??

healthy = healthy_unhealthy.find(1)
unhealthy = healthy_unhealthy.find(2)
outer = boundary.find(3)
inner1 = boundary.find(4)
inner2 = boundary.find(5)
tumor_bound = boundary.find(6)

unhealthy_marker = 2
dx_unhealthy = ufl.Measure('dx', domain=msh, subdomain_data=healthy_unhealthy, subdomain_id=unhealthy_marker)

print("Healthy subdomain:", healthy_unhealthy.find(1))
print("Unhealthy subdomain:", healthy_unhealthy.find(2))
print("Boundary markers:", boundary.values)

facets = mesh.locate_entities_boundary(msh, dim=(msh.topology.dim - 1),marker=lambda x: np.logical_or(np.isclose(x[0], 0.0),np.isclose(x[0], 2.0)))
#dofs = fem.locate_dofs_topological(V, dim, unhealthy)
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)

x = ufl.SpatialCoordinate(msh)

K_u = fem.Function(V)
K_u.interpolate(custom_function)

u = fem.Function(V)
u0 = fem.Function(V) 
v = ufl.TestFunction(V)

u.interpolate(initial_cond) 

d = 0.04

# linear
#F = (dt * d * inner(grad(u), grad(v)) + inner(u, v) - inner(u0, v)) * dx_unhealthy

#nonlinear
F = (dt * d * inner(grad(u), grad(v)) + dt * inner(inner(v, u), (K_u - u))) * dx_unhealthy + inner(u, v) * dx_unhealthy - inner(u0, v) * dx_unhealthy

problem = NonlinearProblem(F, u, bcs=[bc])
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental" #residual
solver.rtol = np.sqrt(np.finfo(default_real_type).eps) * 1e-2

file = XDMFFile(MPI.COMM_WORLD, "output_V1.xdmf", "w") 
file.write_mesh(msh)

u0.x.array[:] = u.x.array

while t < T:
    t += dt 
    r = solver.solve(u)
    print(f"Step {int(t / dt)}: num iterations: {r[0]}")
    u0.x.array[:] = u.x.array 
    file.write_function(u, t)
file.close()

import pyvista
cells, types, x = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(cells, types, x)
print(u.x.array.real)
grid.point_data["u"] = u.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
warped = grid.warp_by_scalar()
plotter.add_mesh(warped)
plotter.show()
2

I believe that you are running into this error because your mesh has two subdomains (healthy and unhealthy) but your problem is integrated only in unhealthy subdomain in F = (dt * d * inner(grad(u), grad(v)) + dt * inner(inner(v, u), (K_u - u))) * dx_unhealthy + inner(u, v) * dx_unhealthy - inner(u0, v) * dx_unhealthy. You will need to create a submesh if you are solving the problem only in one subdomain. Poisson submesh DG coupling · GitHub and mixed_domain_demos/hdg_poisson.py at main · jpdean/mixed_domain_demos · GitHub have examples of problems defined on submeshes.

1 Like