% Numerically calculate soln. to advection eqn. dq/dt + dq/dx = 0 on
% domain 0<x<1 with periodic BCs and smooth and square-wave ICs. Compare 
% Lax-Wendroff (LW) and MC methods at various grid spacings. 
% Plot 2-norm error comparison.

  % User-selectable parameters

  L = 1;  % Domain size
  c = 1;
  tf = 1; % Final time 
  mu = 0.8; % Courant number
  p = 3:6; % dx =  L/ 2.^p 
  
  % Select initial condition (1=sine wave, 2 = square wave)
  ic = 1;
  
  switch ic
    case 1
      q0 = inline('sin(2*pi*x)');  % sine wave
      q0text = 'sin(2\pix)';
    case 2
      q0 = inline('sign(0.25-abs(mod(x,1)-0.4))');   % square wave
      q0text = 'square wave';
  end
  %===========================================
  np = length(p);
  h = zeros(1,np); % Initialize vector of grid spacings
  err_LW = zeros(1,np); % Initialize error matrix for LW method
  err_MC = zeros(1,np); % Initialize error matrix for spectral method

  for ip = 1:np
    nx = 2^p(ip);
    dx = L/nx;
    x = dx*((1:nx) - 0.5); % x cell-centers
    dt = mu*dx/c; % Timestep
    nt = round(tf/dt); % Number of timesteps to take
    qLW = q0(x); % Initial state
    qMC = q0(x); % Initial state

    % Rather than using ghost cells, we define vectors jp1 (= j+1 except 1 at
    % cell nx) and jm1 (= j-1, except nx for cell 1)
    % that account for the periodic BCs and update all cells in one vector
    % operation. 

    jp1 = [2:nx 1];
    j   = 1:nx;
    jm1 = [nx 1:(nx-1)];
  
    for it = 1:nt
      %LW slope
      sigLW = (qLW(jp1)-qLW(j))/dx;
      
      %MC slope
      sig1 = (qMC(jp1)-qMC(jm1))/(2*dx); % centered slope
      sig2 = 2*(qMC(j)-qMC(jm1))/dx;     % 2*left slope
      sig3 = 2*(qMC(jp1)-qMC(j))/dx;     % 2*right slope

      % Calculate sig = minmod(sig1,sig2,sig3)
      
      sigsum = sign(sig1) + sign(sig2) + sign(sig3);
      allpos = sigsum>2;
      allneg = sigsum<-2;
      posneg = abs(sigsum)<2;
      sigMC = zeros(1,nx);
      if(~isempty(sig1(allpos))) 
        sigMC(allpos) = min(sig1(allpos),sig2(allpos));
	    sigMC(allpos) = min(sigMC(allpos),sig3(allpos));
      end
      if(~isempty(sig1(allneg))) 
        sigMC(allneg) = max(sig1(allneg),sig2(allneg));
        sigMC(allneg) = max(sigMC(allneg),sig3(allneg));
      end
      if(~isempty(sig1(posneg))) 
        sigMC(posneg) = 0;
      end

      fh= c*(qLW+ 0.5*dx*sigLW.*(1 - mu));          % LW Flux at j+0.5
      qLW(j) = qLW(j) - (dt/dx)*(fh(j) - fh(jm1));  % Update qLW

      fh= c*(qMC+ 0.5*dx*sigMC.*(1 - mu));          % MC Flux at j+0.5
      qMC(j) = qMC(j) - (dt/dx)*(fh(j) - fh(jm1));  % Update qLW
      
    end
    
    qe = q0(x-c*tf);  % Exact solution at time tf
    err_LW(ip) = sqrt(dx)*norm(qLW-qe);  
    err_MC(ip) = sqrt(dx)*norm(qMC-qe);  
    h(ip) = dx;  
    
    % Plot solutions for coarsest resolution
    
    if(ip==1)
      subplot(2,2,2*ic-1)
      plot(x(j),qLW(j),'rx',x(j),qMC(j),'bx')
      xlabel('x')
      ylabel(['q(x,' num2str(tf) ')'])
      ylim([-1.25 1.25])
      text(0.8,1.1,['nx = ' num2str(nx)])
      title(['q_t+q_x=0, \mu = 0.8, q_0 = ' q0text]) 

      hold on
      xe = (0:0.01:1)*L;
      plot(xe,q0(xe-c*tf),'g--')
      legend('LW','MC','exact')
      
      % Must calculate slopes at tf
      
      %LW slope
      sigLW = (qLW(jp1)-qLW(j))/dx;
      
      %MC slope
      sig1 = (qMC(jp1)-qMC(jm1))/(2*dx); % centered slope
      sig2 = 2*(qMC(j)-qMC(jm1))/dx;     % 2*left slope
      sig3 = 2*(qMC(jp1)-qMC(j))/dx;     % 2*right slope

      % Calculate sig = minmod(sig1,sig2,sig3)
      
      sigsum = sign(sig1) + sign(sig2) + sign(sig3);
      allpos = sigsum>2;
      allneg = sigsum<-2;
      posneg = abs(sigsum)<2;
      sigMC = zeros(1,nx);
      if(~isempty(sig1(allpos))) 
        sigMC(allpos) = min(sig1(allpos),sig2(allpos));
	    sigMC(allpos) = min(sigMC(allpos),sig3(allpos));
      end
      if(~isempty(sig1(allneg))) 
        sigMC(allneg) = max(sig1(allneg),sig2(allneg));
        sigMC(allneg) = max(sigMC(allneg),sig3(allneg));
      end
      if(~isempty(sig1(posneg))) 
        sigMC(posneg) = 0;
      end


      for i = 1:nx
        plot([x(i)-0.5*dx x(i)+0.5*dx],...
             [qLW(i)-0.5*sigLW(i)*dx qLW(i)+0.5*sigLW(i)*dx],'r')
        plot([x(i)-0.5*dx x(i)+0.5*dx],...
             [qMC(i)-0.5*sigMC(i)*dx qMC(i)+0.5*sigMC(i)*dx],'b')
      end
      hold off
    end

  end

  subplot(2,2,2*ic)
  loglog(h,err_LW,'rx',h,err_MC,'bx')
  xlabel('dx')
  ylabel('grid 2-norm of error')
  legend('LW','MC',4)
  hold on

  % Add dx^2 (sine wave) and dx^0.5 (square wave) fit lines
  
  switch ic
    case 1
      csine = 20;   
      plot(h,csine*h.^2,'k--')
      text(.07,.3,'dx^2','Color','k')
    case 2
      csqwave  = 1;
      plot(h,csqwave*h.^0.25,'k--')
      ylim([0.1 1])
      text(.07,.8,'dx^{0.25}','Color','k')
  end

  hold off
