How to impose integral constraint on PDE?

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.

\frac{\partial a}{\partial t} - D\Delta a = f(a, b)\\ b = \frac{M}{S} - \frac{1}{S}\int_{\Omega}\,a\,\mathrm{d}x

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.

\int_{\Omega} (a^{k+1} - a^k)\,u\,\mathrm{d}x + \Delta t D\int_{\Omega}\,\nabla a^{k+1}\cdot \nabla u\,\mathrm{d}x - \int_{\Omega}\,f(a^{k+1}, \frac{M}{S} - m^{k+1})\,u\,\mathrm{d}x\\ + \int_{\Omega}\,(a^{k+1} - m^{k+1})\,r\,\mathrm{d}x = 0

The idea is to set

m^{k+1} = \frac{1}{S}\int_{\Omega}\,a^{k+1}\,\mathrm{d}x\quad\mathrm{and}\\ b^{k+1} = \frac{M}{S} - m^{k+1}

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:

  1. Is the above weak formulation wrong? If yes, can you suggest what would be the correct one?
  2. 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?
  3. 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_
    if useXmf:
        xmfout.write(C0.split()[0], t)#-----> Can't open the output files in Paraview 5.6.0
        pvdout << (C0.split()[0], t)
1 Like

I found the answer to my first two questions. For the first question, I needed to change the line below

a0, m0 = C0.split()


a0, m0 = split(C0)

For the second question, I found that the Lagrange multiplier is the last element in the vector corresponding to C0 and I can access it using C0.vector()[:][-1].
If someone can help me sort out the issue with XDMFFile it will be very helpful.

The extension of your XDMFFile is ‘.xmf’. This is probably not recognized by paraview. Try changing this to .xdmf, then it will probably work.


Thank you! That works.