# Cahn-Hilliard Euler not working
**URL:** https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643
**Category:** General Discussion
**Created:** [August 4, 2023, 9:05pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643 "2023-08-04T21:05:34Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@bATTLEGROUND](https://community.freefem.org/u/bATTLEGROUND)
#### Post date: [August 4, 2023, 9:05pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/1 "2023-08-04T21:05:34Z")
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.
```auto
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;
else
return Outside;
}
mesh Th=square(72,72);
plot(Th,wait=1);
fespace Vh(Th,P1);
fespace Vh2(Th,P1);
real dt=0.01;
int i,k;
Vh2 u,w,phi,oldv,oldz,psi;
u=initialCondition(x,y);
plot (u,wait=1,cmm="Cahn-Hilliard",value=true);
for (int i=0;i<30;i++)
{oldv[]=u[];
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);
}
```
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 5, 2023, 11:03am UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/2 "2023-08-05T11:03:41Z")
Could you upload the code or put on github? Copy/paste has been a problem.
Thanks.
---
### Author:  [@bATTLEGROUND](https://community.freefem.org/u/bATTLEGROUND)
#### Post date: [August 5, 2023, 1:13pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/3 "2023-08-05T13:13:24Z")
Is this better?
>
I’ve never used github before
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 5, 2023, 6:08pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/4 "2023-08-05T18:08:19Z")
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.
```auto
diff gistfile1.edp gistfile1.txt_orig
8,10c8,9
< //real Outside = 0;
< real Pe= 62;
< real Ch= 0.01;
---
> real Pe= 62;
> real Ch=0.01;
22,24c21,23
< int sz=72*3;
< mesh Th=square(sz,sz);
< plot(Th,cmm="mesh",wait=1);
---
>
> mesh Th=square(72,72);
> plot(Th,wait=1);
28c27
< real dt=1e-8;
---
> real dt=0.01;
32c31
< Vh2 phiw;
---
>
37c36
< for (int i=0;i<300;i++)
---
> for (int i=0;i<30;i++)
40,49c39,40
< // 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);
50a42
> 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);
52c44
< 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](https://community.freefem.org/uploads/short-url/fPeAyMYeefwTgA1KfGmod0c419r.edp) (1.4 KB)
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 6, 2023, 1:56pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/5 "2023-08-06T13:56:11Z")
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](https://community.freefem.org/uploads/short-url/3ixuABrIUOE6MmWYjazAbZTc77P.edp) (2.1 KB)
---
### Author:  [@bATTLEGROUND](https://community.freefem.org/u/bATTLEGROUND)
#### Post date: [August 6, 2023, 4:09pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/6 "2023-08-06T16:09:50Z")
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!
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 6, 2023, 5:20pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/7 "2023-08-06T17:20:20Z")
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.
---
### Author:  [@bATTLEGROUND](https://community.freefem.org/u/bATTLEGROUND)
#### Post date: [August 6, 2023, 7:44pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/8 "2023-08-06T19:44:54Z")
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.
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 6, 2023, 8:46pm UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/9 "2023-08-06T20:46:23Z")
if its free full text can you post the link?
Thanks
---
### Author:  [@marchywka](https://community.freefem.org/u/marchywka)
#### Post date: [August 7, 2023, 11:58am UTC](https://community.freefem.org/t/cahn-hilliard-euler-not-working/2643/10 "2023-08-07T11:58:17Z")
I take it this is the article citation you sent but is not full text,
[https://www.sciencedirect.com/science/article/abs/pii/S0009250905008043?via%3Dihub](https://www.sciencedirect.com/science/article/abs/pii/S0009250905008043?via%3Dihub)
although it is cited a few times and curious if this is your overall interest,
[https://www03.core.ac.uk/download/pdf/231814262.pdf](https://www03.core.ac.uk/download/pdf/231814262.pdf)
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.