%{
  FINITE-DIFFERENCE SOLUTION OF THE QUASI-1D EULER EQUATIONS
  FOR AN IDEAL GAS. MAC-CORMACK TIME-STEPPING SCHEME.

  M. Piller.
  Started 23/05/2023.

  How to call the program:

  1. Set parameters for the considered flow scenario either in
     "subsonic_subsonic_isoentropic_flow_analytical_solution"
     or in
     "subsonic_supersonic_isoentropic_flow_analytical_solution"
  2. Set additional parameters in "startup"
  3. Run the program:
  % whichCase = "subSonic-superSonic";
  whichCase = "subSonic-subSonic";
  restart   = 0;
  [t,x,P,T,rho,u,c,Ma]=nozzle(whichCase,restart);

  %}
  function [t,x,P,T,rho,u,c,Ma]=nozzle(whichCase,restart)

    close all; clc;

    %% *********************************************************
    %% SET-UP THE CASE
    %% *********************************************************
    switch whichCase
      case{'subSonic-subSonic'}
      [subsubCase,ierr]=subsonic_subsonic_isoentropic_flow_analytical_solution;
      case{'subSonic-superSonic'}
      [subsubCase,ierr]=subsonic_supersonic_isoentropic_flow_analytical_solution;
  endswitch

  % INITIALIZATION ROUTINE
  [L,x,At,A,N,h,crat,GasR,T0,P0,rho0,a0,...
  A_star,Pstar,Pe,Me,rho,u,P,T,t,tf,analyticalSolution,...
  Cx,CFL_lim,CFL_local,press,dens,soundSpeed,bconds,INT,Istore,nitMax] = startup(subsubCase,restart,ierr);

  hf = figure('color',[0.5 0.5 0.5],'nextplot','replace');

  % START TIME-STEPPING
  STORED = [];
  nit = 0;
  while (min(t)<tf)&&(nit<nitMax)

    nit = nit+1;

    % SET TIME-STEP TO RESPECT CFL CONDITION
    dt = set_time_step(u,T,h,crat,GasR,soundSpeed,INT,CFL_lim,CFL_local);

    t = t + dt;

    % PREDICTOR-CORRECTOR SCHEME
    [P,T,rho,u,c,Ma] = time_stepping(x,u,rho,P,T,A,P0,T0,rho0,Pe,GasR,crat,Cx,dens,press,soundSpeed,bconds,INT,t,dt,nit);

    % STORE VARIABLES AT SELECTED POSITIONS FOR TIME-TRACING
    STORED = [STORED;[t(Istore),P(Istore),T(Istore),rho(Istore),u(Istore),c(Istore),Ma(Istore),rho(Istore).*u(Istore).*A(Istore)]];

    if (mod(nit,1000)==0)||(min(t)>=tf)
      figure(hf);
      subplot(2,3,1);hp=plot(x,u./a0,'bo-');grid on;box on;xl=xlabel('x');yl=ylabel('u'); ht=text(L/2,mean(u),['(Mean) Time: ',num2str(mean(t)),' s']);set(gca,'color',[0.6 0.6 0.6]);
      subplot(2,3,2);hp=plot(x,rho./rho0,'bo-');grid on;box on;xl=xlabel('x');yl=ylabel('rho'); ht=text(L/2,mean(rho),['(Mean) Time: ',num2str(mean(t)),' s']);set(gca,'color',[0.6 0.6 0.6]);
      subplot(2,3,3);hp=plot(x,T./T0,'bo-');grid on;box on;xl=xlabel('x');yl=ylabel('T'); ht=text(L/2,mean(T),['(Mean) Time: ',num2str(mean(t)),' s']);set(gca,'color',[0.6 0.6 0.6]);
      subplot(2,3,4);hp=plot(x,P./P0,'bo-');grid on;box on;xl=xlabel('x');yl=ylabel('P'); ht=text(L/2,mean(P),['(Mean) Time: ',num2str(mean(t)),' s']);set(gca,'color',[0.6 0.6 0.6]);
      subplot(2,3,5);hp=plot(x,Ma,'bo-');grid on;box on;xl=xlabel('x');yl=ylabel('Mach'); ht=text(L/2,mean(Ma),['(Mean) Time: ',num2str(mean(t)),' s']);set(gca,'color',[0.6 0.6 0.6]);
      subplot(2,3,6);hp=plot(x,c./a0,'bo-');grid on;box on;xl=xlabel('x');yl=ylabel('c'); ht=text(L/2,mean(c),['(Mean) Time: ',num2str(mean(t)),' s']);set(gca,'color',[0.6 0.6 0.6]);
      nit
    endif

  endwhile

  % SAVE WHAT YOU NEED FOR RESTARTING THE SOLUTION FROM THE PRESENT STATE
  save "restart.mat" u rho P T t

  figure('color',[0.5 0.5 0.5])
  hold on
  n = 1;
  for j=1:length(Istore)
    hp=plot(STORED(:,j).*((x(end)-x(1))/a0),STORED(:,j+n*length(Istore))./P0,'b-');
    set(hp,'linewidth',2.0);
  endfor
  set(gca,'color',[0.6 0.6 0.6]);
  grid on;
  box on;
  xl = xlabel('tL/a_0');
  yl = ylabel('P/P_0');

  figure('color',[0.5 0.5 0.5])
  hold on
  n = 2;
  for j=1:length(Istore)
    hp=plot(STORED(:,j).*((x(end)-x(1))/a0),STORED(:,j+n*length(Istore))./T0,'b-');
    set(hp,'linewidth',2.0);
  endfor
  set(gca,'color',[0.6 0.6 0.6]);
  grid on;
  box on;
  xl = xlabel('tL/a_0');
  yl = ylabel('T/T_0');

  figure('color',[0.5 0.5 0.5])
  hold on
  n = 3;
  for j=1:length(Istore)
    hp=plot(STORED(:,j).*((x(end)-x(1))/a0),STORED(:,j+n*length(Istore))./rho0,'b-');
    set(hp,'linewidth',2.0);
  endfor
  set(gca,'color',[0.6 0.6 0.6]);
  grid on;
  box on;
  xl = xlabel('tL/a_0');
  yl = ylabel('\rho/\rho_0');

  figure('color',[0.5 0.5 0.5])
  hold on
  n = 4;
  for j=1:length(Istore)
    hp=plot(STORED(:,j).*((x(end)-x(1))/a0),STORED(:,j+n*length(Istore))./a0,'b-');
    set(hp,'linewidth',2.0);
  endfor
  set(gca,'color',[0.6 0.6 0.6]);
  grid on;
  box on;
  xl = xlabel('tL/a_0');
  yl = ylabel('u/u_0');

  figure('color',[0.5 0.5 0.5])
  hold on
  n = 6;
  for j=1:length(Istore)
    hp=plot(STORED(:,j).*((x(end)-x(1))/a0),STORED(:,j+n*length(Istore)),'b-');
    set(hp,'linewidth',2.0);
  endfor
  set(gca,'color',[0.6 0.6 0.6]);
  grid on;
  box on;
  xl = xlabel('tL/a_0');
  yl = ylabel('Ma');

  figure('color',[0.5 0.5 0.5])
  hold on
  n = 7;
  for j=1:length(Istore)
    hp=plot(STORED(:,j).*((x(end)-x(1))/a0),STORED(:,j+n*length(Istore))./(rho0*a0*min(A)),'b-');
    set(hp,'linewidth',2.0);
  endfor
  set(gca,'color',[0.6 0.6 0.6]);
  grid on;
  box on;
  xl = xlabel('tL/a_0');
  yl = ylabel('(\rho u A/(\rho_0 a_0 A_t)');

endfunction

%% *********************************************************************************************************** %%
%% *********************************************************************************************************** %%
%% *********************************************************************************************************** %%
%%                                  SUBROUTINES                                                                %%
%% *********************************************************************************************************** %%
%% *********************************************************************************************************** %%
%% *********************************************************************************************************** %%

% SET TIME-STEP
function dt = set_time_step(u,T,h,crat,GasR,soundSpeed,INT,CFL_lim,CFL_local)

  c  = soundSpeed(T,crat,GasR);

  % local time-step
  dt = CFL_lim*0.5*(h(1:end-1)+h(2:end))./(abs(u(INT))+c(INT));
  if ~CFL_local
    dt = min(dt)*ones(size(dt));
  endif

endfunction

% STARTUP ROUTINE
function [L,x,At,A,N,h,crat,GasR,T0,P0,rho0,a0,A_star,Pstar,Pe,Me,rho,u,P,T,t,tf,analyticalSolution,...
  Cx,CFL_lim,CFL_local,press,dens,soundSpeed,bconds,INT,Istore,nitMax] = startup(subsubCase,restart,ierr)

  if ierr
    disp(['ERROR IN SETTING UP THE CASE']);
    return;
  endif

  % Grid points where the solution has to be stored
  Istore = [15];

  % Geometry
  L      = subsubCase.L;
  x      = subsubCase.x(:);
  At     = subsubCase.At;
  A      = subsubCase.A(:);
  N      = length(x);
  h      = x(2:end)-x(1:end-1);

  % Type of gas
  crat   = subsubCase.g;
  GasR   = subsubCase.R;

  % Reservoir conditions
  T0   = subsubCase.T0;
  P0   = subsubCase.P0;
  rho0 = subsubCase.rho0;
  a0   = subsubCase.a0;

  % Sonic conditions
  A_star = subsubCase.A_star;
  Pstar  = subsubCase.Pstar;

  % Outflow conditions
  Pe = subsubCase.Pe;
  Me = subsubCase.Me;

  % Final time
  tf = subsubCase.tf;

  % Initial conditions
  if restart
    ic = dir('restart.mat');
    if isempty(ic)
      disp(['NO RESTART FILE']);
      return;
    endif
    load("restart.mat");
  else

    rho     = subsubCase.rho_i(x,L,rho0);
    u       = subsubCase.V_i(x,L,T0,a0);
    T       = subsubCase.T_i(x,L,T0);
    if isfield(subsubCase,"P_i")
      P = subsubCase.P_i(x,L,P0,T0,GasR);
    else
      P = rho.*GasR.*T;
    endif

    t = zeros(N-2,1);
  endif

  % Analytical solution
  analyticalSolution.Mach = subsubCase.Mach;
  analyticalSolution.P      = subsubCase.P;
  analyticalSolution.T      = subsubCase.T;
  analyticalSolution.rho    = subsubCase.rho;

  % Numerics
  Cx        = 0;
  CFL_lim   = 0.5;
  CFL_local = 0;
  nitMax    = 20000; % Max Nr of time-steps

  % Boundary conditions
  bconds = subsubCase.bconds;

  % Utilities
  press      = @(rho,T,R) rho.*T.*R;
  dens       = @(P,T,R) P./(R*T);
  soundSpeed = @(T,g,R) sqrt(g*R*T);

  INT = 2:N-1;

endfunction

function [P,T,rho,u,c,Ma] = time_stepping(x,u,rho,P,T,A,P0,T0,rho0,Pe,GasR,crat,Cx,dens,press,soundSpeed,bconds,INT,t,dt,nit)

  N = length(x);

  rho_bar = zeros(N,1);
  u_bar   = zeros(N,1);
  T_bar   = zeros(N,1);
  P_bar   = zeros(N,1);

  %%% **************************************************************** %%%
  %%%                       PREDICTOR STEP                             %%%
  %%% **************************************************************** %%%

  % Scheme for spatial derivatives
  if mod(nit,2)==0
    scheme  = 'b';
  else
    scheme = 'f';
  endif

  % Discretization of the Euler equations to get the (approximated) time rate-of-change of rho, u, T
  [dudt,drhodt,dTdt,S,drho2,du2,dT2] = calc_time_derivatives(x,u,rho,P,T,A,GasR,crat,Cx,INT,scheme);

  % Predicted values (accouting also for an artificial dissipation term)
  rho_bar(INT) = rho(INT) + dt.*drhodt + S.*drho2;
  u_bar(INT)   = u(INT)   + dt.*dudt   + S.*du2;
  T_bar(INT)   = T(INT)   + dt.*dTdt   + S.*dT2;
  P_bar(INT)   = press(rho_bar(INT),T_bar(INT),GasR);

  %% Boundary conditions
  [P_bar,T_bar,rho_bar,u_bar] = enforce_bconds(bconds,P_bar,T_bar,rho_bar,u_bar,P0,T0,rho0,GasR,crat,dens,press,Pe);

  %%% **************************************************************** %%%
  %%%                       CORRECTOR STEP                             %%%
  %%% **************************************************************** %%%

  % Scheme for spatial derivatives (base Mac-Cormack scheme alternates forward- and
  % backward finite-differences)
  if mod(nit,2)~=0
    scheme  = 'b';
  else
    scheme = 'f';
  endif

  % Discretization of the Euler equations to get the (approximated) time rate-of-change of rho, u, T
  [dudt_bar,drhodt_bar,dTdt_bar,S,drho2,du2,dT2] = calc_time_derivatives(x,u_bar,rho_bar,P_bar,T_bar,A,GasR,crat,Cx,INT,scheme);

  % Time-derivatives of dependent variables: Cranck-Nicolson scheme
  drhodt = 0.5*(drhodt + drhodt_bar);
  dudt   = 0.5*(dudt   + dudt_bar);
  dTdt   = 0.5*(dTdt   + dTdt_bar);

  % Corrected values (accouting also for an artificial dissipation term)
  rho(INT) = rho(INT) + dt.*drhodt + S.*drho2;
  u(INT)   = u(INT)   + dt.*dudt   + S.*du2;
  T(INT)   = T(INT)   + dt.*dTdt   + S.*dT2;
  P(INT)   = press(rho(INT),T(INT),GasR);

  %% Boundary conditions
  [P,T,rho,u] = enforce_bconds(bconds,P,T,rho,u,P0,T0,rho0,GasR,crat,dens,press,Pe);

  % Mach number and speed of sound
  c  = soundSpeed(T,crat,GasR);
  Ma = u./c;

endfunction





% Boundary conditions. Different flow scenarios require different boundary conditions.
% The flow scenario is selecte by the input variable "bconds"
function [P,T,rho,u] = enforce_bconds(bconds,P,T,rho,u,P0,T0,rho0,GasR,crat,dens,press,Pe)

  N = length(T);

  switch bconds
    % Subsonic-subsonic nozzle
    case{'sub-sub'}
    % Subsonic inlet
    P(1)   = P0;
    T(1)   = T0;
    rho(1) = dens(P(1),T(1),GasR);
    u(1)   = 2.0*u(2)-u(3);

    % Subsonic outlet
    P(N)   = Pe;
    T(N)   = 2.0*T(N-1)-T(N-2);
    rho(N) = dens(P(N),T(N),GasR);
    u(N)   = 2.0*u(N-1)-u(N-2);
    % Subsonic-supersonic nozzle
    case{'sub-super'}
    % Subsonic inlet
    P(1)   = P0;
    T(1)   = T0;
    rho(1) = dens(P(1),T(1),GasR);
    u(1)   = 2.0*u(2)-u(3);

    % Supersonic outlet
    T(N)   = 2.0*T(N-1)-T(N-2);
    rho(N) = 2.0*rho(N-1)-rho(N-2);;
    u(N)   = 2.0*u(N-1)-u(N-2);
    P(N)   = press(rho(N),T(N),GasR);

endswitch

endfunction

% Discretize Euler equations (explicit approach)
function [dudt,drhodt,dTdt,S,drho2,du2,dT2] = calc_time_derivatives(x,u,rho,P,T,A,GasR,crat,Cx,INT,scheme)

dudx      = fd(x,u,scheme);
drhodx    = fd(x,rho,scheme);
drhouAdx  = fd(x,rho.*u.*A,scheme);
dTdx      = fd(x,T,scheme);
dlnAdx    = fd(x,log(A),scheme);

% Time-derivatives of dependent variables
drhodt = -rho(INT).*dudx(INT)-rho(INT).*u(INT).*dlnAdx(INT)-u(INT).*drhodx(INT);
dudt   = -u(INT).*dudx(INT)-GasR*(T(INT).*drhodx(INT)./rho(INT)+dTdx(INT));
dTdt   = -u(INT).*dTdx(INT)-(crat-1)*T(INT).*(dudx(INT)+u(INT).*dlnAdx(INT));

% Artificial dissipation
du2   = (u(3:end)-2*u(2:end-1)+u(1:end-2));
drho2 = (rho(3:end)-2*rho(2:end-1)+rho(1:end-2));
dT2   = (T(3:end)-2*T(2:end-1)+T(1:end-2));

S = Cx*abs(P(3:end)-2*P(2:end-1)+P(1:end-2))./(P(3:end)+2*P(2:end-1)+P(1:end-2));

endfunction

function A = nozzleShape(x,A_star,L)

flag = ones(size(x));
flag(x./L>0.5) = 0;
A = (1+2.2*(3.0*x./L-1.5).^2).*flag + (1+0.2223*(3.0*x./L-1.5).^2).*(1-flag);
A = A.*A_star;

endfunction

function A = nozzleShape2(x,A_star,L)

A = 1+2.2*(3.0*x./L-1.5).^2;
A = A.*A_star;

endfunction

function [subsubCase,ierr]=subsonic_subsonic_isoentropic_flow_analytical_solution

ierr = 0;

% Finel time
tf = 100;

g     = 1.4;
R     = 287;
P0    = 500000;
T0    = 300;

Pstar = P0*(1+0.5*(g-1))^(-g/(g-1));
Pe    = 0.93*P0;

rho0 = P0/(R*T0);
a0   = sqrt(g*R*T0);


if (Pe<Pstar)
  disp(['OUTLET PRESSURE IS TOO LARGE FOR A SUBSONIC-SUBSONIC ISOENTROPIC FLOW!'])
  ierr = -1;
  return;
endif

Me = sqrt(((Pe/P0)^((1-g)/g) - 1)*(2/(g-1)));

L  = 0.25;
At = 0.0001;
N  = 31;
x  = 0:(L/(N-1)):L;

A      = nozzleShape(x,At,L);
A_star = A(end)*Me^2/(((2/(g+1))*(1+0.5*(g-1)*Me^2))^((g+1)/(g-1)));

% Initial conditions
rho_i = @(x,L,rho0) (1.0-0.023*(3.0*x./L)).*rho0;
T_i   = @(x,L,T0)   (1.0-0.009333*(3.0*x./L)).*T0;
V_i   = @(x,L,T0,a0)   (0.05+0.11*(3.0*x./L)).*a0;
P_i   = @(x,L,P0,T0,R)   rho_i(x,L,P0./(R*T0)).*R.*T_i(x,L,T0);

% Mach number distribution
f = @(M,A,A_star,g) (A/A_star)^2-(1/M^2)*((2/(g+1))*(1+0.5*(g-1)*M^2))^((g+1)/(g-1));
Mach = nan(size(x));
options = optimset('TolX',1e-9);
for j=1:length(x)
  A_j = A(j);
  Mach(j) = fzero(@(M)f(M,A_j,A_star,g),Me,options);
endfor


% Pressure, density, temperature distributions
P   = @(x,P0,g,M) P0*(1+0.5*(g-1)*M.^2).^(-g/(g-1));
rho = @(x,rho0,g,M) rho0*(1+0.5*(g-1)*M.^2).^(-1/(g-1));
T   = @(x,T0,g,M) T0*(1+0.5*(g-1)*M.^2).^(-1);

% What boundary conditions?
bconds = 'sub-sub';

subsubCase.tf     = tf;
subsubCase.x      = x;
subsubCase.L      = L;
subsubCase.At     = At;
subsubCase.A      = A;
subsubCase.A_star = A_star;
subsubCase.g      = g;
subsubCase.R      = R;
subsubCase.T0     = T0;
subsubCase.P0     = P0;
subsubCase.rho0   = rho0;
subsubCase.a0     = a0;
subsubCase.Pstar  = Pstar;
subsubCase.Pe     = Pe;
subsubCase.Me     = Me;
subsubCase.rho_i  = rho_i;
subsubCase.V_i    = V_i;
subsubCase.T_i    = T_i;
subsubCase.P_i    = P_i;
subsubCase.Mach   = Mach;
subsubCase.P      = P;
subsubCase.T      = T;
subsubCase.rho    = rho;
subsubCase.bconds = bconds;

endfunction

function [subsubCase,ierr]=subsonic_supersonic_isoentropic_flow_analytical_solution

ierr = 0;

tf = 30;

g     = 1.4;
R     = 287;
P0    = 500000;
T0    = 300;

Pstar = P0*(1+0.5*(g-1))^(-g/(g-1));
Pe    = 0.9*P0;

rho0 = P0/(R*T0);
a0   = sqrt(g*R*T0);


if (Pe<Pstar)
  disp(['OUTLET PRESSURE IS TOO LARGE FOR A SUBSONIC-SUBSONIC ISOENTROPIC FLOW!'])
  ierr = -1;
  return;
endif

Me = sqrt(((Pe/P0)^((1-g)/g) - 1)*(2/(g-1)));

L  = 0.25;
At = 0.0001;
N  = 31;
x  = 0:(L/(N-1)):L;

% For subsonic/supersonic nozzle, A_star == At
A      = nozzleShape2(x,At,L);
A_star = At;

% Initial conditions
rho_i = @(x,L,rho0) (1.0-0.3146*(3.0*x./L)).*rho0;
T_i   = @(x,L,T0)   (1.0-0.2314*(3.0*x./L)).*T0;
V_i   = @(x,L,T0,a0)   (0.1+1.09*(3.0*x./L)).*sqrt(T_i(x,L,T0)/T0).*a0;
P_i   = @(x,L,P0,T0,R)   rho_i(x,L,P0./(R*T0)).*R.*T_i(x,L,T0);

% Mach number distribution
f = @(M,A,A_star,g) (A/A_star)^2-(1/M^2)*((2/(g+1))*(1+0.5*(g-1)*M^2))^((g+1)/(g-1));
Mach = nan(size(x));
options = optimset('TolX',1e-9);
for j=1:length(x)
  A_j = A(j);
  Mach(j) = fzero(@(M)f(M,A_j,A_star,g),1,options);
endfor


% Pressure, density, temperature distributions at steady state
P   = @(x,P0,g,M) P0*(1+0.5*(g-1)*M.^2).^(-g/(g-1));
rho = @(x,rho0,g,M) rho0*(1+0.5*(g-1)*M.^2).^(-1/(g-1));
T   = @(x,T0,g,M) T0*(1+0.5*(g-1)*M.^2).^(-1);

% What boundary conditions?
bconds = 'sub-super';

subsubCase.tf     = tf;
subsubCase.x      = x;
subsubCase.L      = L;
subsubCase.At     = At;
subsubCase.A      = A;
subsubCase.A_star = A_star;
subsubCase.g      = g;
subsubCase.R      = R;
subsubCase.T0     = T0;
subsubCase.P0     = P0;
subsubCase.rho0   = rho0;
subsubCase.a0     = a0;
subsubCase.Pstar  = Pstar;
subsubCase.Pe     = Pe;
subsubCase.Me     = Me;
subsubCase.rho_i  = rho_i;
subsubCase.V_i    = V_i;
subsubCase.T_i    = T_i;
subsubCase.P_i    = P_i;
subsubCase.Mach   = Mach;
subsubCase.P      = P;
subsubCase.T      = T;
subsubCase.rho    = rho;
subsubCase.bconds = bconds;

endfunction

  function df = fd(x,u,dir)

    df = nan(size(u));

    switch dir
      case{'f'}
      df(2:end-1) = (u(3:end)-u(2:end-1))./(x(3:end)-x(2:end-1));
      case{'b'}
      df(2:end-1) = (u(2:end-1)-u(1:end-2))./(x(2:end-1)-x(1:end-2));
  endswitch

  endfunction