%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                            %
% program for computation of Euler Riemann problems          %
% with conservative mach-uniform pressure-correction         %
% method.                                                    %
%                                                            %
% The code is used to generate figures 4, 7 and 10           %
% of :                                                       %
%                                                            %
%  D.R. van der Heul, C. Vuik and P. Wesseling               %
% 'A Conservative Pressure Correction Method for             %
%  Compressible Flows at All Speeds'                         %
%  Computers and Fluids, 32, 1113-1132, 2003                 %
%  http://ta.twi.tudelft.nl/nw/users/vuik/papers/Heu03VW.pdf %
%                                                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Variables
% mx    = momentum field in x direction
% ux    = velocity field in x direction
% rho   = density
% p     = pressure
% gamma = ratio of specific heat constants

global J lambda
global  PRL  CRL MACHLEFT  gamma  pleft  pright  rholeft  rhoright  uleft...
	uright  tend  lambda
J =  400;                       % number of gridcells
h =   1/J;                      % meshwidth
xwall=-0.5-h+(1:J+1)*h;         % location of momentum 
xcenter=-0.5-h/2+(1:J)*h;       % location of scalars

% initialize the variables

mxold  =   zeros(J+1,1);
uxold  =   zeros(J+1,1);
rhoold =   zeros(J  ,1);
pold   =   zeros(J  ,1);
mxnew  =   mxold       ;
uxnew  =   uxold       ;
rhonew =   rhoold      ;
pnew   =   pold        ;

% initialize help variables

rhocent =  zeros(J+1,1);
uxc     =  zeros(J  ,1);

% initialize matrices 
P       =  zeros(J  ,3);
A       =  zeros(J  ,3);
B	=  zeros(J+1,3);

% define boundary conditions

gamma=1.4;
EULBC;

% fill with initial condition

EULInit;

% settings for timeintegration

dt     = lambda*h;
n      = floor(tend/dt);

% copy new solution to old timelevel

EULShift;

%-----------------------------

% start of time loop

for i=1:n
  
disp(i);

% Build matrix for scalar equation

  EULscalmat;

% Build matrix for momentum equation x-direction

  EULMmat;
     
% compute additional vectors for timeintegration
     
  EULrhof;
  EULMxf;

% advance solution 

  EULstage;
     
% compute additional quantities
     
  EULadst;

% build pressure correction matrix and solve
% for pressure correction
  
  EULDPP;

% update momentum field and pressure

  EULUpdate;


% shift solution

  EULShift;

end

% end of time loop

%--------------------------------------

% redimensionalize all quantities

  EULDim;

% compute additional quantities

  EULAD;
    
% plot solution

  EULplot;
