How to impose integral constraint on PDE?

variational formulation
1

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_
    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)
1 Like
2

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

a0, m0 = C0.split()

to

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.

3

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

2 Likes
4

Thank you! That works.