Hi,
I am very new to Fenics and currently learning to use Fenics in my work. As a starting point, I am trying to replicate some literature. The problem is coupled PDE associated with simple chemical reactions.
I noticed that there are differences in solutions from the literature. It seems like the Robin boundary condition at x = 0 causing differences. I checked with other software too and was able to replicate the literature data. Can anybody know what seems to be a problem? Thanks in advance.
from fenics import *
import matplotlib.pyplot as plt
mesh = UnitIntervalMesh(10)
P1 = FiniteElement(‘P’, mesh.ufl_cell(), 2)
element = MixedElement([P1, P1])
V = FunctionSpace(mesh, element)
v_1, v_2 = TestFunctions(V)
u = Function(V)
u_1, u_2 = split(u)
color = MeshFunction(“size_t”, mesh, 0)
color.set_all(0)
bc1 = ‘x[0] == 0’
CompiledSubDomain(bc1, tol=1e-15).mark(color, 1)
ds = Measure(‘ds’, domain=mesh, subdomain_data=color)
Pe, Bo = Constant(10), Constant(10)
Hw = Constant(-0.05)
def R(u1, u2):
return 4.0u1exp(18*(1-u2**-1))
#end def
u.interpolate(Constant((0.1, 1.0)))
w = as_vector([1.0])
F1 =1/Pedot(grad(u_1), grad(v_1))dx - u_1v_1ds(1) + v_1ds(1)
+ dot(w, grad(u_1))v_1dx
+ R(u_1, u_2)v_1dx
F2 =1/Bodot(grad(u_2), grad(v_2))dx - u_2v_2ds(1) + v_2ds(1)
+ dot(w, grad(u_2))v_2dx
+ Hw*R(u_1, u_2)v_2dx
F = F1 + F2
solve(F == 0, u)