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)