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.

Thank you @dokken for the suggestions. I am new to DG method, I read some tutorials and posts for Poisson-dg in fenics. I implemented some very simple changes as follow (adding F_dg) and do find the method can ensure non-negative!

# Define normal vector and mesh size
n = FacetNormal(mesh)
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2

# Define parameters
alpha = 4.0
gamma = 8.0

# 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) \
      + k_on*u_1*u_2*v_1*dx(2) - k_off*u_3*v_1*dx(2)  \
      + ((u_2 - u_n2) / dt)*v_2*dx(2) + D_B*dot(grad(u_2), grad(v_2))*dx(2) \
      + k_on*u_1*u_2*v_2*dx(2) - k_off*u_3*v_2*dx(2)   \
      + ((u_3 - u_n3) / dt)*v_3*dx(2) + D_AB*dot(grad(u_3), grad(v_3))*dx(2) \
      - k_on*u_1*u_2*v_3*dx(2) + k_off*u_3*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)

F_dg = D_A*(- dot(avg(grad(v_1)), jump(u_1, n))*dS \
      - dot(jump(v_1, n), avg(grad(u_1)))*dS \
      + alpha/h_avg*dot(jump(v_1, n), jump(u_1, n))*dS \
      - dot(grad(v_1), u_1*n)*ds(1) \
      - dot(v_1*n, grad(u_1))*ds(1) \
      + (gamma/h)*v_1*u_1*ds(1)) 

F = F_ring + F_circle + F_dg

From my current understanding, the only difference of the DG variational form for reaction-diffusion systems in this case with the Poisson DG example in Fenics is to times the form by diffusion constant D_A of the DG terms in Poisson equation, since in my problem D_B and D_AB is zero so the DG variational form don’t have D_B and D_AB related diffusion terms. It seems now the solution is ensured to be non-negative but I am worried about the formulation of DG is not correct since it doesn’t consider information about subdomain ring and circle. I don’t have boundary conditions for the outer ring so I guess the boundary terms in DG can be set as zero, however, should I consider some DG terms related with the interface of the circle and ring? It seems that DG form is the same for the ring and circle subdomain so now I only have one DG term for whole ring and circle, does this look correct? Thank you!

Another question which might be related with DG is that I am currently using the backward Euler method, which should be stable for large k_{on} values, for example k_{on} = 10^9. However, the program can not converge with large k_{on} values which seems to be related to time discretization. Could you provide some insight into why this might be happening and how to address it?

What boundary conditions do you expect at this domain?
You have only written down your physics in words, making it hard for anyone to parse the meaning here.
Recall that dS integrals integrates over all interior facets, including those at the internal boundary.

You always have some boundary conditions in some sense, They might me Neumann, Dirichlet, Robin or some non-linear reaction. Please write down the strong formulation of your problem.

Your code is not complete and is diverging from the original question. Please make a new post