MATLAB lessons N.7 and N.8 - Enrico Bertolazzi

:open:

%
% Compute Ader Flux
%
function [F,S] = ADER2Flux( Q,            ...
                            Xnodes,       ...
                            DT,           ...
                            Flux,         ...
                            SlopeLimiter, ...
                            Leading0,     ...
                            Leading1,     ...
                            Source )

  nr = size(Q,1) ;
  DX = Xnodes(2)-Xnodes(1) ;

  % Reconstruct phase
    % construct limited slopes (given N+2 averages return N slopes)
    DQ = feval( SlopeLimiter, Q ) ;

    % extrapolate from the left
    QL = Q(:,2:end-2) + DQ(:,1:end-1)./2 ;

    % extrapolate from the right
    QR = Q(:,3:end-1) - DQ(:,2:end)./2 ;

    % compute slopes
    DELTA = DQ ./ DX ;

    % left slopes
    DELTAL = DELTA(:,1:end-1) ;

    % right slopes
    DELTAR = DELTA(:,2:end) ;

  % End of Reconstruct phase

  % Compute leading term
    D0 = Leading0( QL, QR ) ;
    D1 = Leading1( DELTAL, DELTAR, D0 ) ;
  % End Compute leading term

  % Evaluate ADER2 state
    D0L = D0 - D1 ./2 ;
    D0R = D0 + D1 ./2 ;
    %
    % In order to eliminate the loop 
    %
    %Qader = D0 ;
    %for k=1:size(Qader,2)
    %  A0         = feval( JFlux, D0(:,k) ) ;
    %  Qader(:,k) = Qader(:,k) + (DT/2) * (feval(Source, D0(:,k)) - A0 * D1(:,k)) ;
    %end
    %
    % The second order approxination A0 * D1 ~ Flux( D0+D1/2 ) - Flux(D0-D1/2)
    %
    Qader = D0 + (DT/2) * ( feval( Source, D0 ) - feval( Flux, D0R ) + feval( Flux, D0L ) ) ;
    F     = feval( Flux, Qader) ;
  % End Evaluate numerical flux

  % Evaluate source
    Qs = Q(:,3:end-2) ; % consider only "true" cells
    DD = DELTA(:,2:end-1) ./2 ; % consider only half of the slope of "true" cells
    %
    % In order to eliminate the loop 
    %
    % for k=1:size(Qs,2)
    %   Qs(:,k) = Qs(:,k) + (DT/2) * ( feval( Source, Qs(:,k) ) - feval( JFlux, Qs(:,k) )*DELTA(:,k+1) ) ;
    % end
    %
    % The second order approxination JFlux(Qs) DELTA = Flux( Qs+DELTA/2 ) - Flux(Qs-DELTA/2)
    %
    DL = Qs - DD ;
    DR = Qs + DD ;  
    
    Qs = Qs + (DT/2) * ( feval( Source, Qs ) - feval( Flux, DR ) +  feval( Flux, DL ) ) ;
    S  = feval( Source, Qs ) ;
  % End Evaluate source

end

:open:

%  General Advancing procedure used to test numerical solver for hyperbolic equations
%
%  Name of program: Advancing
%
%  Developed by: Enrico Bertolazzi
%                Department of Mechanical and Structural Engineering
%                enrico.bertolazzi@ing.unitn.it
%                http://www.ing.unitn.it/~bertolaz/
%
%  Last revision: 2 June 2009
%
%  INPUTS:
%
%     Qinit       = initial cell average values
%     Xnodes      = boundary of the cells
%     InitialTime = starting simulation time
%     FinalTime   = end of simulation time
%     StepTime    = the time difference at which the solution is saved
%     CFLratio    = the Courant–Friedrichs–Lewy ratio at which the solver can advance
%
%     maxSpeed    = function with one vector argument which compute maximum speed signal for each cell.
%                   Usage: feval( maxSpeed, Q ).  
%
%     doStep      = function which perform step of the numerical scheme
%                   Usage: feval( doStep, Time, DT, DX, XC, Q ).
%                         Time = Actual time
%                         DT   = advancing time step
%                         DX   = cells size
%                         XC   = cells center
%                         Q    = cells actual state
%
%  OUTPUTS:
%
%     Tframes     = vector with the time values of saved frame
%     Qframes     = cell array with the solution frames computed
%                   Qframes{1}   = contains the solution at time Tframes(1)
%                   Qframes{2}   = contains the solution at time Tframes(2)
%
%                   ...
%
%                   Qframes{end} = contains the solution at time Tframes(end)
%
function [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                           Xnodes,       ...
                                           InitialTime,  ...
                                           FinalTime,    ...
                                           StepTime,     ...
                                           CFLratio,     ...
                                           maxSpeed,     ...
                                           doStep )

  % Initial conditions are copied in the vector Q
  Q = Qinit ;

  % Compute cell size DX, the size of the cells are assumed all equals.
  DX = Xnodes(2) - Xnodes(1) ;

  % save initial frame
  Qframes{1} = Q ;
  Tframes(1) = InitialTime ;

  kframe   = 1 ; % number of saved framed 
  Time     = InitialTime ; % initial time
  NextTime = InitialTime + StepTime ; % Time to save next frame

  % Time marching procedure, a maximum of 10000 steps are allowed
  for NITER=1:10000

    % Courant-Friedrichs-Lewy (CFL) condition imposed
    DT = CFLratio * DX / max( feval( maxSpeed, Q ) ) ;
        
    % reduce DT if necessary for the last step
    if Time+DT >= FinalTime
      DT = FinalTime - Time ;
      % perform last step
      Q = feval( doStep, Time, DT, Xnodes, Q ) ;
      % save last frame
      kframe          = kframe+1 ;
      Qframes{kframe} = Q ;
      Tframes(kframe) = FinalTime ;
      % exit loop
      break ;
    end

    % adjust time step if necessary to match the time grid
    DT1 = NextTime - Time ;
    while DT >= DT1 % if DT contains a grid step perform smaller step

      % perform step
      Qnew = feval( doStep, Time, DT1, Xnodes, Q ) ;

      % Solution frame is saved
      kframe          = kframe+1 ;
      Qframes{kframe} = Qnew ;
      Tframes(kframe) = NextTime ;

      % check if we have reached the final mesh point
      if FinalTime < NextTime 
        break ;
      end ;

      % update next time frame
      NextTime = NextTime + StepTime ;
      DT1      = NextTime - Time ;
    end

    % perform step
    Q = feval( doStep, Time, DT, Xnodes, Q ) ;

    % update time
    Time = Time + DT ;
  end

end

:open:

%  Name of program: buildXYForPlot
%
%  Developed by: Enrico Bertolazzi
%                Department of Mechanical and Structural Engineering
%                enrico.bertolazzi@ing.unitn.it
%                http://www.ing.unitn.it/~bertolaz/
%
%  Last revision: 2 June 2009
%
%  Given a vector CB of cells boundary edges and a vector Q of cell average values
%  build the vectos x and y which can be used to plot the piecewise constant function
%  which represents the solution
%
%  INPUT
%    Xnodes: cells boundary       (row vector)
%    Q:      cell average value
%
%  OUTPUT
%    X:  x-coordinate 
%    Y:
%
%          +------+------+
%          |      |       |      +------+
%   +------+      |       |      |      |
%   |      |      |       +------+      |
%   | Q(1) | Q(2) | Q(3)  | Q(4) |      |
%   |      |      |       |      |      |
% CB(1)  CB(2)  CB(3)
%
%  Last revision: 2 June 2009
%
function [X,Y] = buildXYForPlot( Xnodes, Q )

  % allocate memory
  N = size(Q,2) ;
  X = zeros( 1,         2*N ) ;
  Y = zeros( size(Q,1), 2*N ) ;

  X(1:2:2*N-1) = Xnodes(1:end-1) ;
  X(2:2:2*N)   = Xnodes(2:end) ;

  Y(:,1:2:2*N-1) = Q ;
  Y(:,2:2:2*N)   = Q ;
                                    
end

:open:

%  Given average states compute WENO limited slopes
%
%  Name of program: WENOSlopeLimiter.m
%
%  Developed by: Enrico Bertolazzi
%                Department of Mechanical and Structural Engineering
%                enrico.bertolazzi@ing.unitn.it
%                http://www.ing.unitn.it/~bertolaz/
%
%  This code is inspired on the library
%    NUMERICA
%    A Library of Source Codes for Teaching,
%    Research and Applications,
%    by E. F. Toro
%    Published by NUMERITEK LTD,
%    Website: www.numeritek.coM
%
%  The algorithm are described in the book:
%
%    E. F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics"
%                 Springer-Verlag, 1997, Second Edition, 1999
%
%  Last revision: 27 May 2009
%
%  The scheme can be visualized in the following way
% 
%   +===+  The cell where the average value Qi is "stored".
%
%   +...+  Ghost cell. These cells are used to extrapolate values for the
%          construction of the bounary conditions. Two ghost cell for boundary
%          are enough for all the considered fluxes.
%
%   +xxx+  Uncomputed values. These cells do not cotains valid values.
%          Are used to avoid shifting of indices or complex indicization.
%
%   Averages states
%
%   +...+...+===+===+===+===+===    ===+===+...+...+
%     Q1  Q2  Q3                                QN+4
%
%   Limited slopes, the slopes S1 and SN+4 are not computed (indeed it is not possible)
%
%   +xxx+===+===+===+===+===+===    ===+===+===+xxx+
%         S1  S2  S3                        SN+1
%
%  INPUTS:
%
%     Q = average values
%
%  OUTPUTS:
%
%     SLOPE = limited slopes
%
%     The slopes are built by finite difference of sucessive averages.
%     In this way slopes associated to interfaces are used to copute cell slopes.
%     The slopes associated to a cell is the minimum slope between the two interface slopes.
%
function DELTA = WENOSlopeLimiter( Q )

  % compute interface slopes, SLOPE(k) is the slopes of interface between cell k and k+1
  SLOPE = Q(:,2:end) - Q(:,1:end-1) ;

  % left slopes
  LSLOPE = SLOPE(:,1:end-1);

  % right slopes
  RSLOPE = SLOPE(:,2:end);

  % find minimum slopes
  DELTA = minmod( LSLOPE, RSLOPE ) ;

  % is equivalent to the loop
  %
  % for k=1:lenght(SLOPE)
  %   if abs(LSLOPE(k)) < abs(RSLOPE(k))
  %     SLOPE(k) = LSLOPE(k);
  %   else
  %     SLOPE(k) = RSLOPE(k);
  %   end
  % end

end

%
% given the vector (matrix) a and b return
% res(k) = sign(a)*min(abs(a(k)),abs(b(k))) if a(k)*b(k)>0
% res(k) = 0 otherwise
%
function res = minmod( a, b )
  res = zeros( size(a) ) ;
  I1 = find( a > 0 & b > 0 ) ;
  I2 = find( a < 0 & b < 0 ) ;
  I3 = find( a.*b <= 0 ) ;
  res(I1) = min(a(I1),b(I1)) ;
  res(I2) = max(a(I2),b(I2)) ;
  res(I3) = 0 ;
end

:open:

%  Test Advancing routine by solving the problem
%
%   Q(x,t) + Q(x,t)  = 0    -1 <= x <= 1
%    t        x
%
%            /
%            |  0 if x < -1/2
%   Q(x,0) = |  1 if x > -1/2 and x < 1/2
%            |  0 if x > 1/2
%             \
%               +---+
%               |   |
%      ---------+   +----------
%
%  Cyclic boundary conditions
%
%  Last revision: 28 May 2009
%
function TestAder2LinearTransport

  N     = 100 ; % number of cells
  N4    = N/4 ;

  % setup intial conditions
  Qinit = [ zeros(1,N4) ones(1,N-2*N4) zeros(1,N4) ] ;
  
  % setup the boundary of the cells
  Xnodes = -1:2/N:1 ;
  
  % initial time and time step
  InitialTime = 0 ;
  FinalTime   = 20 ;
  StepTime    = 0.1 ;
  CFLratio    = 0.95 ;

  [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                    Xnodes,       ...
                                    InitialTime,  ...
                                    FinalTime,    ...
                                    StepTime,     ...
                                    CFLratio,     ...
                                    @maxSpeed,    ...
                                    @doStepADER2 ) ;

  [ Qframes1, Tframes1 ] = Advancing( Qinit,        ...
                                      Xnodes,       ...
                                      InitialTime,  ...
                                      FinalTime,    ...
                                      StepTime,     ...
                                      CFLratio,     ...
                                      @maxSpeed,    ...
                                      @doStepGodunov ) ;

                                % plot the solution at final time step
  [XX,YY] = buildXYForPlot( Xnodes, Qframes{end} ) ;
  [X1,Y1] = buildXYForPlot( Xnodes, Qframes1{end} ) ;
  plot( XX, YY, X1, Y1 ) ;
                             
end
%
% Return a vector with the maximum speed of each cell
%
function SPEED = maxSpeed( Q )
  SPEED = ones(size(Q)) ;
end

%
% Flux for linear transport
%
function F = Flux( Q )
  F = Q ;
end

%
% Jacobians of the flux (vectorial form)
%
function A = JFlux( Q )
  A = ones(size(Q)) ;
end
%
% Solve the Riemann problem for the linear transport equation exactly.
%
function USTAR = Riemann(UL, UR)
  USTAR = UL ;
end

function S = Source( Q )
  S = zeros(size(Q)) ;
end

function D0 = Leading0( QL, QR )
  D0 = Riemann(QL,QR) ;
end

function D1 = Leading1( DL, DR, D0 )
  D1 = Riemann(DL, DR) ;
end
%
% Advancing with ADER2 numerical scheme
%
function Qnew = doStepADER2( Time, DT, Xnodes, Q )

  % extrapolation
  Qe = [ Q(:,end-1) Q(:,end) Q Q(:,1) Q(:,2) ] ;

  [F,S] = ADER2Flux( Qe, Xnodes, DT, @Flux, @WENOSlopeLimiter, @Leading0, @Leading1, @Source ) ;

  DX   = Xnodes(2)-Xnodes(1) ;
  Qnew = Q - (DT/DX) .* ( F(:,2:end) - F(:,1:end-1) ) + DT * S ;

end

%
% Advancing with simple Godunov numerical scheme
%
function Qnew = doStepGodunov( Time, DT, Xnodes, Q )
  % extrapolation (with only one ghost cell by side)
  Qe = [ Q(:,end)  Q Q(:,1) ] ;
  % reconstruction no needed

  % compute numerical flux
  QL    = Qe(:,1:end-1) ; % cell on the left
  QR    = Qe(:,2:end)   ; % cell on the right
  QSTAR = Riemann(QL, QR) ;
  F     = Flux(QSTAR) ;
  S     = Source(Q) ;

  DX   = Xnodes(2) - Xnodes(1) ; % cells are assumed all of equal size
  Qnew = Q - (DT/DX) .* ( F(2:end) - F(1:end-1) ) + DT * S ;

end

:open:

%  Test Advancing routine by solving the problem
%
%   Q (x,t) + Q(x,t) Q (x,t) = 0    -pi <= x <= pi
%    t                x
%
%   Q(x,0) = 1+sin(x)
%
%  Cyclic boundary conditions
%
%  Last revision: 28 May 2009
%
function TestAder2Burger

  N = 100 ; % number of cells

  % setup the boundary of the cells
  DX     = 2*pi/N ;
  Xnodes = [-pi:DX:pi] ;

  % setup intial conditions
  Xc    = (Xnodes(1:end-1)+Xnodes(2:end))/2 ; % center of the cells
  Qinit = 1+sin(Xc) ;

  % initial time and time step
  InitialTime = 0 ;
  FinalTime   = 10 ;
  StepTime    = 0.1 ;
  CFLratio    = 0.9 ;

  [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                    Xnodes,       ...
                                    InitialTime,  ...
                                    FinalTime,    ...
                                    StepTime,     ...
                                    CFLratio,     ...
                                    @maxSpeed,    ...
                                    @doStepADER2 ) ;

  [ Qframes1, Tframes1 ] = Advancing( Qinit,        ...
                                      Xnodes,       ...
                                      InitialTime,  ...
                                      FinalTime,    ...
                                      StepTime,     ...
                                      CFLratio,     ...
                                      @maxSpeed,    ...
                                      @doStepGodunov ) ;

                                % plot the solution at final time step
  [XX,YY] = buildXYForPlot( Xnodes, Qframes{end} ) ;
  [X1,Y1] = buildXYForPlot( Xnodes, Qframes1{end} ) ;
  plot( XX, YY, X1, Y1 ) ;
                             
end
%
% Return a vector with the maximum speed of each cell
%
function SPEED = maxSpeed( Q )
  SPEED = abs(Q) ;
end

function A = JFlux( Q )
  A = Q ;
end

function F = Flux( Q )
  F = (Q .* Q) ./ 2 ;
end
%
% Solve the Riemann problem for the inviscid Burgers equation exactly.
%
function USTAR = Riemann(UL, UR)
  % UL         Left data state
  % UR         Right data state
  % S          Shock speed
  % USTAR      Sampled state
  ISW = find ( UL <  UR ) ;
  IRW = find ( UL >= UR ) ;
  % Solution is a shock wave
    % Compute shock speed S
    S(ISW) = 0.5*(UL(ISW) + UR(ISW)) ;
    % Sample the state along the t-axis
    IGE = find( S(ISW) >= 0 ) ;
    ILT = find( S(ISW) <  0 ) ;
    USTAR(IGE) = UL(IGE) ;
    USTAR(ILT) = UR(ILT) ;
  % Solution is a rarefaction wave.
    % There are 3 cases:
    I1 = find( UL >= 0 ) ;
    I2 = find( UR <= 0 ) ;
    I3 = find( UL < 0 & UR > 0 ) ;
    USTAR(I1) = UL(I1) ; % Right supersonic rarefaction
    USTAR(I2) = UR(I2) ; % Left supersonic rarefaction
    USTAR(I3) = 0 ;      % Transonic rarefaction
end

function S = Source( Q )
  S = zeros(size(Q)) ;
end

function D0 = Leading0( QL, QR )
  D0 = Riemann(QL,QR) ;
end

function D1 = Leading1( QL, QR, D0 )
  D1 = zeros(size(QL)) ;
  I1 = find( D0 >  0 ) ;
  I2 = find( D0 <= 0 ) ;
  D1(I1) = QL(I1) ;
  D1(I2) = QR(I2) ;
  %for k=1:size(QL,2)
  %  D1(k) = LinearRiemannProblem( QL(k), QR(k), A(k) ) ;
  %end
end
%
% Advancing with simple Godunov numerical scheme
%
function Qnew = doStepADER2( Time, DT, Xnodes, Q )

  % extrapolation
  Qe = [ Q(:,end-1) Q(:,end) Q Q(:,1) Q(:,2) ] ;

  [F,S] = ADER2Flux( Qe, Xnodes, DT, @Flux, @WENOSlopeLimiter, @Leading0, @Leading1, @Source ) ;

  DX   = Xnodes(2)-Xnodes(1) ;
  Qnew = Q - (DT/DX) .* ( F(:,2:end) - F(:,1:end-1) ) + DT * S ;

end

%
% Advancing with simple Godunov numerical scheme
%
function Qnew = doStepGodunov( Time, DT, Xnodes, Q )
  % extrapolation (with only one ghost cell by side)
  Qe = [ Q(end)  Q Q(1) ] ;
  % reconstruction no needed

  % compute numerical flux
  QL    = Qe(1:end-1) ; % cell on the left
  QR    = Qe(2:end)   ; % cell on the right
  QSTAR = Riemann(QL, QR) ;
  F     = Flux( QSTAR ) ;

  DX   = Xnodes(2) - Xnodes(1) ; % cells are assumed all of equal size
  Qnew = Q - (DT/DX) .* ( F(2:end) - F(1:end-1) ) ;

end

:open:

%  Test Advancing routine by solving the problem
%
%   Q (x,t) + A Q (x,t)  = 0    -1 <= x <= 1
%    t           x
%
%            /
%            |  [0;0] if x < -1/2
%   Q(x,0) = |  [1;2] if x > -1/2 and x < 1/2
%            |  [0;0] if x > 1/2
%             \
%               +---+
%               |   |
%      ---------+   +----------
%
%  Cyclic boundary conditions
%
%  Last revision: 28 May 2009
%
function TestAder2

  global A_jacobian ;
  
  A_jacobian = [ 0 1 ; 1 0 ] ;

  N     = 100 ; % number of cells
  N4    = N/4 ;

  % setup intial conditions
  Qinit = [ zeros(2,N4) [1;2]*ones(1,N-2*N4) zeros(2,N4) ] ;
  
  % setup the boundary of the cells
  Xnodes = -1:2/N:1 ;
  
  % initial time and time step
  InitialTime = 0 ;
  FinalTime   = 6.3 ;
  StepTime    = 0.1 ;
  CFLratio    = 0.9 ;

  [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                    Xnodes,       ...
                                    InitialTime,  ...
                                    FinalTime,    ...
                                    StepTime,     ...
                                    CFLratio,     ...
                                    @maxSpeed,    ...
                                    @doStepADER2 ) ;

  [ Qframes1, Tframes1 ] = Advancing( Qinit,        ...
                                      Xnodes,       ...
                                      InitialTime,  ...
                                      FinalTime,    ...
                                      StepTime,     ...
                                      CFLratio,     ...
                                      @maxSpeed,    ...
                                      @doStepGodunov ) ;

                                % plot the solution at final time step
  [XX,YY] = buildXYForPlot( Xnodes, Qframes{end} ) ;
  [X1,Y1] = buildXYForPlot( Xnodes, Qframes1{end} ) ;
  plot( XX, YY, X1, Y1 ) ;
                                    
end

%
% Return a vector with the maximum speed of each cell
%
function SPEED = maxSpeed( Q )
  SPEED = ones(1,size(Q,2)) ;
end

function F = Flux( Q )
  F = JFlux( Q ) * Q ;
end

function A = JFlux( Q )
  global A_jacobian ;
  A = A_jacobian ;
end

function S = Source( Q )
  S = zeros(size(Q)) ;
end

function D0 = Leading0( QL, QR )
  global A_jacobian ;
  D0 = LinearRiemannProblem( QL, QR, A_jacobian ) ;
end

function D1 = Leading1( QL, QR, D0 )
  global A_jacobian ;
  D1 = LinearRiemannProblem( QL, QR, A_jacobian ) ;
end
%
% Advancing with simple Godunov numerical scheme
%
function Qnew = doStepADER2( Time, DT, Xnodes, Q )

  % extrapolation
  Qe = [ Q(:,end-1) Q(:,end) Q Q(:,1) Q(:,2) ] ;

  [F,S] = ADER2Flux( Qe, Xnodes, DT, @Flux, @WENOSlopeLimiter, @Leading0, @Leading1, @Source ) ;

  DX   = Xnodes(2)-Xnodes(1) ;
  Qnew = Q - (DT/DX) .* ( F(:,2:end) - F(:,1:end-1) ) + DT * S ;

end

%
% Advancing with simple Godunov numerical scheme
%
function Qnew = doStepGodunov( Time, DT, Xnodes, Q )
  % extrapolation (with only one ghost cell by side)
  Qe = [ Q(:,end)  Q Q(:,1) ] ;
  % reconstruction no needed

  % compute numerical flux
  QL    = Qe(:,1:end-1) ; % cell on the left
  QR    = Qe(:,2:end)   ; % cell on the right
  QSTAR = Leading0(QL, QR) ;
  F     = Flux( QSTAR ) ;

  DX   = Xnodes(2) - Xnodes(1) ; % cells are assumed all of equal size
  Qnew = Q - (DT/DX) .* ( F(:,2:end) - F(:,1:end-1) ) ;

end

:open:

%  Test ADER2 by solving the problem
%
%   Q(x,t) + Q(x,t)  = 0    -1 <= x <= 1
%    t        x
%
%   Q(x,0) = 1+sin(x)
%
% With Cyclic boundary conditions
%
%  Last revision: 28 May 2009
%
function TestAder2Order

  NN = [ 50 100 200 400 800 ] ;

  for k=1:4
    
    N = NN(k) ;% number of cells

    % setup the boundary of the cells
    DX     = 2*pi/N ;
    Xnodes = [-pi:DX:pi] ;

    % setup intial conditions
    Xc    = (Xnodes(1:end-1)+Xnodes(2:end))/2 ; % center of the cells
    Qinit = Q0(Xc) ;
  
    % initial time and time step
    InitialTime = 0 ;
    FinalTime   = 20 ;
    StepTime    = 1 ;
    CFLratio    = 0.95 ;
  
    Qesatta = Q0(Xc - 20) ;

    [ Qframes, Tframes ] = Advancing( Qinit,        ...
                                      Xnodes,       ...
                                      InitialTime,  ...
                                      FinalTime,    ...
                                      StepTime,     ...
                                      CFLratio,     ...
                                      @maxSpeed,    ...
                                      @doStepADER2 ) ;
                                  
    % plot the solution at final time step
    [XX,YY] = buildXYForPlot( Xnodes, Qframes{end} ) ;
    [X1,Y1] = buildXYForPlot( Xnodes, Qesatta ) ;
    plot( XX, YY, X1, Y1) ;
  
    err(k) =  max(abs(YY-Y1)) ;
    disp( sprintf( 'Errore %g', err(k) ) ) ;
  end
  
  disp( sprintf( 'Order %g', log(err(1)/err(2))/log(2) ) ) ;  
  disp( sprintf( 'Order %g', log(err(2)/err(3))/log(2) ) ) ;
  disp( sprintf( 'Order %g', log(err(3)/err(4))/log(2) ) ) ;

end

%
%
%
function Q = Q0( Xc ) 
  Q = 1+sin(Xc) ;
end


%
% Return a vector with the maximum speed of each cell
%
function SPEED = maxSpeed( Q )
  SPEED = ones(size(Q)) ;
end

%
% Flux for linear transport
%
function F = Flux( Q )
  F = Q ;
end

%
% Jacobians of the flux (vectorial form)
%
function A = JFlux( Q )
  A = ones(size(Q)) ;
end
%
% Solve the Riemann problem for the linear transport equation exactly.
%
function USTAR = Riemann(UL, UR)
  USTAR = UL ;
end

function S = Source( Q )
  S = zeros(size(Q)) ;
end

function D0 = Leading0( QL, QR )
  D0 = Riemann(QL,QR) ;
end

function D1 = Leading1( DL, DR, D0 )
  D1 = Riemann(DL, DR) ;
end
%
%  Given average states compute unlimited slopes
%
%  INPUTS:
%
%     Q = average values
%
%  OUTPUTS:
%
%     DELTA = slopes
%
%     The slopes are built by finite difference of sucessive averages.
%     In this way slopes associated to interfaces are used to copute cell slopes.
%     The slopes associated to a cell is the minimum slope between the two interface slopes.
%
function DELTA = SlopeNoLimiter( Q )
  % compute interface slopes, SLOPE(k) is the slopes of interface between cell k and k+1
  DELTA = (Q(:,3:end) - Q(:,1:end-2))./2 ;
end

%
% given the vector (matrix) a and b return
% res(k) = sign(a)*min(abs(a(k)),abs(b(k))) if a(k)*b(k)>0
% res(k) = 0 otherwise
%
function res = minmod( a, b )
  res = zeros( size(a) ) ;
  I1 = find( a > 0 & b > 0 ) ;
  I2 = find( a < 0 & b < 0 ) ;
  I3 = find( a.*b <= 0 ) ;
  res(I1) = min(a(I1),b(I1)) ;
  res(I2) = max(a(I2),b(I2)) ;
  res(I3) = 0 ;
end

%
% Advancing with ADER2 numerical scheme
%
function Qnew = doStepADER2( Time, DT, Xnodes, Q )

  % extrapolation
  Qe = [ Q(:,end-1) Q(:,end) Q Q(:,1) Q(:,2) ] ;

  [F,S] = ADER2Flux( Qe, Xnodes, DT, @Flux, @SlopeNoLimiter, @Leading0, @Leading1, @Source ) ;

  DX   = Xnodes(2)-Xnodes(1) ;
  Qnew = Q - (DT/DX) .* ( F(:,2:end) - F(:,1:end-1) ) + DT * S ;

end

MATLAB lessons N.7 and N.8¶