A problem of stability

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,

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(jpi/N)+1)/2;
coordinate(i
(N+1)+j,1)=(cos(ipi/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*(uxolddx(ux)+uxdx(uxold)+uyolddy(ux)+uydy(uxold))
- vy*(uxolddx(uy)+uxdx(uyold)+uyolddy(uy)+uydy(uyold))
- vz*(uxolddx(uz)+uyolddy(uz))
+ p*(dx(vx)+dy(vy)-1ikvz)
- 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))
- nu
k
k
(ux
vx+uy
vy+uz
vz)
+ q*(dx(ux)+dy(uy))-1ikuzq
+ p
q1e-18
)
+ int2d(Th)(1i
shift*(uxvx+uyvy+uz*vz))
+ on(10,ux=0,uy=0,uz=0);

varf a([ux,uy,uz,p],[vx,vy,vz,q])=
int2d(Th)(-vx*(uxolddx(ux)+uxdx(uxold)+uyolddy(ux)+uydy(uxold))
- vy*(uxolddx(uy)+uxdx(uyold)+uyolddy(uy)+uydy(uyold))
- vz*(uxolddx(uz)+uyolddy(uz))
+ p*(dx(vx)+dy(vy)-1ikvz)
- 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))
- nu
k
k
(ux
vx+uy
vy+uz
vz)
+ q*(dx(ux)+dy(uy))-1ikuzq
+ 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*(uxvx+uyvy+uz*vz))
;

//find the first piont

for(int t=0;t< 20;t++)
{
solve Baseflow([dUx,dUy,dP],[Vx,Vy,Q])=
int2d(Th)( Vx*(dUxdx(Ux)+dUydy(Ux))+Vy*(dUxdx(Uy)+dUydy(Uy))
+Vx*(Uxdx(dUx)+Uydy(dUx))+Vy*(Uxdx(dUy)+Uydy(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
dPQ
)
-int2d(Th)(Vx
(Uxdx(Ux)+Uydy(Ux))+Vy*(Uxdx(Uy)+Uydy(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;	

}