% Convection Diffsion Example for Class
%  THP Nov 29, 2005
%
%
%  Solves:  du/dt + a du/dx = mu d^2 u /dx^2 + 2 a 
%        Initial Condition:
%           u(x,0) = 2x + sin(2 pi x)
%        Boundary Conditions:
%           u(0,t) = 0.0   u(1,t) = 2.0
%        Final Steady State Solution:   
%           u = 2 x 
%
%        Sign(a) arbritary  i.e. Use +- splitting on a
%        Uses: 
%            1st order upwind for the convection part
%            2nd order central for the diffusion part
%        The three methods from the lectures are used:
%            Euler Explicit:  u_{n+1} = u_n + h (u')_n
%            Euler Ixplicit:  u_{n+1} = u_n + h (u')_{n+1}
%            Predictor-Corrector :
%                  P:  ut_{n+1/2} = u_n + h/2 (u')_n
%                  C:  u_{n+1}    = u_n + h   (ut')_{n+1/2}
%

clear;

%  Inputs:
M = 50; %  Number of Intervals for 0 < x < xmax
M = Enter('M',M);

a = 1.0; %  Wave speed, can be + or -
a = Enter('a',a);

mu = 1.e-2;  % Diffusion Coefficient
mu = Enter('mu',mu);

xmax = 1.0; % Domain Size
dx = xmax/(M);

x = [0:M]*dx;  % domain

nmax = 100;  % Number of iterations
nmax = Enter('nmax',nmax);

cfl_ee = 0.5;  %  CFL for Euler Explicit
cfl_ee = Enter('cfl for Euler Explicit',cfl_ee);
h_ee = cfl_ee/(abs(a)/dx + mu/dx^2);
% thp h_ee = min(cfl_ee/abs(a)*dx,cfl_ee/mu*dx^2);

cfl_ei = 10; %  CFL for Euler Implicit
cfl_ei = Enter('cfl for Euler Implicit',cfl_ei);
h_ei = cfl_ei/(abs(a)/dx + mu/dx^2);
% thp h_ei = min(cfl_ei/abs(a)*dx,cfl_ei/mu*dx^2);

cfl_pc = 0.5; %  CFL for Predictor-Corrector
cfl_pc = Enter('cfl for Predictor-Corrector',cfl_pc);
h_pc = cfl_pc/(abs(a)/dx + mu/dx^2);
% thp h_pc = min(cfl_pc/abs(a)*dx,cfl_pc/mu*dx^2);

disp([ 'h Euler Explicit = ' num2str(h_ee)]);
disp([ 'h Euler Implicit = ' num2str(h_ei)]);
disp([ 'h Predictor-Corrector = ' num2str(h_pc)]);

%  Split a = a^+ + a^-
ap = 0.5*(a + abs(a));
am = 0.5*(a - abs(a));

R = 0.0*x;  % Initialize R = 0.0

u = 2.0*x + sin(2.0*pi*x); % Initial Condition

u_ss = 2.0*x;  % Steady-State Solution

%  Initialize for each method variable
u_ee = u;
u_ei = u;
u_pc = u;

set (gcf,'doublebuffer','on');

for n = 1:nmax % loop for nmax times

  % Euler Explicit
  [u_ee,resid_ee(n)] = EE(M,a,mu,dx,ap,am,h_ee,u_ee);
  if resid_ee(n) > 1.0e5 
    disp([ 'resid overflow Euler Explicit: ABORT']);
    break;
  end

  % Euler INplicit
  [u_ei,resid_ei(n)] = EI(M,a,mu,dx,ap,am,h_ei,u_ei);
  if resid_ei(n) > 1.0e5 
    disp([ 'resid overflow Euler Implicit: ABORT']);
    break;
  end

  % Predictor-Corrector
  [u_pc,resid_pc(n)] = PC(M,a,mu,dx,ap,am,h_pc,u_pc);
  if resid_pc(n) > 1.0e5 
    disp([ 'resid overflow Predictor-Corrector: ABORT']);
    break;
  end

  % Plotting
  subplot(2,1,1);
  semilogy([1:n],resid_ee,'k',[1:n],resid_ei,'g',[1:n],resid_pc,'b');
  legend([ 'EE CFL = ' num2str(cfl_ee)], ...
      [ 'Ei CFL = ' num2str(cfl_ei)], ...
      [ 'PC CFL = ' num2str(cfl_pc)],'Location','SouthWest');
  xlabel('N');
  ylabel('Residual');

  title([ 'Convection-Diffusion Eq:  {\partial u / \partial t} ' ...
	' + a {\partial u / \partial x} = \mu {\partial^2 u / \partial x^2} + 2 a']);
  subplot(2,1,2);
  plot(x,u_ss,'ro',x,u_ee,'k',x,u_ei,'g',x,u_pc,'b');
  xlabel('x');
  ylabel('u');
  legend('Exact','EE','EI','PC','Location','NorthWest');
  title([ '   N = ' num2str(n)]);
  pause(0.1);

end