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 On Thursday, 23 October 2025 at 13:06:30 UTC+1 Irene Moreno wrote: > 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].
