# Problem with C-H +Stokes + viscosity coupling
**URL:** https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633
**Category:** General Discussion
**Created:** [August 2, 2023, 4:28pm UTC](https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633 "2023-08-02T16:28:51Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@bATTLEGROUND](https://community.freefem.org/u/bATTLEGROUND)
#### Post date: [August 2, 2023, 4:28pm UTC](https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633/1 "2023-08-02T16:28:51Z")
Hi! I’m trying to simulate the mixing of two liquids of different viscosities using Cahn-Hilliard and Stokes equations. Due to the non-linear nature of C-H equation I’m using Newton iterations. If the viscosity ratio (lambda) is 1 the code works without any problems but for any other ratio after around 100 iterations I start getting these minimums and maximums that spread and take over the entire domain as shown in the picture. Any ideas what’s wrong with my code? Thanks!
```auto
real radius = 0.05;
real M=1;
real centerX = 0.35;
real centerX2 = 0.5;
real centerY = 0.5;
real Inside = 1.0;
real Outside = -1.0;
real Pe= 10000;
real Ch=0.01;
func real initialCondition(real x, real y)
{
if (y<0.5)
return Inside;
else
return Outside;
}
mesh Th=square(72,72);
plot(Th,wait=1);
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);
real lambda=10;
real a=0.01;
real dt=0.01;
int i,k;
func real df(real u) { return u^3-u; }
func real ddf(real u) { return 3*u^2-1; }
Vh2 [u,w],[v,z],[oldu,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha; // f'(u)
Vh ddfalpha; // f"(u)
fespace Uh(Th, P1);
Uh U, V;
Uh UU, VV;
fespace Ph(Th, P1);
Ph p, pp;
varf vdJ([v,z],[phi,psi]) = int2d(Th)
((u-oldu)*phi+dt*(dx(u)*U+dy(u)*V)*phi-dt*1/Pe*(dx(oldz)*dx(phi)+dy(oldz)*dy(phi))+oldz*psi-
Ch*Ch*(dx(oldu)*dx(psi)+dy(oldu)*dy(psi))-dfalpha*psi);
varf vhJ([v,z],[phi,psi]) = int2d(Th)
(-v*phi-dt*1/Pe*(dx(z)*dx(phi)+dy(z)*dy(phi))+z*psi-Ch*Ch*(dx(v)*dx(psi)
+dy(v)*dy(psi))-ddfalpha*v*psi);
[u,w]=[initialCondition(x,y),0];
plot (u,wait=1,cmm="Cahn-Hilliard",value=true);
macro SGrad(u,v) [[dx(u),0.5*(dx(v)+dy(u))],[0.5*(dx(v)+dy(u)),dy(v)]] //
macro div(u1,u2) (dx(u1)+dy(u2)) //
macro mu(a) (0.5*(a+1)-1/lambda*0.5*(a-1)) //
for (int i=0;i<150;i++)
{
oldu[]=u[];
solve Stokes( [U,V,p],[UU,VV,pp]) =
int2d(Th)( 2* mu(oldu)*(SGrad(U,V):SGrad(UU,VV)) - div(U,V)*pp)
- int2d(Th)( div(UU,VV)*p)
+ on(3, 2, 4, U=0, V=0)
+ on(1, U=16*(x)^2*(1-x)^2, V=0)
;
for (int k=0;k<10;k++) // boucle de Newton
{
dfalpha = df( oldu ) ;
ddfalpha = ddf( oldu );
auxv[]=vdJ(0,Vh2);
real res= auxv[]'*auxv[];
cout << i << " residu^2 = " << res << endl;
matrix H;
if( res< 1e-12) break;
H= vhJ(Vh2,Vh2,factorize=1,solver=LU);
vv[]=H^-1*auxv[];
v[] -= vv[];
oldu[]=v[];
}
u[]=v[];
plot (u,wait=1,fill=true, cmm="Cahn-Hilliard"+i,value=true);
}
```
---
### Author:  [@cmd](https://community.freefem.org/u/cmd)
#### Post date: [August 2, 2023, 5:19pm UTC](https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633/2 "2023-08-02T17:19:27Z")
Your Peclet number is huge. Try with a smaller Pe value or dramatically increase the resolution of your mesh. Also, you should use an inf-sup stable finite element discretization like P2-P1.
---
### Author:  [@bATTLEGROUND](https://community.freefem.org/u/bATTLEGROUND)
#### Post date: [August 2, 2023, 5:41pm UTC](https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633/3 "2023-08-02T17:41:23Z")
Thank you so much for the reply! I’ll try out those suggestions. If I can ask one more thing, did I add the advection correctly considering I’m using Newton iterations? I don’t have much experience with it so I’m not sure if I should add anything to vhJ?
EDIT: P2-P1 discretization and bigger mesh resolution fixed the problem. Thank you!
---
### Author:  [@Monirul25](https://community.freefem.org/u/Monirul25)
#### Post date: [February 8, 2024, 5:40am UTC](https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633/4 "2024-02-08T05:40:11Z")
//Mesh
mesh Th=square(100,100);
//border C(t=0,2_pi){x=cos(t);y=sin(t);};
//mesh Th=buildmesh (C(50));
plot(Th,wait=1);
macro Grad(u)[dx(u),dy(u)]//
macro uex(t)(exp(-2_t)_(sin(pi_x)^2)_(sin(pi_y)^2)) //
//Fespace
fespace Vh(Th,P1);
fespace Vh2(Th,[P1,P1]);
//Vh uex;
//Parameter
real a=0.01;
real dt=0.10;
int i,k;
//real t;
// Function
func real df(real u){return u^3-u; }
func real ddf(real u){return 3\*u^2-1; }
Vh2 [u,w],[v,z],[oldv,oldz],[auxv,auxz],[vv,zz];
Vh dfalpha;
Vh ddfalpha;
//Vh uex;
//macro uex(t) (exp(-2_t)_(sin(pi_x)^2)_(sin(pi\*y)^2)//
// Problem
varf vdJ([v,z],[phi,psi])=int2d(Th)
((oldv-u)_phi+dt_(dx(oldz)_dx(phi)+dy(oldz)dy(phi))+oldzpsi-
a_(dx(oldv)\*dx(psi)+dy(oldv)_dy(psi))-dfalpha_psi);
varf vhJ([v,z],[phi,psi])=int2d(Th)
(v_phi+dt_(dx(z)_dx(phi)+dy(z)dy(phi))+zpsi-a_(dx(v)_dx(psi)
+dy(v)dy(psi))-ddfalphav_psi);
//u=uex(t);
[u,w]=[(sin(pi_x)^2)_(sin(pi_y)^2),-2_a\*(pi^2)_(sin(pi^2_(x^2+y^2))
-2\*(pi^2)_(sin(pi_x)^2)_(sin(pi_y)^2)_(x^2+y^2))+((sin(pi_x)^2)^3)_((sin(pi_y)^2)^3)
-(sin(pi_x)^2)_(sin(pi\*y)^2)];
//u=uex;
plot(u,wait=1,cmm=“Cahn-Hilliard” ,value=true);
//u=uex(t);
//macro uex(t) (exp(-2_t)_(sin(pi_x)^2)_(sin(pi_y)^2))//
// u=uex(t);
//uex(t)=exp(-2_t)_(sin(pi_x)^2)_(sin(pi_y)^2);
//real t=0;// start from t=0
real t;
[//t=0.1](https://t=0.1);
for (int i=0;i\<1./dt; i++) // Time loop
{
t=t+dt;
//u=uex(t);
//=exp(-2_t)_((sin(pi_x)^2)_(sin(pi_y)^2));
{
oldv[]=u[];
for (int k=0;k\<10;k++)// Newton loop
{
dfalpha=df( oldv );
ddfalpha=ddf( oldv );
auxv[]=vdJ(0,Vh2);
real res= auxv[]'auxv[]; // Residual
cout \<\< i \<\< " residu^2 = " \<\< res \<\< endl;
matrix H;
if (res\< 1e-12) break;
H=vhJ(Vh2,Vh2, factorize=1,solver=LU);
vv[]=H^-1auxv[]; // Newton
v[] -=vv[];
oldv[]=v[];
}
u[]=v[];
//u=uex(t);
//uex=exp_sin(pi_x)^2_sin(pi_y)^2;
plot(u,wait=1,cmm=“Cahn-Hilliard”,value=true);
//uex(t)=exp(-2_t)_sin(pi_x)^2_sin(pi_y)^2;
//plot(uex(t),wait=1,cmm=“Cahn-Hilliard exact solution”,value=true);
```
}
cout << "t" << t << "L^2-error=" << sqrt(int2d(Th)((u-uex(t))^2))<
### Author:  [@Monirul25](https://community.freefem.org/u/Monirul25)
#### Post date: [February 8, 2024, 5:42am UTC](https://community.freefem.org/t/problem-with-c-h-stokes-viscosity-coupling/2633/5 "2024-02-08T05:42:03Z")
Can you tell me how i post a queries in this platform?