How to enforce a discontinuous initial condition for a continuous Lagrange function space

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


'''
![Screenshot 2024-07-08 at 7.08.05 PM|608x500](upload://uoFZiui5vT4NLg4CGQZPrXV4QnH.png)

AB field shown here, the interface between the circle and ring is not smooth and have negative values.

Your problem is not reproducible as you haven’t defined V. However, I would guess that V is a continuous Lagrange space, while

are discontinuous expressions, resulting in Gibb’s phenomenon, see for instance: Projection and interpolation — FEniCS Tutorial @ Sorbonne and references therein.

It would also be beneficial if you simplified the title of the topic, as it is incredibly problem specific at the moment, making it hard for others to relate to. I would probably rephrase it as:
How to enforce a discontinuous initial condition for a continuous Lagrange function space.

Thank you for the suggestions @dokken ! Yes, I am using the continuous Lagrange space. It seems the link Projection and interpolation — FEniCS Tutorial @ Sorbonne tells us even the mesh size becomes smaller, the Gibbs phenomenon still exist for the Lagrange space but when mesh size becomes smaller, the Gibbs phenomenon can go away with first order discontinuous lagrange DG1. Hence we should consider use DG1 as function space and reformulate the problem as DG variational form. Are there any simple ways to change or should I try to formulate the Discrete Galerkin variational form?

[Update1]: might not related with Gibbs, I drew the mesh using Gmsh [not updated in this post yet] and so that the there are more transition from ring to circle. The solution becomes more smooth at the boundary. However there are still negative values. So I guess the Gibbs phenomenon still persists even I used finer mesh!

Please see the original full code as follows.

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, 32)

# 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 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])*(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())

# 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 expressions used in variational forms
T = 1800 # final time, s
num_steps = 1800  # number of time steps
dt = T / num_steps  # time step size
D_A = 50  # diffusion coefficient
D_B = 0
D_AB = 0 
k_on = 10**2    # reaction on rate
k_off =  20         # reaction off rate 

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)

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_2))*dx(2) \
    + ((u_3 - u_n3) / dt)*v_3*dx(2) + D_AB*dot(grad(u_3), grad(v_3))*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_2))*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_3))*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

surface_integral_u1 = []
surface_integral_u2 = []
surface_integral_u3 = []

for n in range(num_steps):
    # Update current time
    t += dt
 
    # Solve variational problem for time step
    solve(F == 0, u)

    # Save solution to file (VTK)
    _u_1, _u_2, _u_3 = u.split()

    surface_integral_u1.append(assemble(_u_1 * dx))
    surface_integral_u2.append(assemble(_u_2 * dx))
    surface_integral_u3.append(assemble(_u_3 * dx))
    print(assemble(_u_1 * dx), assemble(_u_2 * dx), assemble(_u_3 * dx))

    if n%100 == 0: 
        vtkfile_u_1 << (_u_1, t)
        vtkfile_u_2 << (_u_2, t)
        vtkfile_u_3 << (_u_3, t)
   
    # Update previous solution
    u_n.assign(u)
    
    set_log_level(LogLevel.PROGRESS)
    progress += 1
    set_log_level(LogLevel.ERROR)

It seems that I can not make changes to the title and the old post now after a few days. Let me know if I miss any places to change the title. Thank you!

You would have to reformulate your variational problem using DG/interior penalty/Nitsche terms to use DG for your test and trial functions.
An alternative way for resolving this is letting u_n live in a DG space, and then at each time step interpolate the new solution u into u_n. As the DG space is bigger than the continuous space, you would retain an exact solution at future time steps.