Hi all.

I am trying to solve a reaction-diffusion equation with two species A and B with mass conservation. The diffusion coefficient of B is much much larger than that of A so we assume that concentration of B is always uniform throughout but varies with time depending on mass of A. The reaction term f is non-linear. The governing PDE and the integral constraint are as follows. S is the surface area, M is the conserved total mass of A and B and D is the diffusion coefficient of A.

Referring to the answer given by @ChristianWaluga to this question, I came up with a weak form as follows. I create a Mixed Function Space with a P1 element to represent the concentration of A and a global real number element to represent its mean. In the equation below a^{k+1}, u and m^{k+1}, r are the trial and test functions from P1 and R space respectively.

The idea is to set

Assume no-flux boundary condition on all boundaries. I wrote the following code but it does not do anything. The output is the same as initial condition at all time steps. So I have following questions:

- Is the above weak formulation wrong? If yes, can you suggest what would be the correct one?
- I expected
`C.split()[1]`

to give me a single real number but it gives me an array of same size as`C.split()[0]`

. Am I not setting the global integral constraint correctly? - In the code below, when I print output using
`XDMFFile`

, I cannot open the output in Paraview 5.6.0, it gives error`Failed to read data`

. Is there any mistake in my code?

I am using Fenics 2018.1.0 from conda-forge on Arch Linux.

```
from dolfin import *
mesh = UnitSquareMesh(100, 100)
P1 = FiniteElement('Lagrange', mesh.ufl_cell(), 1)
R = FiniteElement('Real', mesh.ufl_cell(), 0)
V = FunctionSpace(mesh, P1*R)
# The initial conditions
C0 = Function(V)
a0, m0 = C0.split()
# The non-linear reaction term
def f(a, b):
return a*(1.0 - a)*(a - 1.0 - b)
A = Constant(2.0) # Conserved total density
D = Constant(0.05) # Diffusion coefficient
dt = Constant(0.1) # Time step
# Setting up the non linear variational problem
a_, m_ = TrialFunctions(V)
u, r = TestFunctions(V)
b = A - m_
C = Function(V)
F = ((a_ - a0)*u + dt*D*dot(grad(a_), grad(u)) - dt*u*f(a_, b) + (a_ - m_)*r)*dx
F = action(F, C)
J = derivative(F, C)
problem = NonlinearVariationalProblem(F, C, J=J)
solver = NonlinearVariationalSolver(problem)
# Set initial condition and evolve in time
cpp_code = '((x[0] - 0.2)*(x[0] - 0.4) + (x[1] - 0.3)*(x[1] - 0.3)) > 0.04? 0.0 : 2.6651'
ic = Expression(cpp_code, degree=2)
int_a0 = interpolate(ic, V.sub(0).collapse())
int_m0 = interpolate(Expression('0.335', degree=0), V.sub(1).collapse())
assign(C0, [int_a0, int_m0])
T = 1.0
t = 0.0
dt_ = 0.01
xmfout = XDMFFile('Out.xmf')
pvdout = File('Out.pvd')
useXmf = False #Switch for XMF or VTU output file format
while t < T:
t += dt_
solver.solve()
C0.assign(C)
if useXmf:
xmfout.write(C0.split()[0], t)#-----> Can't open the output files in Paraview 5.6.0
else:
pvdout << (C0.split()[0], t)
```