% Demonstration of finite differencing of the nonlinear telegrapher's equations
function demo3
% Physical parameters
domain_width = 1; % domain size
c = 1;     % wave speed
tmax = 1  % total time for simulation

% Grid setup
Nx = 32;
Delta = domain_width/Nx;
x = Delta*(0:(Nx-1)) - domain_width/2;

% Timestep setup
SAFETY = 5 % factor by which to exceed the courant condition, h < Delta/c.
hmax = Delta/(c*SAFETY); % biggest timestep that would satisfy SAFETY.
Nt = tmax * ceil(1/hmax); % number of timesteps to reach tmax, with h < hmax,
   % ensuring that all integer values of t are on the time grid.
h = tmax/Nt; % resulting timestep, h
times = linspace(h, tmax, Nt); % timestep array, leaving out the initial t=0.
h_dimensionless = h * c / Delta % dimensionless timestep
epsilon=0.5 % big epsilon to simulate viscosity (sort of)
%epsilon = (1 - 2*h_dimensionless^2 - sqrt(1 - 4*h_dimensionless^2))...
%   / (2*h_dimensionless^2) % exact stability threshold 
%epsilon = h_dimensionless^2 + 4*h_dimensionless^4 
   % parameter required by rk2plus. Rule of thumb based on an approximation
   % of the stability diagram.

% Initial conditions
pulse_width = domain_width/8;       % std dev of Gaussian pulse.
amplitude = 0.4 % amplitude of dimensionless pulse
vel0 = amplitude*exp(-x.^2/(2*pulse_width^2)); % initial wavefunction.
den0 = vel0;     % this results in a rightward traveling pulse.
q0 = [vel0, den0];                   % entire initial state vector.

% Represent PDE in a form that my ODE solvers can understand
ode = @(q, t) nonlinear(q, t, Delta, c);

unix('rm /tmp/nonlinear*.png'); % Clear out movie frames from /tmp


% Solve the PDE using RK2 stabilized by Longcope's trick.
q = q0; % initialize q, which will become state as a function of time.
density = 1 + q( (Nx+1) : (2*Nx) ); % Total density, background + fluctuation
mass0 = sum(density)*Delta; % Initial total mass.
n=0; % zeroth movie frame, just the initial condition
movieframe(n,q,vel0,x,0,'rk2plus');
t0 = cputime;
for n = 1:Nt
   t = times(n);
   q = rk2plus(ode, q, t, h, epsilon);
   movieframe(n,q,vel0,x,t,'rk2plus');
   density = 1 + q( (Nx+1) : (2*Nx) ); % Total density, background + fluctuation
   mass_rk2plus(n) = sum(density)*Delta; % Total mass at this timepoint.
   velocity = q(1:Nx);
   KE_rk2plus(n) = sum( 0.5 .* density .* velocity.^2 );
   P_rk2plus(n) = sum( density .* velocity .^2 );
end
cputime_elapsed_min_rk2plus = (cputime - t0)/60


% Solve the PDE using RK2, saving each timestep as a movie frame.
q = q0; % initialize q, which will become state as a function of time.
n=0; % zeroth movie frame, just the initial condition
movieframe(n,q,vel0,x,0,'rk2');
t0 = cputime;
for n = 1:Nt
   t = times(n);
   q = rk2(ode, q, t, h);
   movieframe(n,q,vel0,x,t,'rk2');
   density = 1 + q( (Nx+1) : (2*Nx) ); % Total density, background + fluctuation
   mass_rk2(n) = sum(density)*Delta; % Total mass at this timepoint.
   velocity = q(1:Nx);
   KE_rk2(n) = sum( 0.5 .* density .* velocity.^2 );
   P_rk2(n) = sum( density .* velocity .^2 );
end
cputime_elapsed_min_rk2 = (cputime - t0)/60

mass_change_rk2plus = max(mass_rk2plus) - min(mass_rk2plus)

% *** Graphics ***

% Test of mass conservation:
figure(4)
hold off
plot(times, (mass_rk2-mass0)/mass0, 'r-', 'linewidth', 2)
hold on
plot(times, (mass_rk2plus-mass0)/mass0, 'b-', 'linewidth', 2)
title('Mass Conservation in Fluid Equations')
xlabel('t')
ylabel('\Delta m / m_0')
legend('rk2','rk2plus')
hold off
print('mass_conservation.pdf', '-dpdf')

% Test of energy conservation (should this be conserved?):
figure(5)
hold off
plot(times, KE_rk2, 'r-', 'linewidth', 2)
hold on
plot(times, KE_rk2plus, 'b-', 'linewidth', 2)
title('Kinetic Energy in Fluid Equations')
xlabel('t')
ylabel('Kinetic Energy')
legend('rk2','rk2plus','location','northwest')
hold off
print('kinetic_energy.pdf', '-dpdf')

% Test of momentum conservation (should this be conserved?):
figure(6)
hold off
plot(times, P_rk2, 'r-', 'linewidth', 2)
hold on
plot(times, P_rk2plus, 'b-', 'linewidth', 2)
title('Momentum in Fluid Equations')
xlabel('t')
ylabel('Momentum')
legend('rk2','rk2plus','location','northwest')
hold off
print('momentum.pdf', '-dpdf')

function movieframe(n,q,u0,x,t,tag)
% Routine for creating a single movie frame with
% frame number n,
% state q = [u, udot],
% initial wave function u0,
% position grid (N-vector) x,
% time t,
% and short description given by tag.
myfont = '/sw/lib/X11/fonts/applettf/Arial.ttf'; 
   % Name of a font file that works for print() to make .png images. 
Nx = numel(x);
u = q(1:Nx);
figure(1)
hold off
plot(x,u,'b-', x,u0,'r-')
hold on
legend( [tag,' solution'], 'initial condition' );
axis([x(1),x(Nx),-0.2,0.6]);
xlabel('x','FontName',myfont)
ylabel('dimensionless velocity v/c','FontName',myfont)
title( sprintf('%s, frame %03i, t = %0.3f',tag,n,t),'FontName', myfont );
print( sprintf('/tmp/nonlinear_%s_%03i.png',tag,n) ); % save movie frame

