% ========================================================================
% 2D Inertialess shear-drive flow for viscoelastic fluids. 
% Find the eigenvalues of TTs using chebfun and integral formulation.
% 
% m-file needed: svdfr
%
% Original equations:
%
%   x'''' + a3*x''' + a2*x'' + a1*x' + a0*x = b1*d1' + b0*d2;
%
%   u = ii*kx*x;
%   v = -x';
%
%   BCs:
%   x1(-1) = x1(+1) = 0; x1'(-1) = x1'(+1) = 0;
%
% written by: Binh Lieu
%
% reference:
%
% B. K. Lieu and M. R. Jovanovic, "Computation of frequency responses
% for linear time-invariant PDEs on a compact interval", Journal of
% Computational Physics, 2013, doi:10.1016/j.jcp.2013.05.010;
% also arXiv:1112.0579
%
% ========================================================================

%% Changing Chebfun's operator preferences
cheboppref('maxdegree',10000)   % changing maximum degree of functions
cheboppref('restol',1e-6)       % changing residual tolerance
 
%% System parameters
kx = 1; % streamwise wavenumber
kx2 = kx*kx;
kx4 = kx2*kx2;
om = 0; % temporal frequency
 
beta = 0.5; % viscosity ratio
 
Weval = [1, linspace(10,100,10)]; % Weissenberg numbers
Wegrd = length(Weval);
 
dom = domain(-1,1);     % domain of your function
y = chebfun('y',dom);
fone = chebfun(1,dom);      % one function
fzero = chebfun(0,dom);     % zero function
 
%% Boundary conditions
Wa0{1} = [1, 0, 0, 0; 0, 1, 0, 0];
Wb0{1} = [1, 0, 0, 0; 0, 1, 0, 0];
 
%% Looping over Weissenberg number
Smax = zeros(Wegrd,1); A0 = cell(1,1); B0 = cell(1,2); C0 = cell(2,1);
for indWe = 1:Wegrd
    We = Weval(indWe);
    
    % coefficients of the operators A0 and B0
    a0 = (kx4*(1 + 1i*om + 1i*kx*We*y).*(beta - ...
        (2*(-1 + beta)*We*We*(1 + 2*We*We))./((1 + 1i*om + ...
        1i*kx*We*y).^3) - ((-1 + beta)*(1 + 2*We*We))./(1 + ...
        1i*om + 1i*kx*We*y)))./(1 + beta*(1i*om + 1i*kx*We*y));
    
    a1 = (2*1i*(-1 + beta)*kx2*kx*We*(1i*om + ...
        1i*kx*We*y).*(1 + 1i*om - 2*We*We + ...
        1i*kx*We*y))./(((1 + 1i*om + 1i*kx*We*y).^2).*(1 + ...
        beta*(1i*om + 1i*kx*We*y)));
    
    a2 = (2*kx2*(-1 - (1i*om + 1i*kx*We*y).*(2 + 1i*om + ...
        1i*om*We*We + 1i*kx*(We + (We^3)).*y + ...
        beta*(1 + (1i*om + 1i*kx*We*y).*(2 + 1i*om - We*We + ...
        1i*kx*We*y)))))./(((1 + 1i*om + ...
        1i*kx*We*y).^2).*(1 + beta*(1i*om + 1i*kx*We*y)));
    
    a3 = (2*(-1 + beta)*kx*We*(-1i*1i*om + kx*We*y))./((1 + 1i*om + ...
        1i*kx*We*y).*(1 + beta*(1i*om + 1i*kx*We*y)));
    
    b1 = -((1 + 1i*om + 1i*kx*We*y))./(1 + beta*(1i*om + 1i*kx*We*y));
    
    b0 = 1i*kx*((1 + 1i*om + 1i*kx*We*y))./(1 + beta*(1i*om + ...
        1i*kx*We*y));
    
    % system operators
    % a0*phi + a1*phi? + a2*phi?? + a3*phi??? + 1*phi???? =
    % % 0*d1 + b1*d1? + b0*d2 + 0*d2?
    A0 = {[a0, a1, a2, a3, fone]};

    B11 = [fzero, b1]; B12 = [b0, fzero];
    B0 = {B11, B12};
    
    C11 = [fzero, fone];
    C21 = [-1i*kx*fone, fzero];
    C0 = {C11; C21};
    
    % determining the principal singular pair of the frequency response 
    % operator
    [Sfun,Sval] = svdfr(A0,B0,C0,Wa0,Wb0,1,1);
    
    % storing the largest singular value for each Weissenberg number
    Smax(indWe) = Sval(1);

end

%% Plotting the largest singular value as a function of We
plot(Weval,Smax,'x-','LineWidth',1.1,'MarkerSize',10);
xlab = xlabel('We', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin');

% Plotting the principal singular functions for We = 100
% principal singular function corresponding to the streamwise velocity
u = Sfun{1};
% principal singular function corresponding to the wall-normal velocity
v = Sfun{2};

% Plotting the most amplified streamwise velocity profile
plot(y,real(u)/norm(real(u)),'-','LineWidth',1.1);
xlab = xlabel('y', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin');

% Plotting the most amplified wall-normal velocity profile
plot(y,real(v)/norm(real(v)),'-','LineWidth',1.1);
xlab = xlabel('y', 'interpreter', 'tex');
set(xlab, 'FontName', 'cmmi10', 'FontSize', 20);
h = get(gcf,'CurrentAxes');
set(h,'FontName','cmr10','FontSize',15,'xscale','lin','yscale','lin');