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

1

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))
2

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