Computing the Jacobian manually for a coupled problem (e.g. Cahn-Hilliard equation demo)

Question

The forms for the Cahn-Hilliard demo (2d) are generated from the UFL file CahnHilliard2D.ufl:

cell = triangle

P1 = FiniteElement("Lagrange", cell, 1)
ME = P1*P1

du   = TrialFunction(ME)
q, v = TestFunctions(ME)

u   = Coefficient(ME)  # current solution
u0  = Coefficient(ME)  # solution from previous converged step

# Split mixed functions
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

lmbda = Constant(cell) # surface energy parameter
dt = Constant(cell) # time step
theta = Constant(cell) # time stepping parameter

# mu_(n+theta)
mu_mid = (1-theta)*mu0 + theta*mu

# Compute the chemical potential df/dc
c = variable(c)
f = 100*c**2*(1-c)**2
dfdc = diff(f, c)

F0 = c*q*dx  - c0*q*dx   + dt*dot(grad(mu_mid), grad(q))*dx
F1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
F  = F0 + F1

J = derivative(F, u, du)

Here, the form of the Jacobian J is calculated automatically by taking the Gateaux derivative through the function derivative on the last line. I want to compute the Jacobian manually.

My motivation is that I am trying to implement a different coupled problem, where I am using "Quadrature" elements to incorporate

• two different stress-type functions
  
  
  • $\boldsymbol{\sigma}(\boldsymbol{\varepsilon}, d)$ and
  • $\boldsymbol{\beta}(\boldsymbol{\varepsilon}, d)$
    which are functions of the two coupled fields $\boldsymbol{\varepsilon} = \frac{1}{2}(\nabla\mathbf{u} + \nabla\mathbf{u}^T)$ and $d$, and which are substituted into two functionals, say $F_0$ and $F_1$, respectively.
• I also have four tangent moduli that are related to these stresses,
  
  
  • $\frac{\partial \boldsymbol{\sigma}}{\partial \boldsymbol{\varepsilon}}$,
  • $\frac{\partial \boldsymbol{\sigma}}{\partial d}$,
  • $\frac{\partial \boldsymbol{\beta}}{\partial \boldsymbol{\varepsilon}}$, and
  • $\frac{\partial \boldsymbol{\beta}}{\partial d}$
    These show up in the Jacobian terms
  • $\frac{\partial F_0}{\partial \mathbf{u}}$,
  • $\frac{\partial F_0}{\partial d}$,
  • $\frac{\partial F_1}{\partial \mathbf{u}}$, and
  • $\frac{\partial F_1}{\partial d}$.

Now, it wouldn't make sense for me to present my whole problem here and expect you to solve it. But I believe that I can work out the solution myself if I could see this done in one of the FEniCS demos, namely the Cahn-Hilliard equation. Is it possible to reformulate the demo so that one can manually assemble the manually computed Jacobian terms, i.e. $\frac{\partial F_0}{\partial c}$, $\frac{\partial F_0}{\partial \mu}$, $\frac{\partial F_1}{\partial c}$, and $\frac{\partial F_1}{\partial \mu}$?

I have the following start

d2fdc2 = diff(dfdc, c)
dF0dc = q * dx
dF0dmu = dt * theta * div(q) * dx
dF1dc = - d2fdc2 * v * dx - lmbda * div(v) * dx
dF1dmu = v * dx

But I don't know how to proceed. Do you have any ideas?

Answer 1

As in the example given by MiroK, it is sufficient to just sum up the directional derivatives

d2fdc2 = diff(dfdc, c)

dL0dc = dc * q * dx
dL0dmu = dt * theta * dot(grad(dmu), grad(q)) * dx
dL1dc = - dc * d2fdc2 * v * dx -  lmbda * dot(grad(dc), grad(v)) * dx
dL1dmu = dmu * v * dx

a = dL0dc + dL0dmu + dL1dc + dL1dmu

This yields the same result as

a = derivative(L, u, du)

Answer 2

Hi, perhaps the following helps with your question

import random
from dolfin import *

# Class representing the intial conditions
class InitialConditions(Expression):
    def __init__(self, **kwargs):
        random.seed(2 + MPI.rank(mpi_comm_world()))
    def eval(self, values, x):
        values[0] = 0.63 + 0.02*(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):
        NonlinearProblem.__init__(self)
        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)

# Model parameters
lmbda  = 1.0e-02  # surface parameter
dt     = 5.0e-06  # time step
theta  = 0.5      # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson

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

# Create mesh and build function space
mesh = UnitSquareMesh(96, 96)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
ME = FunctionSpace(mesh, P1*P1)

# 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
dc, dmu = split(du)
c,  mu  = split(u)
c0, mu0 = split(u0)

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

# Compute the chemical potential df/dc
c = variable(c)
f    = 100*c**2*(1-c)**2
dfdc = diff(f, c)

# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu

# Weak statement of the equations
L0 = c*q*dx - c0*q*dx + dt*dot(grad(mu_mid), grad(q))*dx
L1 = mu*v*dx - dfdc*v*dx - lmbda*dot(grad(c), grad(v))*dx
L = L0 + L1

# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
A = assemble(a)
print '|A|', A.norm('linf')

x, y = split(u)
b = derivative(L0, x, dc) + derivative(L0, y, dmu) + \
    derivative(L1, x, dc) + derivative(L1, y, dmu)
B = assemble(b)
print '|B|', B.norm('linf')

# See that they are the same
A.axpy(-1, B, True)
print '|A-B|', A.norm('linf') 

Comments

So I can define the bilinear form like this:

a = dL0dc + dL0dmu + dL1dc + dL1dmu

where the rhs terms are as I have written in the question body. However, the sum of those values do not result in a bilinear form and I am getting an error, because e.g. dL0dc = q * dx is obviously not a bilinear form.

Should it be though? I don't understand why
$$\frac{\partial L_0}{\partial c} = \int_{\Omega} q\,\mathrm{d}x$$
would not be true? Am I not calculating correctly?

---

Yes, the result should be a bilinear form. I'd need to know what is $L_0$ to answer the question about its derivative.

---

I just realized I am making conceptual mistakes and it doesn't have to do with FEniCS. I am currently working on it and will let you know if I come up with something.