Solution to Cahn-Hilliard equation does not conserve mass

9 days ago by

I am trying to solve the Cahn-Hilliard equation using some code based on the Fenics Cahn-Hilliard demo. The code is able to find a physically reasonable solution, but the total mass of the system is not conserved. I did the same check on the original demo and the concentration is conserved. My equation is slightly different than that of the demo:

$\frac{\partial h}{\partial t}-\nabla\cdot\left(\frac{h^3}{3\eta}\nabla\left(-\sigma\nabla^2h-\Pi\left(h\right)\right)\right)=0$ht ·(h33η (σ2hΠ(h)))=0

I have tried both the default Neumann and periodic boundary conditions, but it doesn't seem to matter. Different size of the time step also does not seem to make a difference. The mass lost appears to grow proportionally with the velocity field (the part to the right of the divergence operator). As time goes on, the fluctuations in the system get larger and so does the velocity and the lost mass. I just cannot seem to find where the mass is going. Any suggestions?

Here's the complete code:

import random
from datetime import datetime
from dolfin import *

Lx = 10000
Ly = 10000
Nx = 20
Ny = 20
T  = 1000
dt = 5.0e-6  # time step
eta = 18E-6*(1E-9) # in Pascal.seconds (normalized to nm)
sigma = 0.01 # in N/m
PI0 = 1.0055E-21*(1E18) # in Joules (normalized to nm)
hc = 3.6 # in nanometers
L = 1.33 # in nanometers

# Class representing the intial conditions
class InitialConditions(UserExpression):
    def __init__(self, **kwargs):
        #random.seed(2 + MPI.rank(MPI.comm_world))
    def eval(self, values, x):
        values[0] = 6 - 0.01*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)

# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
    def __init__(self, a, L):
        self.L = L
        self.a = a
    def F(self, b, x):
        assemble(self.L, tensor=b)
    def J(self, A, x):
        assemble(self.a, tensor=A)

# Sub domain for Periodic boundary condition
class PeriodicBoundary(SubDomain):

    # Left and bottom boundaries are "target domain" G
    def inside(self, x, on_boundary):
        # return True if on left or bottom boundaries
        return bool((near(x[0], 0) or near(x[1],0))  and on_boundary)

    # Map top and right boundaries (H) to bottom and left boundaries (G)
    def map(self, x, y):
        # map top right corner
        if near(x[0], Lx) and near(x[1], Ly):
            y[0] = x[0] - Lx
            y[1] = x[1] - Ly
        # map right edge
        elif near(x[0], Lx):
            y[0] = x[0] - Lx
            y[1] = x[1]
        # map top edge
            y[0] = x[0]
            y[1] = x[1] - Ly

# Model parameters

theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
c = 5.1301993206474563
beta = 0.37

# Form compiler options
parameters["form_compiler"]["optimize"]     = True
parameters["form_compiler"]["cpp_optimize"] = True

# Create mesh and build function space
#mesh = UnitSquareMesh.create(96, 96, CellType.Type.quadrilateral)
mesh = RectangleMesh.create( MPI.comm_world, [Point(0,0), Point(Lx,Ly)], [Nx,Ny], CellType.Type.quadrilateral)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)
#ME = FunctionSpace(mesh, P1*P1, constrained_domain=PeriodicBoundary())

# Trial and test functions of the space ``ME`` are now defined::
# Define trial and test functions
du    = TrialFunction(ME)
q, v  = TestFunctions(ME)

# Define functions
u   = Function(ME)  # current solution
u0  = Function(ME)  # solution from previous converged step

# Split mixed functions
dh, dmu = split(du)
h,  mu  = split(u)
h0, mu0 = split(u0)

# Create intial conditions and interpolate
u_init = InitialConditions(degree=1)

n = FacetNormal(mesh)

# which is then used in the definition of the variational forms::
# Weak statement of the equations
R0 = (1/(3*eta))*h0**3*dot(grad(mu0),grad(q)) - (1/(3*eta))*dot(grad(mu0),grad(h0**3))*q
R = (1/(3*eta))*h**3*dot(grad(mu),grad(q)) - (1/(3*eta))*dot(grad(mu),grad(h**3))*q

L0 = h*q*dx - h0*q*dx + theta*R*dt*dx + (1-theta)*R0*dt*dx

PI = PI0*(hc-h)/(h**3*(h+L))
L1 = mu*v*dx + PI*v*dx - sigma*dot(grad(h), grad(v))*dx

L = L0 + L1

# This is a statement of the time-discrete equations presented as part
# of the problem statement, using UFL syntax. The linear forms for the
# two equations can be summed into one form ``L``, and then the
# directional derivative of ``L`` can be computed to form the bilinear
# form which represents the Jacobian matrix::

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)

# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L)
solver = NewtonSolver()
solver.parameters["linear_solver"] = "lu"
solver.parameters["convergence_criterion"] = "incremental"
solver.parameters["relative_tolerance"] = 1e-6

# Output file
file = File("output/output.pvd", "compressed")
mass_file = open("mass.txt", "w")

# Step in time
t = 0.0

file << (u.split()[0], t)
mass0 = assemble(h*dx)
mass_file.write("#time\t\ttotal mass\t\tdm\t\tdm/m\n")
mass = mass0
mass_file.write("{0}\t\t{1}\t\t{2}\t\t{3}\n".format(t, mass, mass-mass0, (mass-mass0)/mass0))
for l in range(1, T):
    t += dt
    u0.vector()[:] = u.vector()
    solver.solve(problem, u.vector())
    mass = assemble(h*dx)
    mass_file.write("{0}\t\t{1}\t\t{2}\t\t{3}\n".format(t, mass, mass-mass0, (mass-mass0)/mass0))
    print("t= {0}, dm = {1}".format(l, (mass-mass0)/mass0))
    file << (u.split()[0], t)​
Community: FEniCS Project
Can you please write down explicitly what you mean by Neumann and periodic. You need the net flux (the quantity in the divergence dot the normal) to integrate to zero on the boundary.
written 9 days ago by K D  
I am solving the problem as two coupled equations:

  $\frac{\partial h}{\partial t}-\nabla\cdot\left(\frac{h^3}{\eta}\nabla\mu\right)=0$ht ·(h3η μ)=0  



As I understand, the boundary conditions for this problem are

$\left(\frac{h^3}{\eta}\nabla\mu\right)\cdot n=0$(h3η μ)·n=0  on  $\partial\Omega$∂Ω

$\nabla h\cdot n=0$h·n=0  on  $\partial\Omega$∂Ω

Assuming Neumann boundary conditions for the problem, meaning that I am using the default BC settings in fenics without specifying anything at the boundaries, the two conditions above should be met without having to add any additional constraints to the linear form:

$L=\int\frac{\partial h}{\partial t}qdx+\int\frac{h^3}{\eta}\left(\nabla\mu\cdot\nabla q\right)dx-\int\frac{1}{\eta}\left(\nabla\mu\cdot\nabla h^3\right)qdx+\int\mu vdx-\int\sigma\left(\nabla h\cdot\nabla v\right)dx+\int\Pi\left(h\right)vdx$L=ht qdx+h3η (μ·q)dx1η (μ·h3)qdx+μvdxσ(h·v)dx+Π(h)vdx

Even if I add an additional term to the linear form such as  $\int\left(\frac{h^3}{\eta}\nabla\mu\right)\cdot nqdx$(h3η μ)·nqdx , it doesn't make a difference in the mass conservation.

This linear form above follows the same pattern as in the Cahn-Hilliard demo on the web ( ), which I have tested and seems to conserve mass.

For periodic boundary conditions, I followed various tutorials/forum posts to map the the position of the top edge to the bottom, and likewise the right edge to the left.

Maybe I'm not understanding correctly the boundary conditions of the problem to make sure there is no net flux at the boundary?I think I am doing it the same way as in the demo, but maybe I'm missing something.
written 9 days ago by Juan Vanegas  
Hello Juan,

I don't understand this term:
-\int\frac{1}{\eta} \left(\nabla\mu\cdot\nabla h^3\right) q \ dx

Are you sure it should be there?  I don't think so.
written 9 days ago by K D  
Hello K D,

Starting with the equation at the top

$\frac{\partial h}{\partial t}-\nabla\cdot\left(\frac{h^3}{\eta}\nabla\mu\right)=0$ht ·(h3η μ)=0

and using the product rule for the divergence  $\nabla\cdot\left(f\left(h\right)A\right)=f\left(h\right)\nabla\cdot A+A\cdot\nabla f\left(h\right)$·(ƒ (h)A)=ƒ (h)·A+A·ƒ (h) I get

$\frac{\partial h}{\partial t}-\frac{h^3}{\eta}\nabla^2\mu-\nabla\mu\cdot\nabla\left(\frac{h^3}{\eta}\right)=0$ht h3η 2μμ·(h3η )=0

This is where this term is coming from.

Am I doing something wrong here?
written 9 days ago by Juan Vanegas  
Hi Juan,

I would do the following (ignoring some constants for clarity). The weak form of \( \nabla \cdot \left(h^3 \nabla \mu\right)  \) is given by
\int q \nabla \cdot \left( h^3 \nabla \mu \right)= \int h^3 \nabla q \cdot \nabla \mu
where I simply treat \( h^3 \nabla \mu \) as one term in the integration by parts.

In your approach, you should have an extra term with \(h^3 \nabla^2 \mu\) which is a real pain to deal with (though keeping it there may fix the mass conservation).

PS.  Are you the person at U. Vermont?  Your papers with Arroyo are on my reading list.
written 9 days ago by K D  
Hello K D,

I'm glad I made it to your reading list! Thank you very much for your suggestion. If I remove the term you suggest the mass is now conserved in the code. I followed your suggestion for the integration by parts and I get the same formula you wrote above, but I still don't know where I went wrong the first time. You indicate there should be an additional term  $h^3\nabla^2\mu$h32μ, but I just don't see where it would come from. Thanks again!
written 9 days ago by Juan Vanegas  
Hello Juan,

In your post further up where you wrote out the product rule, you have the term \( h^3 \nabla^2 \mu\) in your equation.  When you convert to the weak form, this term becomes 
\int q h^3 \nabla \cdot \nabla \mu
\int \nabla(h^3 q) \cdot \nabla \mu
which is not quite what you have in \( L \) in another post even further up.  Notice that you can't pull \( h \) out of the \( \nabla (\cdot) \).
written 9 days ago by K D  
Hello K D,

My mistake was that I was integrating by parts only the laplacian of mu:

$h^3\int q\nabla^2\mu dx=-h^3\int\nabla\mu\cdot\nabla qdx$h3q2μdx=h3μ·qdx

and not as you did above. Now I see that was clearly wrong.

For completeness, the correct linear form for this problem should be:

$L=\int\frac{\partial h}{\partial t}qdx+\int\frac{h^3}{\eta}\left(\nabla\mu\cdot\nabla q\right)dx+\int\mu vdx-\int\sigma\left(\nabla h\cdot\nabla v\right)dx+\int\Pi\left(h\right)vdx$L=ht qdx+h3η (μ·q)dx+μvdxσ(h·v)dx+Π(h)vdx

written 6 days ago by Juan Vanegas  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »