# Time-dependent nonlinear reaction-diffusion not proceeding in time loop

This is my first time posting (and I am brand new to FeniCSx), so hopefully I am doing it correctly. I have a nonlinear reaction-diffusion problem that I am trying to solve. It involves diffusion from a sphere of homogeneous concentration in a large cube, with a “reacting” sphere some distance d away from the center of the source (the origin in 3D). Here is a shortened working example. For some reason, the time loop does not proceed past the first time step. I am not sure where I am making an error. Any help is greatly appreciated!

``````import numpy as np
from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import fem, mesh, io, plot, nls
import ufl

# Define initial condition
def initial_condition(x):
return 1000.0 * (x[0]**2 + x[1]**2 + x[2]**2 <= 0.250)

# Parameters for reaction-diffusion
D = 800 # µm^2/s
kp = 57
Km = 0.001

#####################################
# -----------------------------------------------------------------------------------------------
#####################################

def rxn_dif(d,rho):

# Define temporal parameters
t = 0 # Start time
T = 10.0 # Final time
num_steps = 500
dt = T / num_steps # time step size

# ------------------------------------------------------------------------------------------------

# MESH

nx, ny, nz = 50, 50, 50
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([-10,-10,-10]), np.array([10, 10, 10])],
[nx, ny, nz], mesh.CellType.tetrahedron)

V = fem.FunctionSpace(domain, ("CG", 2))

# ------------------------------------------------------------------------------------------------

# INITIAL CONDITION

u_n = fem.Function(V)
u_n.name = "u_n"
u_n.interpolate(initial_condition)

# ------------------------------------------------------------------------------------------------

# BOUNDARY CONDITIONS

fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(
domain, fdim, lambda x: np.full(x.shape[1], True, dtype=bool))
bc = fem.dirichletbc(PETSc.ScalarType(0), fem.locate_dofs_topological(V, fdim, boundary_facets), V)

# ------------------------------------------------------------------------------------------------

# Define solution variable and interpolate initial solution

uh = fem.Function(V)
uh.name = "uh"
uh.interpolate(initial_condition)

# ------------------------------------------------------------------------------------------------

# VARIATIONAL FORM

# define functions to eliminate/produce reaction terms
# Define the functions E and E_inv
def E(x):
return 1.0 * ((x[0]-d)**2 + x[1]**2 + x[2]**2 <= 0.250)

def E_inv(x):
return 1.0 * ((x[0]-d)**2 + x[1]**2 + x[2]**2 > 0.250)

x = ufl.SpatialCoordinate(domain)
v = ufl.TestFunction(V)
f = - 1.0 * kp * uh / Km / (1 + uh / Km) - 1.0 * kp * rho * uh / Km / (1 + uh / Km)
F = uh * v * ufl.dx + dt*ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx - (u_n + dt * f) * v * ufl.dx

problem = fem.petsc.NonlinearProblem(F, uh, bcs=[bc])

# ------------------------------------------------------------------------------------------------

# SOLVE

solver = nls.petsc.NewtonSolver(MPI.COMM_WORLD, problem)

for i in range(num_steps):
t += dt
print(t)

# Solve linear problem
solver.solve(uh)
uh.x.scatter_forward()

# Update solution at previous time step (u_n)
u_n.x.array[:] = uh.x.array

#nbar = fem.assemble_scalar(fem.form(1/(4/3 * np.pi * 0.5**3) * u_h.x.array[:] * ufl.dx))

return uh

print(rxn_dif(5,1))
``````

By adding:

``````from dolfinx import log
log.set_log_level(log.LogLevel.INFO)
``````

you can see that the Newton solver diverges. Make sure that you have a sensible initial guess to the non-linear equation, and follow the suggestions in Default absolute tolerance and relative tolerance - #4 by nate