% 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
method=1 % centered difference
%method=2 % third-order upwind biased

times=[0, 0.25, 0.35, 0.75]; % What times to output the solution at

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

% Choose desired resolution:
%n=32;
%n=64;
%n=128;
n=256;

h=1/n; % Grid spacing
x=([1:n]-0.5)*h; % Cell centers
IC = @(x) sin(pi*x).^exponent;
u0 = FVQuad(IC,x,h);

[T,W] = SolveODE(n, u0, times, h, a, d, method);

% Compute an accurate solution for comparison:
if d<=0 % We know the exact solution here
   W_exact=zeros(3,n);
   for i=1:length(times)
      sol = @(x) sin(pi*(x-times(i))).^exponent;
      W_exact(i,:)=FVQuad(sol,x,h);
   end
   u_exact=W_exact;
else % Just use a refined computation with the best method to get the 'exact' solution
   n_max=1024;
   h_min=1/n_max;
   w0 = FVQuad(IC, ([1:n_max]-0.5)*h_min, h_min);
   [T_exact,W_exact] = SolveODE(n_max, w0, times, h_min, a, d, 3);
   for i=1:length(times) % Must coarsen "exact" solution to compare properly
      u_exact(i,:)=Coarsen(W_exact(i,:), n_max/n);
   end
end

% Plot solution:
figure(1); clf

plot(x, u0, 'k-'); hold on; % Exact/numerical solution at t=0

plot(x, W(2,:), 'ro'); hold on; % Numerical solution
plot(x, W(3,:), 'gs'); hold on; % Numerical solution
plot(x, W(4,:), 'bd'); hold on; % Numerical solution

plot(x, u_exact(2,:), 'r-'); hold on; % Exact solution
plot(x, u_exact(3,:), 'g--'); hold on; % Exact solution
plot(x, u_exact(4,:), 'b-.'); hold on; % Exact solution

legend('t=0','t=0.25','t=0.5','t=0.75','exact','exact','exact');

% Plot errors
figure(2); clf

plot(x, W(2,:)-u_exact(2,:), 'r-o'); hold on; % Numerical solution
plot(x, W(3,:)-u_exact(3,:), 'g--s'); hold on; % Numerical solution
plot(x, W(4,:)-u_exact(4,:), 'b-.d'); hold on; % Numerical solution

legend('t=0.25','t=0.5','t=0.75');
