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.0*u1*exp(18*(1-u2**-1))

#end def

u.interpolate(Constant((0.1, 1.0)))

w = as_vector([1.0])

F1 =1/Pe*dot(grad(u_1), grad(v_1)) dx - u_1v_1*ds(1) + v_1

*ds(1)*

+ dot(w, grad(u_1))

+ R(u_1, u_2)

F2 =1/Bodot(grad(u_2), grad(v_2))

+ dot(w, grad(u_1))

*v_1*dx+ R(u_1, u_2)

*v_1*dxF2 =1/Bo

*dx - u_2*v_2

*ds(1) + v_2*ds(1)

+ dot(w, grad(u_2))

*v_2*dx

+ Hw*R(u_1, u_2)

*v_2*dx

F = F1 + F2

solve(F == 0, u)