function retval = marchTest()
% usage: retval = marchTest()
% description: templated for numerically solving the equation
% (convection-diffusion if D != 0) equation for solution u(x,t):
%          u_t = (D(x)*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, whose syntax must be
% "marchMTD(t0, v0, tStep, numSteps)".
% 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(x, 0) and u(x, T) are always the first and last profiles recorded.

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

% globals used by other function files
global gL gxNodes gpFlag gqFlag gfFlag gbFlag;

% initialization -- set local and/or global parameters of problem
% and user must check functions specified below for fit to their problem
% user set choice of using extra terms or not here
gpFlag =false;     % true if P(x,t,u), defined by user as marchP in file marchRHS.m, is nonzero
gqFlag = false;    % true if Q(x,t,u), defined by user as marchQ in file marchRHS.m, is nonzero
gfFlag = false;      % true if F(x,t), defined by user as marchF in file marchRHS.m, is nonzero
gbFlag = true;     % true if left BC B(t), defined by user as marchB in file marchRHS, is to be met
uFlag = true;        % true if solution u(x,t),  defined by user below in marchU, is known
marchMTD = 'ERK4';   		% name of marching method
numProfiles =  5;			% total number of computed profiles
numSteps =  50;			% total number of steps to march after each profile
T =  2     					% right endpoint in time interval [0, T] of the problem
N = 32  					% 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 -- NB: it is user's responsibility to pay attention to
					   % suitability, like a CFL condition dt <= dx^2/(2*Dmax)
uProfile = zeros(N, numProfiles + 1);   % make room for initial profile
uProfile(:,1) = marchIC(gxNodes);        % insert it
uProfile(:, 1) = FFTransform(uProfile(:, 1), 0); %  and FFT it
%uProfile(:, 1) = fft(uProfile(:, 1)); %  if you need more speed
% run the method in frequency space 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
% dealias fourier profiles and return results to physical space (except known initial)
uProfile(:,1) = marchIC(gxNodes);
for jj = 2:(numProfiles + 1)
  uProfile(ceil(N/3):(N - floor(N/3)+1), jj) = 0;
  uProfile(:, jj) = FFTransform(uProfile(:, jj), 1);
  %uProfile(:, jj) = ifft(uProfile(:, jj)); %  if you need more speed, but be consistent
end

% This section should be rewritten by user for numeric information
% such as L2 and infinity norm of error if solution is known:
if (uFlag)
  disp('Absolute L2 norm of error at T:');
  disp(norm(uProfile(:, numProfiles + 1) - marchSoln(gxNodes, T)));
   disp('Relative L2 norm of error at T:');
  disp(norm(uProfile(:, numProfiles + 1) - marchSoln(gxNodes, T))/norm(marchSoln(gxNodes, T)));
 disp('Absolute infinity norm of error at T:');
  disp(norm(uProfile(:, numProfiles + 1) - marchSoln(gxNodes, T),inf));
 disp('Relative infinity norm of error at T:');
  disp(norm(uProfile(:, numProfiles + 1) - marchSoln(gxNodes, T),inf)/norm(marchSoln(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', '-or: Final profile')
% if exact solution is known, plot it in red
if (uFlag)
  plot(gxNodes, marchSoln(gxNodes, T), 'r');
end
retval = 0
end    % marchTest

function retval = marchIC(x)
% usage: retval = marchIC(x)
% description: This function computes the initial condition
% for convection-diffusion equation at time t = 0.  This
% function should match the left boundary condition
% marchB(0) at x = 0 and be periodic of period gL.
 
retval = sin(x);
%retval = sin(x).^2;
%retval = zeros(size(x)) + ((x >= 2) & (x <= 3));
%retval = zeros(size(x));
%retval = exp(-(x-3).^2/0.4)/sqrt(0.2*2*pi);
end    % marchIC

function retval = marchSoln(x,t)
% usage: retval = marchSolnf(x,t)
% description: This function accepts input vector x, scalar
% time t, and computes the supplied solution for the
% convection-diffusion problem.

% default is constant value zero, handy if soln not known and size of last profile desired
%retval = zeros(size(x));
retval = exp(-t)*sin(x);
%retval = exp(-t/2)*sin(x).^2;
%retval = t*sin(x);
%retval = t*cos(x);
%retval = t*x.*(2*pi - x);
%retval = t*x.*x.*(2*pi - x);
end    % marchSoln

