% Solve the Burger's equation problem:;
% u_t + u u_x = 0     u(x,0) = 1 + sin(x)    u(x,t) = u(x + 2 pi,t);

% This program uses a nonconservative upwind scheme

% Put in the same directory with file burgex.m, run Matlab, enter: burgnc

% Plot results using:
%     ~dommelen/bin/xyplot burgnc.plt 2 ' ' ' ' ' ' 2 '0,100,0,100,.4'
% which creates a printable postscript file named dump.ps

% To use your own plotting software, change the format of file burgnc.plt
% below as needed.

format compact
format short

% since in matlab, arrays must have positive indices, define index offsets
offj=2;
offn=1;

% open files

% input file
infile=fopen('burg.in','r');

% output file
outfile=fopen('burgnc.out','w');

% plot file
plotfile=fopen('burgnc.plt','w');

% read input file and write to output

u0=fscanf(infile,'%f');
comment=fgetl(infile);
fprintf(outfile,'average initial velocity: %f\n',u0);

j_max=fscanf(infile,'%d');
comment=fgetl(infile);
fprintf(outfile,'number of mesh cells j_max: %d\n',j_max);
if j_max<2,
  error('*** number of mesh cells should be at least 2');
end

C=fscanf(infile,'%f');
comment=fgetl(infile);
fprintf(outfile,'maximum Courant number u Delta t/Delta x: %f\n',C);
if C<=0.,
  error('*** Courant number must be positive');
end

t_out=fscanf(infile,'%f');
comment=fgetl(infile);
fprintf(outfile,'output time interval: %f\n',t_out);
if t_out<=0.,
  error('*** output time interval must be positive');
end

outmx=fscanf(infile,'%d');
comment=fgetl(infile);
fprintf(outfile,'number of output time intervals to compute: %f\n',outmx);
if outmx<=0,
  error('*** number of output times must be positive');
end

fclose(infile);

% initialize the parameters

% mesh size
dx=2.*pi/j_max;

% output count
out=1;

% print value of n
nprint=1;

% initial velocities

% set all initial velocities and x-positions (starting j from -1 instead of
% 0 makes it easier to handle periodicity later)
jvals=[-1:j_max]';
x = dx*jvals;
u = [u0+sin(x), u0+sin(x)];

% initialize time
t=0.;

% main loop over all time steps

for n=0:123456789,

  % where to store the velocities at n and n+1
  n_s=mod(n,2);
  np1_s=mod(n+1,2);

  % perform output
  while t>=out*t_out,

    % write computed velocities to the plot file
    fprintf(plotfile,'%d\n',-j_max-1);
    for j=0:j_max,
      fprintf(plotfile,'%f,%f\n',x(offj+j),u(offj+j,offn+n_s));
    end

    % write the exact solution to the plot file as a 201 point curve
    fprintf(plotfile,'%d\n',201);
    for j=0:200,
      fprintf(plotfile,'%f,%f\n',0.005*2.*pi*j,burgex(0.005*2.*pi*j,t,u0));
    end

    out=out+1;
  end

  % stop if no more output is needed
  if out>outmx, break; end

  % compute the time step
  dt=1.e30;
  for j=0:j_max-1,
    if u(offj+j,offn+n_s)~=0., dt=min(dt,C*dx/abs(u(offj+j,offn+n_s))); end
  end

  % increment time
  t=t+dt;

  % set the j=-1 value of u using periodicity
  u(offj-1,offn+n_s)=u(offj+j_max-1,offn+n_s);

  % compute the velocities
  for j=0:j_max-1,
    if u(offj+j,offn+n_s)>=0.,
      u(offj+j,offn+np1_s)=u(offj+j,offn+n_s) ...
        -dt/dx*u(offj+j,offn+n_s)*(u(offj+j,offn+n_s)-u(offj+j-1,offn+n_s));
    else
      u(offj+j,offn+np1_s)=u(offj+j,offn+n_s) ...
        -dt/dx*u(offj+j,offn+n_s)*(u(offj+j+1,offn+n_s)-u(offj+j,offn+n_s));
    end
  end

  % set the boundary velocity at j_max using periodicity
  u(offj+j_max,offn+np1_s)=u(offj+0,offn+np1_s);

  % print diagnostics only occasionally
  if n+1~=nprint, continue; end
  nprint=nprint*2;

  % compute the rms error and minimum and maximum velocity
  errrms=0.;
  u_min=u(offj+0,offn+np1_s);
  u_max=u(offj+0,offn+np1_s);
  for j=0:j_max-1,
    err=abs(u(offj+j,offn+np1_s)-burgex(x(offj+j),t,u0));
    errrms=errrms+err^2;
    u_min=min(u_min,u(offj+j,offn+np1_s));
    u_max=max(u_max,u(offj+j,offn+np1_s));
  end
  errrms=sqrt(errrms/j_max);
  fprintf('n = %d   t = %f  u-range: %f-%f   error: %f\n', ...
   n+1,t,u_min,u_max,errrms);
  fprintf(outfile,'n = %d   t = %f  u-range: %f-%f   error: %f\n', ...
   n+1,t,u_min,u_max,errrms);

end

% also show the results on the screen (in addition to plot file burgnc.plt)

% finite difference results
plotx=x(offj+0:offj+j_max);
ploty=u(offj+0:offj+j_max,offn+n_s);

% exact results
dxe=2.*pi/50.;
plotj=[0:50]';
plotxe=dxe*plotj;
plotye=plotxe;
for j=1:51,
  plotye(j)=burgex(plotxe(j),t,u0);
end

% plot finite difference as points, exact as a line
plot(plotx,ploty,'o',plotxe,plotye,'-');
