%------------------------------------------------------------------
%
%     1-D TIME-MARCHING EULER SOLVER FOR VARIABLE DUCT AREA.
%     BASED ON CELL-CENTERED UPWIND-BIASED FINITE-VOLUME SCHEME
%     AND EXPLICIT RUNGE-KUTTA TIME-INTEGRATION.
%
%
%     WRITTEN BY:      Lars-Erik Eriksson
%                      Chalmers University of Technology
%                      May 2006
%
%     TRANSLATED TO MATLAB BY: 
%                      Martin Olausson
%                      Chalmers University of Technology
%                      Mars 2007
%
%     rev1: added cfl-based timestep, 2009-01-14
%     rev2: minor changes to speed up plotting, 2015-12-07 /Marcus Lejon
%
%------------------------------------------------------------------

clear
close all

disp(' ')
disp('**************')
disp('* G1DEUL     *')
disp('**************')

%--------------------------------------------
%     Define grid nodes and corresponding
%     cross section areas and circumferential 
%     lengths. Store in XM,SX,SCFL arrays. 
%     NI = total number of grid nodes.
%--------------------------------------------
r  = 5e-3;                          % Radius of tube
X1 = 0.0;                           % Left end of mesh
X2 = 1.0;                           % Right en of mesh
MSH.NI = 101;                       % Number of mesh points
MSH.XM = linspace(X1,X2,MSH.NI)';   % Mesh point vector

% Assume circular pipe
MSH.SX = r^2*pi*ones(MSH.NI,1);     % Cross section area vector
MSH.SCFL = sqrt(4*MSH.SX*pi);       % Circumferential lengths vector

% Plot SX
figure(2)
plot(MSH.XM,MSH.SX,'k--')
legend('SX (area)')
xlabel('XM')

%-------------------------------------------
%     Compute grid cell volumes and
%     store inverse of these in VOLI array.
%-------------------------------------------
MSH.VOLI = mtr(MSH);                % Inverse of cell volume


%-------------------------------------------
%     Compute grid cell centers
%-------------------------------------------
MSH.XCP = .5 * ...
    (MSH.XM(1:MSH.NI-1,1) + ...
     MSH.XM(2:MSH.NI,1));

 
%-------------------------------------------
%     Define gas properties; GCNT=gamma=cp/cv
%                            RGAS=Runiv/M
%     For air at standard conditions:
%     GCNT=1.4, RGAS=287   (SI units)
%-------------------------------------------
GAS.GCNT = 1.4;                    % gamma
GAS.RGAS = 287;                    % R


%-------------------------------------------
%     Define friction loss coefficient
%     and wall heat transfer:
%     FFRIC = f (dimensionless)
%             (typical value ~0.005)
%     QWALL = heat flux/unit area
%             (W/m2 in SI units)
%-------------------------------------------

DIFF.FFRIC = 0.;                   % Friction loss coeff
DIFF.QWALL = 213934.0;             % Heat flux/unit area


%--------------------------------------------
%     Define initial solution and store in
%     RO,RU,ET arrays.
%--------------------------------------------

%--------------------------------------------------------------
%     Create constant initial solution.
%--------------------------------------------------------------
ROINIT=0.214;
UINIT =400;
PINIT =0.73e5;

SOLV.RO = ones(MSH.NI-1,1)*ROINIT;
SOLV.RU = ones(MSH.NI-1,1)*ROINIT*UINIT;
SOLV.ET = ones(MSH.NI-1,1)*PINIT/(GAS.GCNT-1.)+.5*ROINIT*UINIT^2;


%--------------------------------------------------------------
% Compute auxiliary flow variables
%--------------------------------------------------------------
disp('auxvr:')
SOLV = auxvr(MSH,GAS,SOLV,1);

SOLV.TIME = 0.0;      % Set start time to zero


%-------------------------------------------
%     Define boundary conditions.
%     'LB' stands for 'Left Boundary'
%     'RB' stands for 'Right Boundary'
%     IT defines type of b.c.:
%       IT=0  gives non-reflecting b.c.
%             RO,U,P are given
%       IT=1  gives subsonic inflow b.c.
%             T0,P0 are given
%       IT=2  gives subsonic outflow b.c.
%             P is given
%       IT=3  gives supersonic inflow b.c.
%             RO,U,P are given
%       IT=4  gives supersonic outflow b.c.
%             No data is given
%       IT=5  gives solid wall b.c.
%             No data is given
%       IT=10 gives periodic b.c.
%             No data is given
%-------------------------------------------

BNDC.ITLB = 1;  % Type of boundary condition, left side (IT=1; subsonic inflow)
BNDC.ROLB = nan;            % Density   (not used if IT=1)
BNDC.ULB  = nan;            % Velocity  (not used if IT=1)
BNDC.PLB  = nan;            % Pressure  (not used if IT=1)
BNDC.T0LB = 275.2;          % Total temperature
BNDC.P0LB = 104136.4;       % Total pressure

BNDC.ITRB = 2;  % Type of boundary condition, right side (IT=2; subsonic outflow)
BNDC.RORB = nan;            % Density   (not used if IT=2)
BNDC.URB  = nan;            % Velocity  (not used if IT=2)
BNDC.PRB  = 72733.4;        % Pressure  
BNDC.T0RB = nan;            % Total temperature  (not used if IT=2)
BNDC.P0RB = nan;            % Total pressure     (not used if IT=2)


%---------------------------------------------------
%     Define coefficients for numerical flux scheme.
%       Standard 3rd-order upwind scheme:
%         -1/6, 5/6, 2/6, 0
%       Centered 2nd-order scheme:
%            0, 1/2, 1/2, 0
%       Centered 4th-order scheme:
%         -1/12, 7/12, 7/12, -1/12
%       First-order upwind scheme:
%            0, 1, 0, 0
%---------------------------------------------------
COEF.COEF1 = -1/6;
COEF.COEF2 =  5/6;
COEF.COEF3 =  2/6;
COEF.COEF4 =    0;


%----------------------------------------------------
%     Define coefficients for extra pressure-
%     and density-controlled numerical dissipation.
%     For smooth solution:  0, 0
%     For compression shocks:   0.5, 0
%     For shocks and strong contacts:  0.5, 2
%----------------------------------------------------
DAMP.FCNTP = 0;
DAMP.FCNTC = 0;


%-----------------------------
%     Define time step and
%     number of time steps.
%-----------------------------
%TSTP.DT    = 5.e-6;      % Time step size (set by CFL-number)
TSTP.CFL   = 0.9;        % Maximum CFL-number
TSTP.NSTEP = 10000;      % Number of time steps


%---------------------------------
% STOP or PAUSE iterating on 
% users request
%---------------------------------
BUTTON.doSTOP         =0;
BUTTON.doPAUSE        =0;
BUTTON.changePLOTSTEP =0;
BUTTON.PLOTSTEP       =20;
BUTTON.PLOTSTEPINC    =10;
figure(1)
BUTTON.hSTOP         = uicontrol('String','STOP' ,                          'Callback','BUTTON.doSTOP =1;',                                                 'Position', [5   5 60  20]);
BUTTON.hPAUSE        = uicontrol('String','PAUSE',                          'Callback','BUTTON.doPAUSE=1;',                                                 'Position', [70  5 60  20]);
BUTTON.hdecPLOTSTEP  = uicontrol('String',['-' num2str(BUTTON.PLOTSTEPINC)],'Callback','BUTTON.changePLOTSTEP=BUTTON.changePLOTSTEP-BUTTON.PLOTSTEPINC;',   'Position', [150 5 30  20]);
BUTTON.hincPLOTSTEP  = uicontrol('String',['+' num2str(BUTTON.PLOTSTEPINC)],'Callback','BUTTON.changePLOTSTEP=BUTTON.changePLOTSTEP+BUTTON.PLOTSTEPINC;',   'Position', [185 5 30  20]);
BUTTON.htextPLOTSTEP = uicontrol('Style','Text','String',['Plot every ' num2str(BUTTON.PLOTSTEP) ' timestep'],                                              'Position', [220 5 200 15]);


%------------------------------
%     MAIN LOOP
%------------------------------

ROREF=SOLV.RO;

for ISTEP=1:TSTP.NSTEP
    
    % Calculate timestep
    TSTP.DT = TSTP.CFL/max((SOLV.C+abs(SOLV.U))./(MSH.XM(2:MSH.NI,1)-MSH.XM(1:MSH.NI-1,1)));
    
    % Step, calculate solution at next timestep
    SOLV = stepexpl(MSH,GAS,DIFF,SOLV,BNDC,COEF,DAMP,TSTP);

    % Calculate and print maximum change in density
    [RMAX(ISTEP),IMAX]=max(abs(SOLV.RO - ROREF));
    ROREF=SOLV.RO;
    cflmax=max(TSTP.DT*(SOLV.C+abs(SOLV.U))./(MSH.XM(2:MSH.NI,1)-MSH.XM(1:MSH.NI-1,1)));
    disp(['ISTEP=' num2str(ISTEP) ...
        ', RMAX=' num2str(RMAX(ISTEP)) ...
        ', IMAX=' num2str(IMAX) ...
        ', TSTP.DT=' num2str(TSTP.DT) ...
        ', CFLMAX=' num2str(cflmax)]);
    
    % Plot solution during iterations
    if mod(ISTEP,BUTTON.PLOTSTEP) == 0 || ISTEP == TSTP.NSTEP
        % Calculate additional variables to plot
        T=SOLV.P./(SOLV.RO*GAS.RGAS);       % Static temperature
        M=abs(SOLV.U)./SOLV.C;              % Mach number
        T0=T.*(1+0.5*(GAS.GCNT-1)*M.^2);    % Total temperature
        P0=SOLV.P.*(1+0.5*(GAS.GCNT-1)*M.^2).^(GAS.GCNT/(GAS.GCNT-1)); % Total pressure
        figure(1)
        % density
        subplot(331)
        plot(MSH.XCP,SOLV.RO,'b')
        title('RO (density)')
        xlabel('X')
        % velocity
        subplot(332)
        plot(MSH.XCP,SOLV.U,'g')
        title('U (velocity)')
        xlabel('X')
        % speed of sound
        subplot(333)
        plot(MSH.XCP,SOLV.C,'c')
        title('C (speed of sound)')
        xlabel('X')
        % pressure
        subplot(334)
        plot(MSH.XCP,SOLV.P,'r')
        title('P (pressure)')
        xlabel('X')
        % temperature
        subplot(335)
        plot(MSH.XCP,T,'m')
        title('T (temperature)')
        xlabel('X')
        % mach
        subplot(336)
        plot(MSH.XCP,M,'k')
        title('M (Mach nr)')
        xlabel('X')
        % P0
        subplot(337)
        plot(MSH.XCP,P0,'r--')
        title('P0 (total pressure)')
        xlabel('X')
        % T0
        subplot(338)
        plot(MSH.XCP,T0,'m--')
        title('T0 (total temperature)')
        xlabel('X')
        % residual
        subplot(339)
        semilogy(RMAX,'k:')
        title(['Rmax, Time=' num2str(SOLV.TIME)])
        xlabel('ISTEP')
        drawnow
    end
    
     % Stop iterating
    if BUTTON.doSTOP == 1
        disp(' ')
        disp('---- STOP ----')
        disp(' ')
        break;
    end
    % Pause execution
    if BUTTON.doPAUSE == 1
        disp(' ')
        disp('---- PAUSE ----')
        disp(' ')
        disp(' hit a key in ''Command Window'' to continue iterating!')
        disp(' ')
        pause
        BUTTON.doPAUSE = 0;
    end
    % Change plotstep
    if BUTTON.changePLOTSTEP ~= 0
        BUTTON.PLOTSTEP = BUTTON.PLOTSTEP + BUTTON.changePLOTSTEP;
        if mod(BUTTON.PLOTSTEP,BUTTON.PLOTSTEPINC) == 1
            BUTTON.PLOTSTEP = BUTTON.PLOTSTEP - 1;
        end
        if BUTTON.PLOTSTEP <= 0
            BUTTON.PLOTSTEP = 1;
        end
        BUTTON.changePLOTSTEP = 0;
        set(BUTTON.htextPLOTSTEP,'String',['Plot every ' num2str(BUTTON.PLOTSTEP) ' timestep'])
    end
end

% Compute auxiliary flow variables
SOLV = auxvr(MSH,GAS,SOLV,1);

% Save solution in "example_3p13.mat"
savefile = 'example_3p13';

eval(['save ' savefile ' SOLV MSH'])
disp(' ')
disp(['Solution saved in ''' savefile '''.'])



