I have a modified Cahn-Hilliard style model to implement, this follows on from a previous question about initial conditions. I have implemented the model, but it will not run for more than around 5 time-steps before failing to converge. Examining the solutions obtained, the chemical potential starts well-behaved before oscillating wildly near the left boundary.
It's entirely possible that there are very basic errors in the formulation or implementation of the problem, but I have been wrestling with this for a couple of weeks now and would appreciate any insights!
The model 1 (equation 19) can be summarised as
$$\mathrm{p} \frac{\partial \phi_{0}}{\partial t} - \mathrm{div}\left(\hat{C} \nabla \phi_{0}\right) = \lambda \mathrm{div}\left(\hat{D} \nabla f\left(\phi_{0}\right)\right) - \frac{\lambda}{p}\mathrm{div}\left(\hat{D}\nabla\left(\mathrm{div}\left(\hat{D} \nabla \phi_{0}\right)-\tilde{g}_{0}\right)\right)$$
This is split by defining an effective chemical potential, obtaining the second order system. We also assume that $\tilde{g}_{0}$ is zero for simplicity.
$$\mu = p f\left(\phi_{0}\right) -\mathrm{div}\left(\hat{D}\nabla\phi_0\right)$$
$$p \frac{\partial \phi_{0}}{\partial t} - \mathrm{div}\left(\hat{C}\nabla\phi_{0}\right) = \frac{\lambda}{p} \mathrm{div}\left(\hat{D}\,\mathrm{grad}\left(\mu\right)\right)$$
This is solved in a channel with no-flux on both $\phi_{0}$ and $\mu$ on the top and bottom walls. Inlet and outlet BCs specify both $\phi_{0}$ and $\mu$ in a consistent manner. I have initially neglected the tensors as I'm not clear on how to handle them in the form at the moment.
My current code is
def macro_solver(mesh, boundaries, Q, D, C, area):
# Class representing the intial conditions
class InitialConditions(Expression):
def eval(self, values, x):
# Initial conditions come from steady state tanh behaviour for this model,
# boundaries match BCs
values[0] = -1.1739*np.tanh(3*(x[0]-1-0.05*cos(6.28*x[1])))+0.1739
values[1] = 2 - x[0]/25
def value_shape(self):
return (2,)
# Class for interfacing with the Newton solver
class CahnHilliardEquation(NonlinearProblem):
def __init__(self, a, L, bcs):
NonlinearProblem.__init__(self)
self.L = L
self.a = a
self.bcs = bcs
def F(self, b, x):
assemble(self.L, tensor=b)
[bc.apply(b,x) for bc in self.bcs]
def J(self, A, x):
assemble(self.a, tensor=A)
[bc.apply(A) for bc in self.bcs]
# Model parameters
dt = 2.0e-06 # time step
theta = 0.5 # time stepping family, e.g. theta=1 -> backward Euler, theta=0.5 -> Crank-Nicolson
p = 1.0/area # Porosity
lmbda = 1.0 # Eelastic relaxation time
# 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
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 2)
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. Changes in this expression affect the boundary conditions below
c = variable(c)
f = 3.0*(c**2-1.0)**2
dfdc = diff(f, c)
# mu_(n+theta)
mu_mid = (1.0-theta)*mu0 + theta*mu
#NOTE HERE - q -> c, v -> mu
# Weak statement of the equations - implied no flux for both mu and c on walls. Overwritten
# by Dirichlet BCs later
L0 = p*c*q*dx(mesh) - p*c0*q*dx(mesh) - lmbda*dt*dot(D*grad(mu_mid),grad(q))*dx(mesh) + dt*dot(C*grad(c), grad(q))*dx(mesh)
L1 = mu*v*dx(mesh) - p*f*v*dx(mesh) - dot(D*grad(c), grad(v))*dx(mesh)
L = L0 + L1
# Compute directional derivative about u in the direction of du (Jacobian)
a = derivative(L, u, du)
# Set constant values for BCs, first for chemical potential
inletmu = Constant(2.0)
outletmu = Constant(0)
# Set constant values for BCS, now for phase field. These depend on both the mu
# conditions and the exact form of the expression for f(c) above
inletc = Constant(np.sqrt(1.0+np.sqrt(2.0/3.0)))
outletc = Constant(-1.0)
# Impose the BCs on labelled facets. This could be done more tersely with renumbered
# boundaries, but for the moment we label them manually
bc1m = DirichletBC(ME.sub(1), inletmu, boundaries, 2)
bc2m = DirichletBC(ME.sub(1), inletmu, boundaries, 3)
bc3m = DirichletBC(ME.sub(1), outletmu, boundaries, 4)
bc4m = DirichletBC(ME.sub(1), outletmu, boundaries, 5)
bc1c = DirichletBC(ME.sub(0), inletc, boundaries, 2)
bc2c = DirichletBC(ME.sub(0), inletc, boundaries, 3)
bc3c = DirichletBC(ME.sub(0), outletc, boundaries, 4)
bc4c = DirichletBC(ME.sub(0), outletc, boundaries, 5)
bcs = [bc1m, bc2m, bc3m, bc4m, bc1c, bc2c, bc3c, bc4c]
# Create nonlinear problem and Newton solver
problem = CahnHilliardEquation(a, L, bcs)
solver = NewtonSolver()
# To solve with MPI we use parallelised solver and preconditioner
solver.parameters["linear_solver"] = "gmres"
solver.parameters["relative_tolerance"] = 1e-6
solver.parameters["preconditioner"] = "hypre_amg"
# Output file
file = File("macro.pvd", "compressed")
# Step in time
t = 0.0
T = 50*dt
file << (u, t)
while (t < T):
t += dt
u0.vector()[:] = u.vector()
solver.solve(problem, u.vector())
file << (u, t)
plot(u.split()[0])
interactive()
Note that the inlet is marked as boundaries 2 and 3, whilst the outlet is marked as 4 and 5. The other walls are 6 and 7, and this is confirmed by plotting the boundary meshfunction.
Edit:
I have also tried this approach with flux boundary conditions, and the problem is the same - the solution 'blows up' near the inlet.
