%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  PHYS 210  Fall 2009
%
%  Solution of 1D advection equation as discussed in class.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Turn off the paging of standard output so that messages appear
% immediately on the screen instead of being buffered.
more off

% Set parameters for a level 6 solution

% Note that tmax is relatively small so that pulse doesn't have 
% enough time to hit the right boundary, after which point we 
% don't expect the convergence of the solution to be O(h^2)
tmax = 0.5;
level = 6;
lambda = 0.5;
al = 1.0;
x0l = 0.75;
deltal = 0.05;

fprintf('Running advect at level %d\n',level);
[x t u] = advect(tmax, level, lambda, al, x0l, deltal);

% Make surface plot and display it 
fprintf('Generating surface plot of soln\n');
clf;
hold on;
xlabel('{\bf x}');
ylabel('{\bf t}');
zlabel('{\bf u}');
view(39,47);
mesh(x, t, u);
input('Type ''Enter'' to continue: ');

% Check convergence of code using 4 discretization levels.
fprintf('Checking code convergence\n');
fprintf('Running advect at level 6\n');
[x6 t6 u6] = advect(tmax, 6, lambda, al, x0l, deltal);
fprintf('Running advect at level 7\n');
[x7 t7 u7] = advect(tmax, 7, lambda, al, x0l, deltal);
fprintf('Running advect at level 8\n');
[x8 t8 u8] = advect(tmax, 8, lambda, al, x0l, deltal);
fprintf('Running advect at level 9\n');
[x9 t9 u9] = advect(tmax, 9, lambda, al, x0l, deltal);
fprintf('Running advect at level 10\n');
[x10 t10 u10] = advect(tmax, 10, lambda, al, x0l, deltal);

% Compute differences between solution values on adjacent levels
% and using common set of mesh points.
[nt6, nx6] = size(u6);
[nt7, nx7] = size(u7);
[nt8, nx8] = size(u8);
[nt9, nx9] = size(u9); 
[nt10, nx10] = size(u10); 
du67 = u6 - u7(1:2:nt7,1:2:nx7);
du78 = u7 - u8(1:2:nt8,1:2:nx8);
du89 = u8 - u9(1:2:nt9,1:2:nx9);
du910 = u9 - u10(1:2:nt10,1:2:nx10);

% Run on levels 9 and 10 for tmax=10 to illustrate effect of 
% boundary condition.  Note the spurious reflection of the "outgoing"
% pulse, which is evident from the surface plot, and completely 
% absent from the continuum/exact solution.   Also note that the 
% amplitude of the reflected pulse decreases as the mesh spacing 
% decreases. 
fprintf('Will now run until tmax = 2.5 to illustrate boundary effects\n');
tmax = 2.5
fprintf('Running advect at level 9\n');
[x9 t9 u9] = advect(tmax, 9, lambda, al, x0l, deltal);
fprintf('Running advect at level 10\n');
[x10 t10 u10] = advect(tmax, 10, lambda, al, x0l, deltal);

[nt9, nx9] = size(u9); 
[nt10, nx10] = size(u10); 

fprintf('Generating level 9 surface plot of soln\n');
clf;
hold on;
xlabel('{\bf x}');
ylabel('{\bf t}');
zlabel('{\bf u_9}');
view(39,47);
x9 = x9(1:8:nx9);
t9 = t9(1:32:nt9);
u9 = u9(1:32:nt9, 1:8:nx9);
mesh(x9, t9, u9);
input('Type ''Enter'' to continue: ');

fprintf('Generating level 10 surface plot of soln\n');
clf;
hold on;
xlabel('{\bf x}');
ylabel('{\bf t}');
zlabel('{\bf u_{10}}');
view(39,47);
x10 = x10(1:16:nx10);
t10 = t10(1:64:nt10);
u10 = u10(1:64:nt10, 1:16:nx10);
mesh(x10, t10, u10);
input('Type ''Enter'' to continue: ');

fprintf('Generating surface plot of u_9 - u_10\n');
clf;
hold on;
xlabel('{\bf x}');
ylabel('{\bf t}');
zlabel('{\bf u_9 - u_{10}}');
view(39,47);
du910 = 4 * (u9 - u10);
mesh(x9, t9, du910);
input('Type ''Enter'' to continue: ');

% Compare run times of loop-less and loop-ful versions of 
% implementation.  
% 
% Use 'tic' and 'toc' functions (MATLAB/octave) which start 
% and stop a "wall-clock" timer: i.e. NOT a cpu timer, but 
% a real-time timer, so we want to ensure that the calculations
% being benchmarked take a reasonable amount of time, so run 
% with level = 9, tmax = 0.5

fprintf('Running timed run of advect at level 9\n');
level = 9
tmax = 0.5
input('Type ''Enter'' to continue: ');
tic
[x9 t9 u9] = advect(tmax, level, lambda, al, x0l, deltal);
toc

fprintf('Running timed run of advectloop at level 9\n');
level = 9
tmax = 0.5
input('Type ''Enter'' to continue: ');
tic
[x9 t9 u9] = advectloop(tmax, level, lambda, al, x0l, deltal);
toc

% Finally, run with lambda > 1 to illustrate the resulting 
% (numerical) instability.  The "saw tooth" nature of the blow-up
% is characteristic of the instability that one generically sees 
% when the Courant number exceeds some limit (in this case 1).
% Such a limit is often referred to as the CFL limit, after the 
% mathematicians Courant, Friedrichs and Levy, who studied finite
% difference approximations in the 1920's as a means to establish
% existence and uniqueness of solutions to the wave equations.
% They identified the analogous limit for the wave equation
% as being a necessary condition for convergence of the difference
% approximation, and a general rule of thumb (which is a theorem
% in many cases) for finite-difference discretizations of 
% time-dependent PDEs is that stability and  convergence are 
% equivalent: i.e. a FDA converges if and only if it is stable, 
% and a FDA is stable if and only if it converges.

% Note that the value of lambda adopted below was carefully chosen
% to make the instability manifest, while still allowing you 
% to see the "remnants" of the continuum solution.  Moreover, 
% for CFL-violating instabilities, as the mesh is refined, the 
% solutions tends to blow up even faster with respect to 
% time.
fprintf('Will now run with lambda > 1 to illustrate instability\n');
input('Type ''Enter'' to continue: ');

tmax = 0.5;
lambda = 1.05;

fprintf('Running advect at level 6\n');
[x6 t6 u6] = advect(tmax, 6, lambda, al, x0l, deltal);

fprintf('Generating surface plot of soln\n');
clf;
hold on;
xlabel('{\bf x}');
ylabel('{\bf t}');
zlabel('{\bf u_6}');
view(39,47);
mesh(x6, t6, u6);
input('Type ''Enter'' to continue: ');

fprintf('\nDone!!\n\n');
