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()
```