% Finite volume spatial discretization to solve:
% u_t+a*u_x=d*u_xx
% 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 short; format compact

global exponent;
exponent = 2; % Smooth solution
%exponent = 100; % Not so smooth solution

%method=0 % upwind (first order accurate for any d)
%method=1 % centered difference (2nd-order for d=0 or large d, 1st order in L_inf for small d)
method=2 % third-order upwind biased (3rd order for d=0, 2nd for large d, problems for small d)

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

% Time to compare at:
t=0.35;

base=2
n_refinements = 4
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] = SolvePDE(n, t, a, d, method);

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

   n=2^(base+i)
   [u,x,h(i)] = SolvePDE(n, t, a, d, method);
   
   if d<=0 % We know the exact solution here
      sol = @(x) sin(pi*(x-t)).^exponent;
      u_exact = FVQuad(sol,x,h(i));
   else % Compare to finer grid as no exact solution known
      u_exact = Coarsen(u_finer, 2)';
      %figure(3); clf; plot(x, u_exact,'k--s'); hold on; plot(x, u, 'g-d');
   end
   
   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

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)).^3, 'k-');
loglog(h, (h/h(1)).^2, 'r-');
loglog(h, (h/h(1)).^1, 'g-');
legend('L1','L2','Linf','3rd','2nd','1st','Location','southeast');
