%
%     Numerical scheme for the Conservation Law
%
%     Name of program: Solver.m
%
%     Purpose: ........
%
%     Programer: Enrico Bertolazzi 
%               (based on the Fortran NUMERICA lib of prof. E. F. Toro )
%
%     Last revision: 11 April 2008
%
%     Theory is found in Chaps. 5, 13 and 14 of Reference 1
%     and in original references therein
%
%     1. E. F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics"
%        Springer-Verlag, 1997, Second Edition, 1999
%
%     This program is inspired from the NUMERICA library:
%
%     NUMERICA
%     A Library of Source Codes for Teaching,
%     Research and Applications,
%     by E. F. Toro
%     Published by NUMERITEK LTD,
%     Website: www.numeritek.coM
%
%  INPUT:
%     PRB A Structure describing the problem with fields
%       Uinit     = Initial solution
%       Xnodes    = Boundary cell nodes
%       StepTime  = solution is saved every StepTime
%       FinalTime = solution is computed up to FinalTime
%       CFL       = advancing time step CFL
%
%     MODEL A function with variable arguments.
%       First argument simulate method name in Object Oriented Paradigm
%       MODEL must support the following call:
%       Uextr = MODEL( 'extrapolate', U, left_bc, right_bc )
%               Extrapolate 4 cell of states(2 on left 2 on right) for boudary condition
%       S     = MODEL( 'max_speed', U )
%               Evaluate maximum speed signal of the flux for CFL computation
%       F     = MODEL( 'flux', U )
%               Ealuate model flux
%  OUTPUT:
%     UFRAMES   = matrix with the frames of the computed solution.
%                 UFRAMES(:,:,k) is the solution at time level k-th
%     TFRAMES   = vector with the time frames. TFRAMES(i) the time of the computed frame UFRAMES(i,:).
%
function [ UFRAMES, TFRAMES ] = Solver(PRB,MODEL,NFLUX)

  % Initial conditions are set up
  U  = PRB . Uinit ;

  % Compute cell size
  DX = PRB . Xnodes(2) - PRB . Xnodes(1) ; % equal spaced nodes!!

  % save initial frame
  UFRAMES(:,:,1) = U ;
  TFRAMES(1)     = 0 ;

  kframe   = 1            ; % number of saved framed 
  Time     = 0            ; % initial time
  NextTime = PRB.StepTime ; % Time to save next frame

  % Time marching procedure
  for NITER=1:10000
  
    % Extrapolate for Boundary conditions
    Uextra = feval( MODEL, 'extrapolate', U, PRB.left_bc, PRB.right_bc ) ;

    % Courant-Friedrichs-Lewy (CFL) condition imposed
    DT   = DX * PRB.CFL / feval( MODEL, 'max_speed', U ) ;
    
    % update time
    Time = Time + DT ;

    % adjust time step if necessary to match the time grid
    while Time >= NextTime
      % reduce DT in order to match NextTime
      DT1   = DT - ( Time - NextTime ) ;

      % Compute intercell flux FLUX
      FLUX = feval( NFLUX, MODEL, DT1, DX, Uextra ) ;

      kframe = kframe+1 ;
      % Solution is updated according to conservative formula and saved in the frame
      UFRAMES(:,:,kframe) = U + (DT/DX) * (FLUX(:,1:end-1) - FLUX(:,2:end))  ;
      TFRAMES(kframe)     = NextTime ;

      % update next time frame
      NextTime = NextTime + PRB.StepTime ;
    end

    % Compute intercell flux FLUX
    FLUX = feval( NFLUX, MODEL, DT, DX, Uextra ) ;
    % Solution is updated according to conservative formula
    U = U + (DT/DX) * ( FLUX(:,1:end-1) - FLUX(:,2:end) ) ;

    if Time >= PRB.FinalTime ; break ; end
  end
  
  % save last frame
  kframe              = kframe+1 ;
  UFRAMES(:,:,kframe) = U ;
  TFRAMES(kframe)     = Time ;
end
