% exemples de numérique

clear all; format compact; clf

% paramètres
N=50 % nombre de points de grille
L=1 % hauteur du domaine

%%%% matrices de dérivation
scal=2/(L);
[y,DM] = chebdif(N,2); 
D1=scal*DM(:,:,1); D2=scal^2*DM(:,:,2); 
y=(y)/scal; 

I=eye(N); Z=0*I;



thecase=5
%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%
switch thecase
  case 1
    disp('advection-diffusion en 1D')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    u=1; mu=0.1;
    A=-u*D1+mu*D2; E=I;

    % conditions limites
    loc=[1,N]; 
    C=I(loc,:);
    E(loc,:)=0; A(loc,:)=C;
    
    % modes propres
    [U,S]=eig(A,E);
    s=diag(S); rem=abs(s)>1e3|isnan(s); s(rem)=[]; U(:,rem)=[];
    [t,o]=sort(-real(s)); s=s(o); U=U(:,o); 
    
    % graphiques
    subplot(1,2,1)
    for ind=1:length(s)
        h=plot(real(s(ind)),imag(s(ind)),'*');
        ttt=['subplot(1,2,2);num=' num2str(ind) ';plot(y,U(:,num));grid on'];
        set(h,'buttondownfcn',ttt);
        hold on
    end
    grid on; hold off
    



  case 2
    disp('diffusion en 2D')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    mu=0.1; alpha=1;
    x=linspace(0,2*pi/alpha,50);
    [X,Y]=meshgrid(x,y);
    
    A=mu*(-alpha^2*I+D2); E=I;

    % conditions limites
    loc=[1,N]; 
    C=I(loc,:);
    E(loc,:)=0; A(loc,:)=C;
    
    % modes propres
    [U,S]=eig(A,E);
    s=diag(S); rem=abs(s)>1e3|isnan(s); s(rem)=[]; U(:,rem)=[];
    [t,o]=sort(-real(s)); s=s(o); U=U(:,o); 
    
    % graphiques
    subplot(1,2,1)
    for ind=1:length(s)
        h=plot(real(s(ind)),imag(s(ind)),'*');
        t0=['subplot(1,2,2);num=' num2str(ind) ';'];
        t1='mesh(X,Y,real(U(:,num)*exp(i*alpha*x)));';
        t2='xlabel(''x''); ylabel(''y'');grid on; axis equal; axis tight';
        set(h,'buttondownfcn',[t0 t1 t2]);
        hold on
    end
    grid on; hold off

    
    
    
  case 3
    disp('stokes en 2D')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    alpha=1; rho=1;mu=1; 
    x=linspace(0,2*pi/alpha,20);
    [X,Y]=meshgrid(x,y);
    
    dy=D1; dyy=D2; 
    dx=i*alpha*I; dxx=-alpha^2*I;
    
    % matrices
    delta=dxx+dyy; 
    E=[rho*I, Z, Z; Z, rho*I, Z; Z, Z, Z]; 
    A=[mu*delta, Z, -dx; Z, mu*delta, -dy; dx, dy, Z]; 
    
    % conditions limites
    III=eye(3*N); 
    loc=[1, N, N+1, 2*N]; 
    C=III(loc,:); 
    E(loc,:)=0; A(loc,:)=C; 

    % modes propres
    [U,S]=eig(A,E);
    s=diag(S); rem=abs(s)>1e3|isnan(s); s(rem)=[]; U(:,rem)=[];
    [t,o]=sort(-real(s)); s=s(o); U=U(:,o); 
    
    % graphiques
    subplot(1,2,1)
    for ind=1:length(s)
        h=plot(real(s(ind)),imag(s(ind)),'*');
        t0=['subplot(1,2,2);num=' num2str(ind) ';'];
        t1=['quiver(X,Y,real(U(1:N,num)*exp(i*alpha*x)),real(U(N+1:2*N,num)*exp(i*alpha*x)));'];
        t2='xlabel(''x''); ylabel(''y'');grid on; axis equal; axis tight';
        set(h,'buttondownfcn',[t0 t1 t2]);
        hold on
    end
    grid on; hold off
    
    
    
  case 4
    disp('fluide avec écoulement de base')
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    alpha=1; Re=200;
    x=linspace(0,2*pi/alpha,20);
    [X,Y]=meshgrid(x,y);
    
    %%%% écoulement de base: couche de cisaillement
    u=0.5*tanh(2*y)+1.5; % base flow 
    up=1./cosh(2*y).^2;% first derivative of base flow

% $$$     %%%% écoulement de base: poiseuille
% $$$     u=1-y.^2; % base flow 
% $$$     up=-2*y;% first derivative of base flow
    
    dy=D1; dyy=D2; 
    dx=i*alpha*I; dxx=-alpha^2*I;
    delta=dxx+dyy; 
    
    % matrices
    S=-diag(u)*dx+delta/Re; 
    A=[S, -diag(up), -dx; Z, S, -dy; dx, dy, Z]; 
    E=blkdiag(I,I,Z); 

    % conditions limites
    III=eye(3*N); 
    loc=[1, N, N+1, 2*N]; 
    C=III(loc,:); 
    E(loc,:)=0; A(loc,:)=C; 
    
    %modes propres
    [U,S]=eig(A,E);
    s=diag(S); rem=abs(s)>1e3|isnan(s); s(rem)=[]; U(:,rem)=[];
    [t,o]=sort(-real(s)); s=s(o); U=U(:,o); 
    
    % graphiques
    subplot(1,2,1)
    for ind=1:length(s)
        h=plot(real(s(ind)),imag(s(ind)),'*');
        t0=['subplot(1,2,2);num=' num2str(ind) ';'];
        t1=['quiver(X,Y,real(U(1:N,num)*exp(i*alpha*x)),real(U(N+1:2*N,num)*exp(i*alpha*x)));'];
        t2='xlabel(''x''); ylabel(''y'');grid on; axis equal; axis tight';
        set(h,'buttondownfcn',[t0 t1 t2]);
        hold on
    end
    grid on; hold off
    
    
    
  case 5
    disp('Rayleigh-Benard')
    %%%%%%%%%%%%%%%%%%%%%%
    Pr=10 % Prandtl number
    Ra=2000 % Rayleigh number
    rho=1,  % top fluid density and viscosity
    g=10 % gravity
    k=1 % thermal diffusivity
    mu=Pr*rho*k % fluid viscosity
    dt=1 % difference de temperature haut/bas
    di=Ra*mu*k/(rho*g*dt*L^3) % thermal dilatation
    gt=-dt/L % base field temperature gradient in y
    alpha=3 % wavenumber

    x=linspace(0,2*pi/alpha,20);
    [X,Y]=meshgrid(x,y);

    dy=D1; dyy=D2; dx=i*alpha*I;  dxx=-alpha^2*I;  lap=dxx+dyy;
    z=zeros(1,2*N); 

    % matrices
    A=[mu*lap, Z, -dx, Z; ...
       Z, mu*lap, -dy, rho*g*di*I; ...
       dx,dy,Z,Z; ...
       Z,-gt*I,Z, k*lap];     
    
    E=[rho*I,Z,Z,Z; ...
       Z,rho*I,Z,Z; ...
       Z,Z,Z,Z; ...
       Z,Z,Z,I]; 
    
    % conditions limites
    loc=[1 N N+1 2*N 2*N+1 3*N 3*N+1 4*N]; 
    II=blkdiag(I,I,I,I);
    C=[II(loc,:)];   
    A(loc,:)=C; E(loc,:)=0;
   
    % modes propres
    [U,S]=eig(A,E); 
    s=diag(S); rem=abs(s)>1e3|isnan(s); s(rem)=[]; U(:,rem)=[];
    [t,o]=sort(-real(s)); s=s(o); U=U(:,o); 
 
    % graphiques
    subplot(1,2,1)
    for ind=1:length(s)
        h=plot(real(s(ind)),imag(s(ind)),'*');
        t0=['subplot(1,2,2);num=' num2str(ind) ';'];
        t1=['quiver(X,Y,real(U(1:N,num)*exp(i*alpha*x)),real(U(N+1:2*N,num)*exp(i*alpha*x)));'];
        t2=['hold on; contour(X,Y,real(U(3*N+1:4*N,num)*exp(i*x*alpha))); xlabel(''x''); ylabel(''y'');'];
        t3='grid on; hold off; axis equal; axis tight';
        set(h,'buttondownfcn',[t0 t1 t2 t3]);
        hold on
    end
    grid on; hold off

end
