Temporal hydrodynamic stability

General Discussion
1

Here is the code to calculate temporal hydrodynamic stability for couette flow
real Mj=0.9; // mach number

complex I = 1i;
real ga=1.4; //gamma
real M2=(1/(gaMjMj));
real Re;
cout<<“Enter Reynold’s number”<<endl;
cin>>Re;
real om,alp,bet;
// defining the type of stability problem

int a,b;
cout<<“Enter 1 for spatial stability and 2 for temporal stability\n”<<endl;
cin>> a;

if(a==1)
{
cout<<“Enter wave frequency\n”<<endl;
cin>>om;}
else{
cout<<“Enter streamwise wave number\n”<<endl;
cin>>alp;
cout<<“Enter spanwise wave number\n”<<endl;
cin>>bet;}
mesh Th;
int m=200;
real h=1./m;
Th=square(1,m,[x*h,y]);
plot(Th,wait=1);

fespace Vh(Th,P2);
fespace Mh(Th,P1);
fespace Gh(Th,[P2,P2,P2]);
Vh W,u,v,w,a2,a3,Wy,a1;

W=y;

Wy=1;
plot(W,wait=1);
real lambda=(10^18);
matrixA1;
matrix B1;

              varf A([u,v,w],[a1,a2,a3])=int2d(Th)(a1*(I*bet*W+(alp^2+bet^2)/(Re))*u-(1/Re)*dy(a1)*dy(u)+a1*lambda*((alp^2)*u-I*alp*dy(v)+(alp*bet)*w)+
                                                         a2*(I*bet*W+(alp^2+bet^2)/(Re))*v-(1/Re)*dy(a2)*dy(v)-a2*I*lambda*(alp*dy(u)+bet*dy(w))-lambda*dy(a2)*dy(v)+
                                                         a3*v*Wy+a3*(I*bet*W+(alp^2+bet^2)/(Re))*w-(1/Re)*dy(a3)*dy(w)+a3*lambda*((alp*bet)*u-I*bet*dy(v)+(bet^2)*w)
                                                         
                                                     )+on(1,3,u=0,v=0,w=0);
                                                         
                                                         
              varf B([u,v,w],[a1,a2,a3])=int2d(Th)(a1*I*u+a2*I*v+a3*I*w); 
              
              
              
                A1=A(Gh,Gh);
               B1=B(Gh,Gh);
              ofstream file("A1.csv");

for (int i = 0; i < A1.m; i++) {
for (int j = 0; j < A1.n; j++) {
if(imag(A1(i,j))>0){file<<real(A1(i,j))<<“+”<<imag(A1(i,j))<<“i”;}
else{file<<real(A1(i,j))<<imag(A1(i,j))<<“i”;}
if (j < A1.m - 1) {
file<< “,”;
}
}
file << endl;
}

       ofstream file1("B1.csv");

for (int i = 0; i < B1.n; ++i) {
for (int j = 0; j < B1.m; ++j) {
if(imag(B1(i,j))>0){file1<<real(B1(i,j))<<“+”<<imag(B1(i,j))<<“i”;}
else{file1<<real(B1(i,j))<<imag(B1(i,j))<<“i”;}
if (j < B1.m - 1) {
file1<< “,”;
}
}
file1<< endl;
}
cout<<A1.m<<endl;
cout<<A1.n<<endl;
Is the way i defined the mesh and define varf is correct?