Thanks for getting back to me. I’m quite new to the finite element method so I’m not too sure how to implement stabilization for my problem. Currently, the problem I’m trying to solve is as follows:

Taking C_f to be u1 and B to be u2 I implement the following :

```
#Define function space for system of concentrations
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
element = MixedElement([P1,P1])
V = FunctionSpace(mesh, element)
#Define test functions
v1,v2 = TestFunction(V)
#Define functions for concentration
u = Function(V)
u_n = Function(V)
u1,u2 = split(u)
u_n1,u_n2 = split(u_n)
#Define terms
Da = 50000
D_T = 2.041*pow(10,-4)
T = 0.5
B_M = 0.356*pow(10,-6)
K_d = 0.0326*pow(10,-6)
ka = (D_T*Da)/(B_M*(T*T))
kd = ka*K_d
k1 = 0.009221
A1 = 0.24
Z_w = 966.21*pow(10,3)
tol = 1E-9
sim_time = 60.0
time_steps = 600
dt = sim_time/time_steps
#Define Initial Conditions
u_n1 = interpolate(Constant(0),V.sub(0).collapse())
u_n2 = interpolate(Constant(0),V.sub(1).collapse())
#Define Boundary Conditions
Flux = Expression("t<30+1e-9 ? ((k1*A1)/Z_w)*exp(-1*k1*t) : 0", degree=2, k1=k1, A1=A1, Z_w=Z_w, t=0)
u1_D = Constant(0)
boundaryconditions = {
1 : {"Neumann" : Flux},
2 : {"Dirichlet" : u1_D}
}
ds = Measure('ds', domain=mesh, subdomain_data=mf)
bcs = []
for i in boundaryconditions:
if "Dirichlet" in boundaryconditions[i]:
# bc = DirichletBC(V, boundaryconditions[i]["Dirichlet"], boundaries, i)
bc = DirichletBC(V.sub(0), boundaryconditions[i]["Dirichlet"], boundaries, i) # Dirichlet is applied only to 1 subspace according to the paper
bcs.append(bc)
integrals_N = []
for i in boundaryconditions:
if "Neumann" in boundaryconditions[i]:
if boundaryconditions[i]["Neumann"] != 0: #if not zero flux
g = boundaryconditions[i]["Neumann"]
integrals_N.append(v1*g*ds(i))
#Define Variational Problem
F = (((u1-u_n1 + u2-u_n2)*v1)/dt)*dx + D_T*dot(grad(u1),grad(v1))*dx + (((u2-u_n2)*v2)/dt)*dx - sum(integrals_N) \
- v2*ka*u1*B_M*dx + ka*u1*v2*u2*dx + kd*v2*u2*dx
#Solve problem
t = -1*dt
J = derivative(F, u)
for n in range(time_steps+1):
t += dt
Flux.t = t
# solve(F == 0, u, bcs, J=J)
solve(F == 0, u, bcs, J=J,
solver_parameters={"newton_solver":{"absolute_tolerance":1e-30,
"linear_solver": "gmres"}})
u1,u2 = u.split()
vtkfile_u1 << (u1,t)
vtkfile_u2 << (u2,t)
u_n.assign(u)
```

The example I’m familiar with with regards to SUPG as seen here deals with velocity. I’m not sure how to implement stabilization to my problem. I suppose its not accurate to refer to this as a advection problem due to the lack of a velocity field; however, I do intend on introducing one down the line.

Are there any examples I can study to stabilize my current problem? I’m not too familiar with stabilization methods and how to employ them. Any advice would be greatly appreciated and thanks for taking the time.

Also, I’m **not too sure if I’ve implemented the flux (time dependent) correctly**. The magnitude of the concentration in the mesh drops significantly after 30s when the flux changes to zero. Which should not be the case.