# A problem of stability

**URL:** https://community.freefem.org/t/a-problem-of-stability/139
**Category:** General Discussion
**Created:** [September 17, 2019, 4:59am UTC](https://community.freefem.org/t/a-problem-of-stability/139 "2019-09-17T04:59:12Z")
**Posts on this page:** 3
**Page:** 1

<div class="post-metadata">

### Author: ![xlzheng](https://avatars.discourse-cdn.com/v4/letter/x/2acd7d/32.png) [@xlzheng](https://community.freefem.org/u/xlzheng)
#### Post date: [September 17, 2019, 4:59am UTC](https://community.freefem.org/t/a-problem-of-stability/139/1 "2019-09-17T04:59:13Z")

</div>

Hi everyone,  
I write a program to simulate the cavity flow’s linear temporal growth rates and the neutral curve, but the results aren’t ideal,

 ![growth%20rates](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/b9e5228df2a9022bcde35d0f85b613e0220971ab.jpeg) ![neutral%20curve](https://canada1.discourse-cdn.com/flex030/uploads/freefem/original/1X/a4dd84a903413398ccb326fc4855687857f1480f.jpeg)

I put the program of linear temporal growth rates as flowing, looking forward your help! Thanks very much!

load “msh3”  
load “tetgen”  
load “medit”  
int N=70;  
real [int,int] coordinate((N+1)_(N+1),2);  
for(int i=0;i\<=N;i++)  
{  
for(int j=0;j\<=N;j++)  
{  
coordinate(i_(N+1)+j,0)=(cos(j_pi/N)+1)/2;  
coordinate(i_(N+1)+j,1)=(cos(i_pi/N)+1)/2;  
cout\<\<coordinate(i_(N+1)+j,0)\<\<" "\<\<coordinate(i\*(N+1)+j,1)\<\<endl;  
}  
}  
mesh Th2d=triangulate(coordinate(:,0),coordinate(:,1));

int[int] r1=[1,10],r2=[1,20];  
mesh Th=change(Th2d,label=r1);  
savemesh(Th,“Th2d.mesh”);

fespace Xh(Th,P2); //for velocity  
fespace Mh(Th,P1); //for pressure and stress  
fespace XhXhXhMh(Th,[P2,P2,P2,P1]);  
fespace XhXhXh(Th,[P2,P2,P2]);  
Xh Ux,Uy,Vx,Vy,dUx,dUy;  
Mh P,dP,Q;  
Ux=(y==1);  
Uy=(x==0); Uy=-Uy; //The boundary conditions  
Xh ux,uy,uz,ux0,uy0,uz0;  
Xh uxold,uyold,vx,vy,vz;  
Mh p,p0;  
Mh dp,pp,pold,q;

// Physical parameter  
real re=500, rebegin=0, reend=1000, Ha=0, k=0.9;  
complex shift=0.+0i, omega;  
real nu= 1./re,errold,abeigen,abeigenold;  
real eps=1e-6;  
real err=0;  
complex[int] Bu1(XhXhXhMh.ndof),u(XhXhXhMh.ndof);  
complex[int] Au1(XhXhXhMh.ndof),Au2(XhXhXhMh.ndof);

varf op([ux,uy,uz,p],[vx,vy,vz,q])=  
int2d(Th)(-vx\*(uxold_dx(ux)+ux_dx(uxold)+uyold_dy(ux)+uy_dy(uxold))  
- vy\*(uxold_dx(uy)+ux_dx(uyold)+uyold_dy(uy)+uy_dy(uyold))  
- vz\*(uxold_dx(uz)+uyold_dy(uz))  
+ p\*(dx(vx)+dy(vy)-1i_k_vz)  
- nu\*(dx(ux)_dx(vx)+dy(ux)dy(vx)+dx(uy)dx(vy)+dy(uy)dy(vy)+dx(uz)dx(vz)+dy(uz)dy(vz))  
- nukk(uxvx+uyvy+uz_vz)  
+ q\*(dx(ux)+dy(uy))-1i_k_uz_q  
+ p_q_1e-18  
)  
+ int2d(Th)(1i_shift\*(ux_vx+uy_vy+uz\*vz))  
+ on(10,ux=0,uy=0,uz=0);

varf a([ux,uy,uz,p],[vx,vy,vz,q])=  
int2d(Th)(-vx\*(uxold_dx(ux)+ux_dx(uxold)+uyold_dy(ux)+uy_dy(uxold))  
- vy\*(uxold_dx(uy)+ux_dx(uyold)+uyold_dy(uy)+uy_dy(uyold))  
- vz\*(uxold_dx(uz)+uyold_dy(uz))  
+ p\*(dx(vx)+dy(vy)-1i_k_vz)  
- nu\*(dx(ux)_dx(vx)+dy(ux)dy(vx)+dx(uy)dx(vy)+dy(uy)dy(vy)+dx(uz)dx(vz)+dy(uz)dy(vz))  
- nukk(uxvx+uyvy+uz_vz)  
+ q\*(dx(ux)+dy(uy))-1i_k_uz_q  
+ p_q\*1e-18  
)  
+ on(10,ux=0,uy=0,uz=0);

varf b([ux,uy,uz,p],[vx,vy,vz,q])=  
int2d(Th)(1i\*(ux_vx+uy_vy+uz\*vz))  
;

//find the first piont

for(int t=0;t\< 20;t++)  
{  
solve Baseflow([dUx,dUy,dP],[Vx,Vy,Q])=  
int2d(Th)( Vx\*(dUx_dx(Ux)+dUy_dy(Ux))+Vy\*(dUx_dx(Uy)+dUy_dy(Uy))  
+Vx\*(Ux_dx(dUx)+Uy_dy(dUx))+Vy\*(Ux_dx(dUy)+Uy_dy(dUy))  
+nu\*(dx(dUx)_dx(Vx)+dy(dUx)dy(Vx))  
+nu(dx(dUy)dx(Vy)+dy(dUy)dy(Vy))  
-dP(dx(Vx)+dy(Vy))  
-Q(dx(dUx)+dy(dUy))  
+1e-8_dP_Q  
)  
-int2d(Th)(Vx_(Ux_dx(Ux)+Uy_dy(Ux))+Vy\*(Ux_dx(Uy)+Uy_dy(Uy))  
+nu\*(dx(Ux)\*dx(Vx)+dy(Ux)_dy(Vx))  
+nu_(dx(Uy)_dx(Vy)+dy(Uy)dy(Vy))  
-P(dx(Vx)+dy(Vy))  
-Q_(dx(Ux)+dy(Uy))  
)  
+on(10,dUx=0,dUy=0)  
;  
cout\<\<“t=”\<\<t\<\<" dU.max “\<\<dUx[].max\<\<” “\<\<dUy[].max\<\<” "\<\<"dP.max= “\<\<dP[].max\<\<endl;  
err= dUx[].linfty + dUy[].linfty + dP[].linfty;  
Ux[]-=dUx[];  
Uy[]-=dUy[];  
P[]-=dP[];  
cout \<\< t \<\< " err = " \<\< err \<\< " " \<\< eps \<\< " rey =” \<\< re \<\< endl;  
if(err \< eps) break; // converge  
if( t\>10 && err \> 10.) break; // Blowup ???  
}  
if(err \< eps)  
{  
uxold=Ux;  
uyold=Uy;

```
//stability		
int kkmax=0;
real erreigen;
		
matrix<complex> OP= op(XhXhXhMh,XhXhXhMh,solver=sparsesolver); 
matrix<complex> A= a(XhXhXhMh,XhXhXhMh,solver=sparsesolver);
matrix<complex> B= b(XhXhXhMh,XhXhXhMh,solver=CG); 
int nev=40;
	
complex[int] ev(nev); //to store nev eigen value
XhXhXhMh<complex>[int] [eux,euy,euz,ep](nev); //to store nev eigen vector
int ncv=4*nev;
	
int kk=EigenValue(OP,B,sigma=shift,value=ev,vector=eux, tol=1e-20,maxit=0,ncv=ncv);
kk=min(kk,nev); // some time the number of converged eigen value
					// can be greater than nev;
abeigenold=-1;	
for (int i=0;i<kk;i++)
{
	//kkmax=(imag(ev[i])>imag(ev[i+1])?i:(i+1));
	//kkmaxold=(imag(ev[kkmax])>imag(ev[kkmaxold])?kkmax:kkmaxold);
		
	ux=eux[i]; 
	uy=euy[i];
	uz=euz[i];
	p=ep[i];
	u=[ux[],uy[],uz[],p[]];
	complex v=ev[i];
	Bu1=B*u;
	Bu1=v*Bu1;
	Au1=A*u;
	u=Au1-Bu1;
	erreigen=norm(u.sum);
	abeigen=imag(ev[i]);
	cout<<"re= "<<re<<" i= "<<i<<" ev[i]= "<<ev[i]<<" erreigen= "<<erreigen<<endl;
	if(abeigen>abeigenold)
	{
		abeigenold=abeigen;kkmax=i;
		cout <<" kkmax= "<<kkmax <<endl;				
	}
}
omega=ev[kkmax];
ux=eux[kkmax]; 
uy=euy[kkmax];
uz=euz[kkmax];
p=ep[kkmax];
cout<<"re= "<<re<<" omega= "<<omega<<endl;
//shift=omega;		
if(imag(omega)>0&&abs(imag(omega))>1e-5)
{
	k=k-0.05;
	//shift=omega; cout<<"restart re"<<endl;
}
else if(imag(omega)<0&&abs(imag(omega))>1e-5)
{
	k=k+0.05;
	//shift=omega;	cout<<"restart re"<<endl;
}
else
{cout<<"sucesse to find the first point"<<endl;
 cout<<"k="<<k<<" re="<<re<<endl;}

```

}  
else  
{ // if blowup, increase nu (more simple)  
cout \<\< " restart re = " \<\< re \<\<" converge wrong"\<\<endl;  
}

real alp,alpbegin=0,alpend=5,err3;  
int itmax=150, iter;  
k=k+0.1;  
while (k\<alpend)  
{  
real erreigen;int kkmax=0;  
shift=0+0i;  
cout\<\<“shift=”\<\<shift\<\<endl;  
matrix OP= op(XhXhXhMh,XhXhXhMh,solver=sparsesolver);  
matrix A= a(XhXhXhMh,XhXhXhMh,solver=sparsesolver);  
matrix B= b(XhXhXhMh,XhXhXhMh,solver=CG);  
int nev=4;

```
complex[int] ev(nev); //to store nev eigen value
XhXhXhMh<complex>[int] [eux,euy,euz,ep](nev); //to store nev eigen vector
int ncv=4*nev;
	
int kk=EigenValue(OP,B,sigma=shift,value=ev,vector=eux, tol=1e-20,maxit=0,ncv=ncv);
kk=min(kk,nev); // some time the number of converged eigen value
					// can be greater than nev;
abeigenold=-1;	
for (int i=0;i<kk;i++)
{
	//kkmax=(imag(ev[i])>imag(ev[i+1])?i:(i+1));
	//kkmaxold=(imag(ev[kkmax])>imag(ev[kkmaxold])?kkmax:kkmaxold);
		
	ux=eux[i]; 
	uy=euy[i];
	uz=euz[i];
	p=ep[i];
	u=[ux[],uy[],uz[],p[]];
	complex v=ev[i];
	Bu1=B*u;
	Bu1=v*Bu1;
	Au1=A*u;
	u=Au1-Bu1;
	erreigen=norm(u.sum);
	abeigen=imag(ev[i]);
	cout<<"re= "<<re<<" i= "<<i<<" ev[i]= "<<ev[i]<<" erreigen= "<<erreigen<<endl;
	if((abeigen>abeigenold)&&(abs(real(v))<1e-5))
	{
		abeigenold=abeigen;kkmax=i;
		cout <<" kkmax= "<<kkmax <<endl;				
	}
}
omega=ev[kkmax];
ux=eux[kkmax]; 
uy=euy[kkmax];
uz=euz[kkmax];
p=ep[kkmax];
cout<<"k= "<<k<<" re= "<<re<<" omega= "<<omega<<endl;
k=k+0.1;	

```

}

---

<div class="post-metadata">

### Author: ![D\_N\_Sarkar](https://avatars.discourse-cdn.com/v4/letter/d/0ea827/32.png) [@D\_N\_Sarkar](https://community.freefem.org/u/D_N_Sarkar)
#### Post date: [October 17, 2020, 10:44am UTC](https://community.freefem.org/t/a-problem-of-stability/139/2 "2020-10-17T10:44:07Z")

</div>

I am also working on similar problems. Did you find a solution to your problem? Any help would be appreciated. Looking forward to your kind reply.

---

<div class="post-metadata">

### Author: ![D\_N\_Sarkar](https://avatars.discourse-cdn.com/v4/letter/d/0ea827/32.png) [@D\_N\_Sarkar](https://community.freefem.org/u/D_N_Sarkar)
#### Post date: [October 20, 2020, 12:26pm UTC](https://community.freefem.org/t/a-problem-of-stability/139/3 "2020-10-20T12:26:16Z")

</div>

I am trying to carry out 3-D hydrodynamic stability analysis of 2-dimensional lid-driven square cavity flow at Re = 200.

I have calculated the 2-dimensional base flow at Re=200 and validated the base flow results with existing literature.

I would like to calculate the eigenvalues for the above-mentioned problem. I am not able to implement the stability analysis calculation in FreeFem.

There is a hydrodynamic stability problem implemented in ‘SLEPc-complex’ avaiable in FreeFem website.

I would like to implement my hydrodynamic stability code without any packages like SLEPc and PETSc. I am not able to find enough guidelines on implementation of hydrodynamic stability problem in FreeFem.

I would be extremely grateful if you please guide me on the implementation of hydrodynamic stability problem in FreeFem.

Looking forward to your kind reply.