% script: convectTest.m
% description: templated for numerically solving the   
% convection-diffusion equation for solution u(x,t):
%          u_t = (D*u_x)_x + P(u)_x + Q(u) + F(x,t)
% on the domain 0 <= x <= L, 0 < t <= T with boundary conditions
%         u(0,t) = B(t),
% and initial conditions
%         u(x,0) = IC(x).
% by a single step time marching method.  
% The solution is computed with user-supplied fixed time step and 
% profiles are saved at equally spaced user-defined intervals.  Treatment
% of profile error and display are defined by user after computations.
% u(0) and u(T) are always the first and last profiles recorded.
% Any global variables that might be needed for f(u,t) and/or marchIC
% are specified here. Syntax of the marching method used must be
% "marchMTD(t0, v0, tStep, numSteps)".

% programmer: t. shores
% last rev.: 3/1/16
% Octave compatible and Matlab compatible file

global gL gxNodes;			% globals used by other convect*.m files

% initialization -- set local and/or global parameters of problem
% and user must check functions specified below for fit to their problem
marchRHS = 'convectRHS';      % name of function f(u,t) to be used
marchMTD = 'ERK4';   		% name of marching method
marchIC = 'convectIC';		% name of initial condition generating function 
marchSoln = 'convectSoln';	% name of solution if known, empty string otherwise
numProfiles =  5;			% total number of computed profiles
numSteps =  25;			% total number of steps to march after each profile
T =  2     					% right endpoint in time interval [0, T] of the problem
N = 16  					% number of samples along x-axis -- must be power of 2:
gL = 2*pi;					% right endpoint in space interval [0, L] of the problem 

% done with parameter input for this test, so calculate profiles
gxNodes = linspace(0,gL,N+1); % and resulting nodes
gxNodes = gxNodes(1:N)';  % trim off right endpoint
profileIndex = [0: numSteps: numSteps*numProfiles];  % time steps at which to save a profile
tStep = T/profileIndex(end); 		  % time stepsize
uProfile = zeros(N, numProfiles + 1);   % make room for initial profile
uProfile(:,1) = feval(marchIC, gxNodes);        % and insert it
% run the method  to get remaining profiles
for jj = 2:(numProfiles + 1)
  t0 = profileIndex(jj-1)*tStep; 	% initial time for this run
  u0 = uProfile(:,jj -1); 	% initial profile for this run
  uProfile(:, jj)  = feval(marchMTD,  marchRHS, t0, u0, tStep, numSteps);
end

% This section should be rewritten by user for numeric information
% such as L2 and infinity norm of error if solution is known:
if length(marchSoln) > 0
  disp('L2 norm of error at T:');
  disp(norm(uProfile(:, numProfiles + 1) - feval(marchSoln, gxNodes, T)));
  disp('Infinity norm of error at T:');
  disp(norm(uProfile(:, numProfiles + 1) - convectSoln(gxNodes, T),inf));
end

% Some simple plotting--  initial profile has * marker and final
% profile has o marker, for plotting just use real parts
uProfile = real(uProfile);
close;
hold on;
plot(gxNodes, uProfile(:,1), '-*');
pause(1);
for jj = 2:numProfiles 
  plot(gxNodes, uProfile(:, jj));
  pause(1);
end
plot(gxNodes, uProfile(:,numProfiles+1), '-o');
legend('-* : Initial profile', '-o : Final profile')
if (length(marchSoln) > 0)
  plot(gxNodes, feval(marchSoln, gxNodes, T), 'r');
end

