Hello all, I’m attempting to simulate reaction-diffusion A+B\leftrightarrow AB within two circular domains centered at the same center: the first chemical A in the ring with initial concentration 0.05 with inner radius0 and outer radius radius=40, it diffuses into the inner circle with the second chemical B with initial concentration 100 with radius0=5, where the reaction occurs to form the third chemical AB reversibly. However, I’ve found the interface of AB concentration profile always appear to be not smooth near the interface, as the following figures shows, there are some positive values and negative values at the boundary of inner circle. When the wave of AB propagates inside the circle, the wavefront becomes smooth and symmetric however near the interface it is not smooth. Ideally the concentration field should be radially symmetric and positive near the interface as well. I thought this problem might be due to mesh or subdomain assignment. But even I did increase the mesh to 64 we can still see the non-smoothness and negative of the field AB or field A or field B at the ring-circle interface, as following figure shows (field AB). Any ideas of how to resolve this would be greatly appreciated. The file used to generate the mesh and subdomains are included as follows.
[Update1] The problem might not be meshing as I refine the inner mesh to radius of15 still the non-smooth happens in the boundary. It might be due to the discontinuity in the initial condition.
from fenics import *
from dolfin import *
from mshr import *
from matplotlib import pyplot as plt
import numpy as np
from mpi4py import MPI
import ufl
# Define the domain: Rectangle minus Circle
center = Point(0,0)
radius0 = 5
radius = 40
# Generate mesh and define spatial coordinates
domain = Circle(center, radius)
mesh = generate_mesh(domain, 64)
# Define refinement criteria
cell_markers = MeshFunction("bool", mesh, mesh.topology().dim())
for cell in cells(mesh):
p = cell.midpoint() # Get the midpoint of the cell
if (p[0] - center.x())**2 + (p[1] - center.y())**2 < (radius0 + 5)**2:
cell_markers[cell] = True # Mark for refinement
# Refine mesh
mesh = refine(mesh, cell_markers)
plot(mesh)
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)
# You can print subdomain markers if needed
print(subdomains.array().shape)
# Counting occurrences of 1 and 2
count_1 = np.count_nonzero(subdomains.array() == 1)
count_2 = np.count_nonzero(subdomains.array() == 2)
# Printing the counts
print(f"Number of 1s: {count_1}")
print(f"Number of 2s: {count_2}")
# Define the initial conditions for each component
u_1_initial = Expression("(x[0])*(x[0]) + (x[1])*(x[1]) <= 25.0 + DOLFIN_EPS ? 0: 0.5", degree=2)
u_2_initial = Expression("(x[0])*(x[0]) + (x[1])*(x[1]) <= 25.0 + DOLFIN_EPS ? 100 : 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())
u_1_proj = interpolate(u_1_initial, V.sub(0).collapse())
u_2_proj = interpolate(u_2_initial, V.sub(1).collapse())
u_3_proj = interpolate(u_3_initial, V.sub(2).collapse())
'''
