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
On Thursday, 23 October 2025 at 13:01:57 UTC+1 Irene Moreno wrote:
> 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.
>
> [image: ICs_phi.png][image: snapshot_phi_0.006737946999085467.png][image:
> snapshot_phi_0.19907123629213125.png][image:
> snapshot_phi_449.7430549622601.png]
> 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
> #######################################################################
--
To unsubscribe from this group, send email to [email protected]
View this message at https://list.nist.gov/fipy
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
