% Script advect_MC.m
%
% Use MC flux-limited finite volume method to solve the advection
% equation with periodic BCs and a square-wave IC:
%  
%   psi_t + c*psi_x = 0, 0 < x < L
%   psi(L,t) = psi(0,t)
%   psi(x,0) = H(0.25L-|x-0.5L|)  (H = Heaviside step function, 1 for
%                                   arg>0, 0 for arg < 0, 0.5 for arg = 0.
%
% with grid spacing dx = L/M and Courant number mu = c*dt/dx.  Output is
% plot of solution at times 0, T.
%=======================================================================
  
  % Input parameters
  
  c = 1;
  L = 1;
  nx = 20;
  mu = 0.5;
  T = 1;
  
  % Derived parameters

  dx = L/nx;
  dt = mu*dx/c;
  nt = round(T/dt);
  
  % Initial condition
  
  phi = zeros(1,nx);
  x = dx*((1:nx) - 0.5);
  isone = abs(x-0.5*L)<0.25*L;
  phi(isone) = 1;

  subplot(2,1,1)
  plot(x,phi,'b--')
  xlabel('x')
  ylabel('\phi')
  ylim([-0.2 1.2])
  
  % Implement MC method with mean value phi(j) and slope sig(j)
  % in each cell j = 1:nx
  % 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. 
  
  for it = 1:nt
    jp1 = [2:nx 1];
    j   = 1:nx;
    jm1 = [nx 1:(nx-1)];

    sig1 = (phi(jp1)-phi(jm1))/(2*dx); % centered slope
    sig2 = 2*(phi(j)-phi(jm1))/dx;     % 2*left slope
    sig3 = 2*(phi(jp1)-phi(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;
    sig = zeros(1,nx);
    if(length(sig1(allpos)) > 0) 
	sig(allpos) = min(sig1(allpos),sig2(allpos));
	sig(allpos) = min(sig(allpos),sig3(allpos));
    end
    if(length(sig1(allneg)) > 0) 
      sig(allneg) = max(sig1(allneg),sig2(allneg));
      sig(allneg) = max(sig(allneg),sig3(allneg));
    end
    if(length(sig1(posneg)) > 0) 
      sig(posneg) = 0;
    end

    fh= c*(phi+ 0.5*dx*sig.*(1 - mu));            % Flux at j+0.5
    phi(j) = phi(j) - (dt/dx)*(fh(j) - fh(jm1));  % Update phi
  end
  hold on
  plot(x,phi,'r-')
  legend('t=0',['t=' num2str(T)])
  title(['MC method on \psi_t + ' num2str(c) '\psi_x = 0, periodic BCs'])
  hold off
