Cahn-Hilliard for the coupled evolution of two variables

[Irene Moreno]

Hi,

I have recently been working on an application of Cahn-Hilliard to a ternary system, and in the process of implementing the model, I have stumbled upon two main questions.

Firstly, I would like to know if the known problems of the coupled implementation will prevent me from developing the model. My system of equations is formed by two coupled sets of two equations each, so 4 eqs in total, with two of them solving for (following the example's notation) phi1, psi1, and the other set solving for phi2, psi2. However, the equations for phi1 and psi1 also depend on phi2 and psi2, respectively, and vice versa.

I have seen that Pysparse seems to take too much memory (is this the default solver?) and Trillinos cannot solve the coupled system, so I understand my best shot is PETSc (since my equations aren't in vector form). How would the implementation of this method be different from the code below?

############
c=0
while elapsed < duration:
    res = 1.0
    phi1.updateOld()
    phi2.updateOld()
    phi3.value = 1.0 - phi1 - phi2   # Unimportant variable, does not appear in the equations at all

    psi1.updateOld()
    psi2.updateOld()

    dt = min(100, numerix.exp(dexp))
    elapsed += dt
    dexp += 0.01

    while res > 1.0e-3:
        res = eq.sweep(dt=dt)
        print(res)
############

And secondly, about the evolution of phi in a binary system: I decided to work from the tutorial (link) towards my equation, as implementing my model directly led to issues. One of the first steps I took was, for the system in 1D, changing the ICs from gaussian noise with the mean in 0.5 to two-layered ones, with phi = 0.9999 at one side of the mesh, and phi = 0.2 in the other.  Full code for this described change is at the end.

Since these are very close to the minima displayed by the free energy density (same one as in the tutorial), I expected the system to just smooth out the interface a bit, with the bulk concentrations staying mainly where they were. However, what I have seen instead is, for some reason, the concentrations increasing/decreasing along the interface, going into values >1 before stabilising in the expected profile (see images below). If one of the ICs defines phi as 0.0001, the phi in that section of the mesh goes into values <0 before stabilising.


where the three first images correspond to the ICs and the very beginning of the simulation, and the last image is the final profile reached.

Since the free energy density used in the tutorial is a polynomial, it doesn't have any singularities, and therefore I understand that phi technically can go below 0 or above 1. In this case, is there a way to constrain the values the solver can accept in the solution? Would this be an acceptable fix, or would it interfere with the dynamics? Would a better approach be to use a log form in the free energy? 

Thank you very much in advance.

Kind regards,
Irene


#######################################################################
from fipy import CellVariable, Grid1D, GaussianNoiseVariable, DiffusionTerm, TransientTerm, ImplicitSourceTerm, Viewer
from fipy.tools import numerix
import numpy as np

nx = 20
mesh = Grid1D(nx=nx, dx=0.25)
lx = 0.25 * nx
height = mesh.cellCenters[0]

D = 1.
a = 1.
epsilon = 0.5

phi = CellVariable(name=r"$\phi$", mesh=mesh)
psi = CellVariable(name=r"$\psi$", mesh=mesh)

phi[:] = 1.0
phi[height > lx/2] = 0.2

psi.value[:] = a**2 * phi.value[:] * (1 - phi.value[:]) * (1 - 2 * phi.value[:])

viewer_phi = Viewer(vars=(phi,), datamin=-0.05, datamax=1.05)
viewer_psi = Viewer(vars=(psi,), datamin=-0.05, datamax=3.0)

viewer_phi.plot(f'ICs_phi.png')
viewer_psi.plot(f'ICs_psi.png')

dfdphi = a**2 * phi * (1 - phi) * (1 - 2 * phi)
d2fdphi2 = a**2 * (1 - 6 * phi * (1 - phi))

eq1 = (TransientTerm(var=phi) == DiffusionTerm(coeff=D, var=psi))
eq2 = (ImplicitSourceTerm(coeff=1., var=psi)
       == ImplicitSourceTerm(coeff=d2fdphi2, var=phi) - d2fdphi2 * phi + dfdphi
       - DiffusionTerm(coeff=epsilon**2, var=phi))

eq = eq1 & eq2

dexp = -5
elapsed = 0.
duration = 500.
c=0

while elapsed < duration:
    dt = min(100, numerix.exp(dexp))
    elapsed += dt
    dexp += 0.01
    eq.solve(dt=dt)
    if c % 25 == 0:
        viewer_phi.plot(f'snapshot_phi_{elapsed}.png')
        viewer_psi.plot(f'snapshot_psi_{elapsed}.png')
    c +=1
#######################################################################

[Irene Moreno]

Sorry, I forgot to include the link for the tutorial - I was looking at the 2D coupled CH one, https://pages.nist.gov/fipy/en/latest/generated/examples.cahnHilliard.mesh2DCoupled.html#module-examples.cahnHilliard.mesh2DCoupled

Kind regards,

Irene

[Irene Moreno]

Sorry again, forgot that the ICs of the tutorial example are ~1 and 0.2. To clarify, when I talked about the interface smoothing out since the compositions were near the free energy density minima, I was thinking about a system with ICs in 0.9999 and 0.0001, which I also have tried, and has also led to phi>1 and <0. For this system in 0.2 though, I understand why that layer ends in a phi value near 0, as it should.

Kind regards,

Irene