Cahn-Hilliard Euler not working

Hi! I’m trying to simulate Cahn-Hilliard equation using explicit Euler scheme. However, no matter what changes I make and what parameters I choose I can’t get it to work. Can I use this scheme or does the non-linearity of the equation strictly require Newton iterations or some similar method? I also tried Euler semi-implicit method and implicit but nothing works. Any help is greatly appreciated.

real radius = 0.05;
real M=1;
real centerX = 0.35;
real centerX2 = 0.5;
real centerY = 0.5;
real Inside = 1;
real Outside = -1;
real Pe= 62;
real Ch=0.01;

func real initialCondition(real x, real y)
    if (sqrt(pow(x - centerX, 2) + pow(y - centerY, 2)) <= radius || 
       sqrt(pow(x - centerX2, 2) + pow(y - centerY, 2)) <= radius || 
       sqrt(pow(x - centerX, 2) + pow(y - 0.8, 2)) <= radius || 
       sqrt(pow(x - centerX2, 2) + pow(y - 0.8, 2)) <= radius)
    return Inside;
    return Outside;

mesh Th=square(72,72);
fespace Vh(Th,P1);
fespace Vh2(Th,P1);

real dt=0.01;
int i,k;  

Vh2 u,w,phi,oldv,oldz,psi;

plot (u,wait=1,cmm="Cahn-Hilliard",value=true);

for (int i=0;i<30;i++)

    solve CahnHill1(w,phi) = 
     int2d(Th)(Ch*Ch*(dx(oldv)*dx(phi)+dy(oldv)*dy(phi)) )+ int2d(Th)((oldv^3-oldv)*phi) -int2d(Th)(w*phi);

    solve CahnHill2(u,phi)= int2d(Th)(u*phi) + int2d(Th)(dt*1/Pe*(dx(w)*dx(phi)+dy(w)*dy(phi)) ) -int2d(Th)(oldv*phi);

      plot (u,wait=1,fill=true, cmm="Cahn-Hilliard_"+i,value=true);


Could you upload the code or put on github? Copy/paste has been a problem.

Is this better?

I’ve never used github before

Probably the biggest problem is the time step is too big. There may be others.
The attached seems to run just slightly out of bounds as min and max exceed
}!}. You can do w and u si,multaneously unless your analysis suggested they
may be better in sequence. Of course generally you want to do more things
at once. Uinsg “old” is a concession to the linearlity constrtain. You may be able
to replace some of the “oldv” things with “u” but again you should look at the
discretized form details.

diff gistfile1.edp gistfile1.txt_orig 
< //real Outside = 0;
< real Pe=  62;
< real Ch= 0.01;
> real Pe= 62;
> real Ch=0.01;
< int sz=72*3;
< mesh Th=square(sz,sz);
< plot(Th,cmm="mesh",wait=1);
> mesh Th=square(72,72);
> plot(Th,wait=1);
< real dt=1e-8;
> real dt=0.01;
< Vh2 phiw;
< for (int i=0;i<300;i++)
> for (int i=0;i<30;i++)
< //    solve CahnHill1(w,phi) = 
< //     int2d(Th)(Ch*Ch*(dx(oldv)*dx(phi)+dy(oldv)*dy(phi)) )+ int2d(Th)((oldv^3-oldv)*phi) -int2d(Th)(w*phi);
< //    solve CahnHill2(u,phi)= int2d(Th)(u*phi) + int2d(Th)(dt*1/Pe*(dx(w)*dx(phi)+dy(w)*dy(phi)) ) -int2d(Th)(oldv*phi);
<     solve CahnHill1(u,w,phi,phiw) = 
<      int2d(Th)(Ch*Ch*(dx(oldv)*dx(phiw)+dy(oldv)*dy(phiw)) )+ int2d(Th)((oldv^3-oldv)*phiw) -int2d(Th)(w*phiw)
<     + int2d(Th)(u*phi) + int2d(Th)(dt*1/Pe*(dx(w)*dx(phi)+dy(w)*dy(phi)) ) -int2d(Th)(oldv*phi);
>     solve CahnHill1(w,phi) = 
>      int2d(Th)(Ch*Ch*(dx(oldv)*dx(phi)+dy(oldv)*dy(phi)) )+ int2d(Th)((oldv^3-oldv)*phi) -int2d(Th)(w*phi);
>     solve CahnHill2(u,phi)= int2d(Th)(u*phi) + int2d(Th)(dt*1/Pe*(dx(w)*dx(phi)+dy(w)*dy(phi)) ) -int2d(Th)(oldv*phi);
<       plot (u,wait=0,fill=true, cmm="Cahn-Hilliard_"+i,value=true);
>       plot (u,wait=1,fill=true, cmm="Cahn-Hilliard_"+i,value=true);

gistfile1.edp (1.4 KB)

See if this is any more of what you expect I changes the initial conditions to a
random thing and you can watch it segregate somewhat. I thought there
may be an issue with length scale and see also the slightly different
formulation which produce similar “u” as your’s except I probably dropped
some scale factors.

gistfile1.edp (2.1 KB)

Hi! Thank you for taking a look at my code and making the changes you did, it’s really working now! I tried putting everything into one solve command instead of two and fixing the min and max bounds but I still got a solution that diverged so I though the main problem was somewhere in the formulation. Now I see from your suggestions that it was all of those issues + time step that was too large+ some other issues you fixed. Thank you once again for your help!

If you are doing any serious work you should probably actually
try to find some analytical cases to test. Should the solution stay bounded?
When do the time derivatives go to zero? etc.

Once I add Stokes equation to this one I can compare the results to an article I found online. So I’ll see how that goes. But maybe I can try to find one with just the C-H.

if its free full text can you post the link?

I take it this is the article citation you sent but is not full text,

although it is cited a few times and curious if this is your overall interest,

Just skimming it one comment was about the interface thickness
and this comes up in a lot of places like with electrochemistry
and etching a solid in a liquid as a lot goes on in a few nm
of the interface.