%%% Linear Computational Instability for Diffusion Equation
%
%  d_phi/dt = kappa d^2 phi/dx^2                          (1)
%
%  The general finite difference equation is
%
%  (phi_j - phi_j^0)/dt = kappa {                         (2)
%         theta [ phi_(j+1) - 2 phi_j + phi_j-1) ]/ dx^2
%    + (1-theta) [ phi_(j+1)^0 - 2 phi_j^0 + phi_j-1)^0 ]/ dx^2
%                               }
%  where phi_j are at future time, phi_j^0 are known values
%
%  Take spatial Discrete Fourier Transform (assume infinite domain)
%
%  PHI(k) = sum_{j=-inf}^inf  [ exp(ik jdx) phi_j ] dx
%
%  Take DFT of (2) and shift dummy indices j+1 -> m, j-1-> n
%    (recall that the sum is from -inf to inf) to recover standard DFT
%     from j+1 and j-1 terms.
%
%   PHI(k) -PHI^0(k) = (kappa dt/dx^2) {
%      theta [ PHI(k) exp(-ikdx)   - 2 PHI(k)   + PHI(k) exp(ikdx)  ]
% + (1-theta)[ PHI_0(k) exp(-ikdx) - 2 PHI^0(k) + PHI^0(k) exp(ikdx)]
%                                      }
%  Collecting terms
%
%           PHI(k) = PHI^0(k) F(k)
%
%   where F(k) is the transfer function for a time step dt
%   for wavenumber k
%
%  F(k) = [ 1 - 2(1-theta)(kappa dt/dx^2) (1 - cos(k dx) ) ] / ...
%         [ 1  +  2 theta (kappa dt/dx^2) (1 - cos(k dx) )  ] 
%
%  The Nyquist wavenumber k_N = pi/dx
%  is highest wavenumber representable on a uniform grid
% -------------------------------------------------------
%
%   Also demonstrate linear computational instability in time domain
%
%                                       Ed Waddington    2008-05-11
% --------------------------------------------------------------------
%
%  set up diffusivity k 
    kappa = 1;
%
%  set up grid length Lx 
    Lx = 20;
%
%  set up discretization dx, and nodes vector x_vec 
    dx    = 1.0;
    x_vec = 0:dx:Lx;
    nx    = length(x_vec);
%
%
%  set duration Lt of time stepping
  %  Lt = 10;
%
%  set time step
    dt = input('OK, what is dt for this test?   ');
 %   nt = fix( Lt/dt ) + 1;
%
%  print out stability parameter  beta = kappa dt/dx^2
%  ----------------------------------------------------
    beta = 2*kappa*dt/dx^2;
    stability_string = num2str( 2*kappa*dt/dx^2);
%
    message_string ='beta = [ 2 kappa dt/dx^2 ] = ';
    panic_message = [ message_string  stability_string ];
  %
    disp( panic_message )
%
%  Set time-step mixing fraction theta for spatial terms
%      i.e. theta*phi + (1-theta)*phi_0
%      theta = 1    Fully Implicit
%            = 0.5  Crank-Nicolson
%            = 0    Explicit
% ------------------------------
%
%   theta should exceed (beta-1)/2 beta 
%  ----------------------------------------
    theta_test_string = ('For stability: choose theta > (beta - 1) / (2 beta) =  ');
    theta_test_message = [ theta_test_string num2str( (beta - 1)/(2*beta) ) ];
    disp( theta_test_message );
%
    theta = input('OK what value for time-partition parameter theta?   ');
%
    theta_label_string ='theta = ';
    theta_value_string =num2str( theta );
    theta_message = [ theta_label_string  theta_value_string ];

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Now set up wavenumber parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    k_N = pi/dx;    % Nyquist wavenumber 
%
    n_k = 51;      % number of points in wavenumber domain
%
    dk = k_N/(n_k - 1);
%
    k_vec = 0: dk: k_N;
%
%  Evaluate transfer function for continuous diffusion equation
%  -------------------------------------------------------------
    F_analytic = exp(- kappa * k_vec.^2 * dt );
%
%  Evaluate transfer function for discrete diffusion equation
%  ------------------------------------------------------------
    F = ( 1 - (1-theta)*(2*kappa*dt/dx^2)*(1 - cos(k_vec*dx) ) ) ./ ...
        ( 1 + theta*(2*kappa*dt/dx^2)*(1 - cos(k_vec*dx) ) );
%
%  show kappa*dt/dx^2
%  -------------------
   %  disp('2*kappa*dt/dx^2  = '); disp(2*kappa*dt/dx^2)
%
%  Plot transfer functions F_analytical(k) and F(k)
%  -------------------------------------------------
%
    figure(1)
    plot(k_vec, F_analytic, 'Color', 'k' )
    hold on;grid on;box on
    xlabel('Wavenumber  k')
    ylabel('Transfer function F(k) for time-step dt')
%
%  put white patch to cover old dt and theta messages
    patch( [ 0.05*pi 0.05*pi 0.2*pi 0.2*pi ], ...
           [ 0.1 0.25 0.25 0.1 ], 'w' ) 
%
%  write current dt and theta onto plot 
    dt_string = ['dt = ' num2str(dt) ];
    text(0.06*pi,0.22, dt_string )
% 
    theta_string = ['\theta = ' num2str(theta) ];
    text(0.06*pi,0.12, theta_string )
%
    title('Black - Analytical Transfer function for time step dt')
%
    color = input('What color of line to plot _this_ time?    ');    
    plot(k_vec, F, 'Color', color )
%   
%
%  Plot transfer function for repeated steps to give
%    total elapsed time t=1
%    When F(k) < 0, and 1/dt is not an integer,
%    the product of transfer functions needed to represent passage
%    of t=1 may be complex.  However, I plot only the magnitude.
%  ---------------------------------------------------------------
    figure(2)
    plot(k_vec, F_analytic.^(1/dt), 'Color', 'k' )
    hold on;grid on;box on
    xlabel('Wavenumber  k')
    title('Multi-step Transfer function F(k) to simulate t=1')
    ylabel('Transfer function  F(k)')
    text(0.7*pi,0.85, 'repeated to simulate t=1' )
%
    plot(k_vec, abs( F.^(1/dt) ), 'Color', color )
%
%
%  put white patch to cover old dt and theta messages
    patch( [ 0.02*pi 0.02*pi 0.2*pi 0.2*pi ], ...
           [ 0.1 0.4 0.4 0.1 ], 'w' ) 
%
%  write current dt and theta onto plot 
    dt_string = ['dt = ' num2str(dt) ];
    text(0.06*pi,0.32, dt_string )
% 
    theta_string = ['\theta = ' num2str(theta) ];
    text(0.06*pi,0.12, theta_string )
%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  end of wavenumber calculation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Now find finite-difference solution in time domain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  How long to run simulation?
%  set duration Lt of time stepping and number nt of time steps
%  ------------------------------------------------------------
    Lt = input('How long to run this simulation?   ');
    nt = fix( Lt/dt ) + 1;
%
%  set string for plot-ends time
    time_ends = [ ' Run ends at '  num2str(Lt) ];
%
%
%  set up initial condition  (sum of a few sinusoids)
%   The b_N term puts in a little ripple at Nyquist wavenumber.
%   Does that sinificantly "feed" runaway growth?
%  ---------------------------------------------------
    a_0 = 10;
    a_1=  1.0;
    b_1 = 2.0;
    b_2 = 2.0;
    b_3 = 2.0;
    b_4 = 2.0;
    b_N = 0.0;
    phi_P_0 = a_0 + a_1 * x_vec/Lx ...
                            + b_1 * sin( pi*x_vec/Lx ) ...
                            + b_2 * sin(2*pi*x_vec/Lx) ...
                            + b_3 * sin(3*pi*x_vec/Lx) ...
                            + b_4 * sin(4*pi*x_vec/Lx) ...
                            + b_N * cos((nx-1)*pi*x_vec/Lx);
%
%
%  Boundary Conditions 
%  --------------------
%   Set to hold steady from initial condition.
%
%   specified value at x = 0
     phi_0 = a_0 + b_N;
%
%   specified gradient at x = Lx
     phi_prime_Lx = a_1/Lx ...
                    + b_1 *   (pi/Lx) * cos(  pi*x_vec(end)/Lx ) ...
                    + b_2 * (2*pi/Lx) * cos(2*pi*x_vec(end)/Lx ) ...
                    + b_3 * (3*pi/Lx) * cos(3*pi*x_vec(end)/Lx ) ...
                    + b_4 * (4*pi/Lx) * cos(4*pi*x_vec(end)/Lx ) ...
                    + b_N * (nx*pi/Lx) * cos(nx*pi*x_vec(end)/Lx );
%
%
%  set up shifted phi arrays for phi_E and phi_W
    phi_D_0 = [ phi_P_0(2:end)   0  ];
    phi_U_0 = [ 0  phi_P_0(1:end-1) ];
%
%
%  set up coefficient vectors
    a_D = repmat( (theta*kappa*dt/dx^2), nx, 1)';
    a_U = a_D;
    a_P = ( 1 + a_D + a_U );
%
    a_E = zeros( size(a_P) );
    a_W = a_E;
%
    a_D_0 = repmat( ( (1-theta)*kappa*dt/dx^2), nx, 1)';
    a_U_0 = a_D_0;
    a_P_0 = 1 - a_D_0 - a_U_0;
%
    b = a_P_0 .*phi_P_0 + a_D_0 .*phi_D_0 + a_U_0 .*phi_U_0;
%
%
%  Introduce Boundary Conditions into equations
%
%  fixed value at x_vec(1)
    a_P(1) = 1;
    a_U(1) = 0;
    a_D(1) = 0;
%
    a_P_0(1) = 0;
    a_U_0(1) = 0;
    a_D_0(1) = 0;
    b(1)   = phi_0;
%
%  fixed gradient at x_vec(end)
    a_P(end) =  1/dx;
    a_U(end) = -1/dx;
    a_D(end) = 0;
%
    a_P_0(end) = 0;
    a_U_0(end) = 0;
    a_D_0(end) = 0;

    b(end)   = phi_prime_Lx;
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   end of setup for time domain
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%S
%
%  Start plot with initial-condition profile
    fig_no = input('Which figure (>=3) to use for finite-difference solution?   ');
    figure(fig_no)
%
%   set y scale limits on plot
    phi_max = 25;
%
    plot(x_vec, phi_P_0, 'k')
    hold on;grid on;box on;
    xlabel('Depth  z')
    ylabel('Solution  \phi(z)')
    text( 0.1*Lx, 0.95*phi_max, panic_message )
    text( 0.1*Lx, 0.90*phi_max, theta_message )
    text( 0.1*Lx, 0.85*phi_max, theta_test_message);
    text( 0.1*Lx, 0.80*phi_max, time_ends )
    title('Stability tests with Finite Differences')
    set(gca,'ylim', [ 0, phi_max] )
%
    pause
%
%  Here it is!   Big time-step loop!
%  ---------------------------------
    for j = 1:nt-1
   %
   % First, update RHS with current "old" solution
   %  set up shifted phi arrays for phi_0_E and phi_0_W
       phi_D_0 = [ phi_P_0(2:end) 0 ];
       phi_U_0 = [ 0 phi_P_0(1:end-1) ];
   %
   % Update RHS vector, which depends on phi_0
       b = a_P_0 .*phi_P_0 + a_D_0 .*phi_D_0 + a_U_0 .*phi_U_0;
   %
   
   %  Re-introduce Boundary conditions into RHS vector
       b(1)   = phi_0;
       b(end) = phi_prime_Lx;
   %
   %
       a_P(end) = 1/dx;
       a_U(end) = 1/dx;
       a_D(end) = 0;
%
       a_P_0(end) = 0;
       a_U_0(end) = 0;
       a_D_0(end) = 0;


   %
   %  Lots of equations? Let's solve the suckers all together!
        phi = solver( a_E, a_W, a_U, a_D, a_P, b );
   %
   %  update "old" time step
        phi_P_0 = phi;
   %    
   %  Add current time step to plot  
        figure(fig_no)
        plot(x_vec, phi, 'Color', color )
   %
   %  show time step and elapsed time
        time_message = ['Current Time = ' num2str( j*dt) ];
        patch( Lx*[0.1 0.1 0.4 0.4], phi_max*[0.73 0.77 0.77 0.73 ], 'w' )
        text(0.1*Lx, 0.75*phi_max, time_message)
   %
   %
      pause
        
   %
    end
%
%
% print maximum value of phi - is it going wild?
   disp('Maximum of |phi| =    '); disp( max( abs(phi) ) );
% 
%
%