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].
