I’m new to Fenics, it looks really poweful! I’m attempting to simulate reaction-diffusion within two circular domains centered at the same center: the first chemical A in the ring with initial concentration 0.05, it diffuses into the inner circle with the second chemical B with initial concentration 1, where the reaction occurs to form the third chemical AB reversibly. However, I’ve encountered an issue where the program runs for a few iterations but the solver fails to converge for longer iterations and stopped. I tried both including or not including boundary conditions bcs, both would have this bug. Also I tried relaxing the convergence tol and including relaxation, still not working. I tried to decrease the timestep to be 0.1 s and 0.01s both will stop after some runs at similar simulation time. Maybe the problem lies in mesh or variational formulation of the problem? I tried to reduce k_on to 10 and the program can converge for all iterations which seems weird. Any insights into resolving the problem or improving the Fenics code would be greatly appreciated.

```
from fenics import *
from dolfin import *
from mshr import *
from matplotlib import pyplot as plt
import numpy as np
from mpi4py import MPI
# Initialize MPI
comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()
# Define the domain: Rectangle minus Circle
center = Point(25, 25)
radius0 = 5
radius = 25
# Generate mesh and define spatial coordinates
domain = Circle(center, radius)
mesh = generate_mesh(domain, 64)
x = SpatialCoordinate(mesh)
# Define subdomains
class CircleSubdomain(SubDomain):
def inside(self, x, on_boundary):
return (x[0] - center.x())**2 + (x[1] - center.y())**2 < radius0**2
# Initialize subdomain marker
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim())
subdomains.set_all(2) # Initialize all subdomains to 2
# Mark subdomain inside the circle as 1
circle_subdomain = CircleSubdomain()
circle_subdomain.mark(subdomains, 1)
# Define function space for system of concentrations and basis functions on this space
P1 = FiniteElement('P', triangle, 1)
element = MixedElement([P1, P1, P1])
V = FunctionSpace(mesh, element)
# Define test functions
v_1, v_2, v_3 = TestFunctions(V)
u = Function(V)
u_n = Function(V)
# Split system functions to access components
u_1, u_2, u_3 = split(u)
u_n1, u_n2, u_n3 = split(u_n)
# Define the initial conditions for each component
u_1_initial = Expression("(x[0]-25)*(x[0]-25) + (x[1]-25)*(x[1]-25) <= 25.0 + DOLFIN_EPS ? 0: 0.05", degree=2)
u_2_initial = Expression("(x[0]-25)*(x[0]-25) + (x[1]-25)*(x[1]-25) <= 25.0 + DOLFIN_EPS ? 1 : 0", degree=2)
u_3_initial = Constant(0)
# Project initial conditions
u_1_proj = project(u_1_initial, V.sub(0).collapse())
u_2_proj = project(u_2_initial, V.sub(1).collapse())
u_3_proj = project(u_3_initial, V.sub(2).collapse())
# Assign initial conditions to the function
assign(u_n.sub(0), u_1_proj)
assign(u_n.sub(1), u_2_proj)
assign(u_n.sub(2), u_3_proj)
# Define Dirichlet boundary conditions (used or not)
u_D1 = Constant(0.05) # Example Dirichlet condition for component 1
u_D2 = Constant(0.0) # Example Dirichlet condition for component 2
u_D3 = Constant(0.0) # Example Dirichlet condition for component 3
# Define Dirichlet boundary conditions function with DOLFIN_EPS
def boundary(x, on_boundary):
tol = DOLFIN_EPS
return on_boundary and abs(sqrt((x[0] - center.x())**2 + (x[1] - center.y())**2) - radius) < tol
# Define boundary conditions
bc1 = DirichletBC(V.sub(0), u_D1, boundary)
bc2 = DirichletBC(V.sub(1), u_D2, boundary)
bc3 = DirichletBC(V.sub(2), u_D3, boundary)
# Collect boundary conditions
bcs = [bc1, bc2, bc3]
# Define expressions used in variational forms
T = 1800 # final time
num_steps = 18000 # number of time steps
dt = T / num_steps # time step size
D_A = 0.1 # diffusion coefficient
D_B = 0
D_AB = 0
k_on = 10**5
k_off = 10**(-6)
D_A = Constant(D_A)
D_B = Constant(D_B)
D_AB = Constant(D_AB)
k_on = Constant(k_on)
k_off = Constant(k_off)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
# Variational form with subdomains
F_ring = ((u_1 - u_n1) / dt)*v_1*dx(2) + D_A*dot(grad(u_1), grad(v_1))*dx(2) \
+ ((u_2 - u_n2) / dt)*v_2*dx(2) + D_B*dot(grad(u_2), grad(v_1))*dx(2) \
+ ((u_3 - u_n3) / dt)*v_3*dx(2) + D_AB*dot(grad(u_3), grad(v_1))*dx(2)
F_circle = ((u_1 - u_n1) / dt)*v_1*dx(1) + D_A*dot(grad(u_1), grad(v_1))*dx(1) \
+ k_on*u_1*u_2*v_1*dx(1) - k_off*u_3*v_1*dx(1) \
+ ((u_2 - u_n2) / dt)*v_2*dx(1) + D_B*dot(grad(u_2), grad(v_1))*dx(1) \
+ k_on*u_1*u_2*v_2*dx(1) - k_off*u_3*v_2*dx(1) \
+ ((u_3 - u_n3) / dt)*v_3*dx(1) + D_AB*dot(grad(u_3), grad(v_1))*dx(1) \
- k_on*u_1*u_2*v_3*dx(1) + k_off*u_3*v_3*dx(1)
# Combine the variational forms
F = F_ring + F_circle
# Create VTK files for visualization output
vtkfile_u_1 = File('u_1.pvd')
vtkfile_u_2 = File('u_2.pvd')
vtkfile_u_3 = File('u_3.pvd')
progress = Progress('Time-stepping', num_steps)
# Time-stepping
t = 0
counter = 0
for n in range(num_steps):
# Update current time
t += dt
solve(F == 0, u)
_u_1, _u_2, _u_3 = u.split()
vtkfile_u_1 << (_u_1, t)
vtkfile_u_2 << (_u_2, t)
vtkfile_u_3 << (_u_3, t)
u_n.assign(u)
set_log_level(LogLevel.PROGRESS)
progress += 1
set_log_level(LogLevel.ERROR)
```