function retval = marchTest8()
% usage: retval = marchTest8()
% 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/29/16
% Octave compatible and Matlab compatible file

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

% 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 = true;      % 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 = false;     % 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 = 'ERK23';   		% name of marching method
numProfiles =  6;			% total number of computed profiles
numSteps =  10;			% total number of steps to march after each profile
T = 6;          				% right endpoint in time interval [0, T] of the problem
N = 32;  					% number of samples along x-axis -- must be power of 2:
gL = pi;					% right endpoint in space interval [0, L] of the problem 

% done with parameter input for this test, so calculate profiles
gxCount = 0;	% this is the count for rhs evaluations
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,  'marchRHS8', 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

% Delete the following line if complex solutions are relevant
uProfile = real(uProfile);

% This section should be rewritten by user for numeric information such as
% L2 and infinity norm of error if solution is known, function evaluation in any case:
sString = sprintf('Final time for this run is T = %g.', T);
disp(sString);
sString = sprintf('Sample number along x-axis for this run is N = %d.', N);
disp(sString);
sString = sprintf('Step size in time for this run is tStep = %g.', tStep);
disp(sString);
sString = sprintf('Function evaluations for this run with method %s is %g.', marchMTD, gxCount);
disp(sString);
if (uFlag)
  % account for errors at profile time T
   solnVec = marchU(gxNodes, T);
   errVec = uProfile(:, numProfiles + 1) - solnVec;
   % vTmp0 counts relative error at points where solnVec = 0 as 0
   vTmp0 = errVec;
   vTmp0(find(solnVec == 0)) = 0;
   % vTmp1 counts solution at points where solnVec = 0 as 1
   vTmp1 = solnVec;
   vTmp1(find(solnVec == 0)) = 1;
   relErrVec = vTmp0./vTmp1;
   sdX = sqrt(gL/numProfiles);
   sString = sprintf('L2 norm of absolute errors at time T is %g.', norm(errVec)*sdX);
   disp(sString);
   sString = sprintf('L2 norm of error by L2 norm of solution at time T is %g.', norm(errVec)/norm(solnVec));
   disp(sString);  
   sString = sprintf('Infinity norm of absolute errors at time T is %g.', norm(errVec,inf));
   disp(sString);
   sString = sprintf('Largest absolute error occurs at x = %g.', gxNodes(find(abs(errVec) == norm(errVec,inf))));
   disp(sString);
   sString =   sprintf('Infinity norm of relative errors at time T is %g.', norm(relErrVec,inf));
   disp(sString);
   sString = sprintf('Largest relative error occurs at x = %g.\n', gxNodes(find(abs(relErrVec) == norm(relErrVec,inf))));
   disp(sString);
end

% Some simple plotting--  initial profile has dashes and final
% profile has o marker, for plotting just use real parts
uProfile = real(uProfile);
close;
hold on;
pause(1);
plot(gxNodes, uProfile(:,1), '--');
pause(1);
for jj = 2:(numProfiles + 1)
  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, marchU(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);
end    % marchIC

function retval = marchU(x,t)
% usage: retval = marchU(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 =sin(x).*exp(t*cos(x));
end    % marchSoln

