%{
# The Orr-Sommerfeld and Squire equation

Here we do just as in [poiseuille_uvp.m](), except that instead of using the primitive variables formulation $(u,v,p)$, we use the Orr-Sommerfeld and Squire equations, with the variables $v$ and $w$. We test it on the Poiseuille flow.
%}

clear all; clf

N=200; % number of grid points
alpha=1; % wavenumber in x
beta=0; % wavenumber in z
Re=10000; % Reynolds number

% differentiation matrices 
[y,DM] = chebdif(N,4); 
d.y=DM(:,:,1);    
d.yy=DM(:,:,2);    
d.yyyy=DM(:,:,4);    

I=eye(N);Z=zeros(N,N);II=blkdiag(I,I);
d.x=i*alpha*I;
d.z=i*beta*I;

% base flow and derivatives
u=1-y.^2;
up=-2*y;
upp=-2;

% laplacian
k2=alpha^2+beta^2;
lap=(d.yy-k2*I); 
lap2=(d.yyyy-2*k2*d.yy+k2^2*I); 

% system matrices
LOS=-diag(u)*lap*d.x+diag(upp)*d.x+lap2/Re;     
LSQ=-diag(u)*d.x+lap/Re;
LC= -diag(up)*d.z;

A=[LOS,Z;LC,LSQ]; 
E=[lap,I];

%{
# Boundary conditions

They are: zero at the walls for $v$ and $w$ (no-slip) and zero derivative at the walls for $v$ (that comes from the no-slip condition on $u$). The contraint matrix is thus
$$
C=\left(\begin{array}{c}
I|_{-1}&0\\
I|_{1}&0\\
0&I|_{-1}\\
0&I|_{1}\\
\partial_y|_{-1}&0\\
\partial_y|_{1}&0\\
\end{array}\right)
$$
such that the boundary conditions are implemented
$$
Cq=0
$$
where $q$ is the state vector
$$
q=\left(\begin{array}{c}
v\\
w
\end{array}\right)
$$
%}
% boundary conditions
v=1:N; w=v+N; % location vectors
loc=[v([1,2,N-1,N]),w([1,N])];
C=[II([v([1,N]),w([1,N])],:); ... % Dirichlet on v and w
    d.y([1,N],:), Z([1,N],:)];    % Neumann on v

A(loc,:)=C;
E(loc,:)=0;
    
% computing eigenmodes 
[U,S]=eig(A,E);
s=diag(S);  [t,o]=sort(-real(s)); s=s(o); U=U(:,o);
rem=abs(s)>1000; s(rem)=[]; U(:,rem)=[];

% loading the scanned figure for comparison and plotting
a=imread('poiseuille_spectra.png'); ss=size(a);
xx=linspace(0,1,ss(2));
yy=linspace(0,-1,ss(1));
image(xx,yy,a); axis xy
hold on

plot(-imag(s)/alpha,real(s),'ro')
grid on
axis([0,1,-1,0.1]);
xlabel('wave speed'); ylabel('exponential growth rate')

set(gcf,'paperpositionmode','auto');
print('-dpng','-r80','orr_sommerfeld_squire.png');

%{
![The spectrum](orr_sommerfeld_squire.png)

# Exercices/Contributions

* Please write down the equations and how to obtain them in the comments
* Please test it for a case with non-zero $\beta$ (an oblique wave)
* Please test it for the Blasius boundary layer
* Please do a quantitative comparison of the accuracy with the primitive variables formulation [poiseuille_uvp.m]().


%}