% Finite volume spatial discretization to solve:
% u_t+a*u_x=d*u_xx
% with either periodic BCs or
% with Dirichlet BCs on left and Neumann on right:
% u(0,t)=sin(-pi*t)^exponent;
% u_x(1,t)=0 if d>0
% u(x,0)=sin(pi*x)^exponent
% If there were only advection the solution would be u(x,t)=sin(pi*(x-t))^exponent

clear
format long; format compact

global periodic;
periodic=1;

% Choose advection and diffusion coefficients
a=1;
d=0 % Advection only
%d=0.0001
%d=0.001
%d=0.01
%d=0.1
%d=1
d=0.05 % Test stability

exponent = 2; % Smooth solution
%exponent = 100; % Not so smooth solution
% Initial condition:
SOL = @(x,t) sin(pi*(x-a*t)).^exponent;
%SOL = @(x,t) sign(x-a*t-0.4)-sign(x-a*t-0.6); % Test square wave for limiting
IC = @(x) SOL(x,0);
DBC = @(t) SOL(0,t);
D = @(x) d*(2+cos(2*pi*x));
%D = @(x) d*ones(size(x)); % Constant diffusion case

% Time to compare at:
if(periodic)
   t=1.00
else
   t=0.35
end
nu=0.25 % Desired advective CFL number

base=3;
n_refinements = 5;
colors=['k','r','g','b','m','c'];

error_L1=zeros(n_refinements,1);
error_L2=zeros(n_refinements,1);
error_Linf=zeros(n_refinements,1);
h=zeros(n_refinements,1);

% Obtain finest level solution
n=2^(base+n_refinements+1)
[u_finer, x_finer, h_finer] = AdvDiff(a, D, 1.0, t, nu, n, IC, DBC);

figure(2); clf; 
figure(3); clf;
for i=n_refinements:-1:1

   n=2^(base+i)
   [u,x,h(i)] = AdvDiff(a, D, 1.0, t, nu, n, IC, DBC);
   
   if d<=0 % We know the exact solution here
      u_exact = SOL(x,t); % Because only up to second order we can just pretend this is finite difference
   else % Compare to finer grid as no exact solution known
      u_exact = Coarsen(u_finer, 2);
   end
   figure(3);
   if(i==n_refinements) plot(x, u_exact,[colors(i),'--']); hold on; end
   plot(x, u, [colors(i),'o-']); hold on; 
   
   figure(2);
   plot(x, (u-u_exact)/norm(u(1:3*n/4)-u_exact(1:3*n/4),'inf'), [colors(i),'-o']);
   hold on;

   error_L1(i)=h(i)*norm(u-u_exact,1);
   error_L2(i)=sqrt(h(i))*norm(u-u_exact,2);
   error_Linf(i)=norm(u-u_exact,'inf');
   
   u_finer=u;
   x_finer=x;
end

error_L1
error_L2
error_Linf

figure(1); clf;
loglog(h, error_L1/error_L1(1), 'ro'); hold on;
loglog(h, error_L2/error_L2(1), 'ks'); hold on;
loglog(h, error_Linf/error_Linf(1), 'gd'); hold on;
loglog(h, (h/h(1)).^2, 'r-');
loglog(h, (h/h(1)).^1, 'g-');
legend('L1','L2','Linf','2nd','1st','Location','southeast');
