%
%     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
%     Ncells    = number of cell discretizing [-1,1] the domain of the problems
%     CFL       = time step advancing CFL ratio
%     StepTime  = solution is saved every StepTime
%     FinalTime = solution is computed up to FinalTime
%
%  OUTPUT:
%     X         = the cell boundary of the subdivided interval [-1,1]; X(1) = -1,..., X(N) = 1
%     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)

  % Initial conditions are set up
  U  = feval( MODEL, 'init', PRB.Ncells ) ;
  DX = feval( MODEL, 'get_dx' ) ;

  kframe   = 1 ;
  Time     = 0 ;
  NextTime = PRB.StepTime ;

  % save initial frame
  UFRAMES(:,:,kframe) = U ;
  TFRAMES(kframe)     = Time ;

  % Time marching procedure
  for NITER=1:10000
  
    % Extrapolate for Boundary conditions
    [ ULL, UL, UR, URR ] = feval( MODEL, 'extrapolate_for_bc', U ) ;
    
    Uold = [ ULL, UL, U, UR, URR ] ;

    % Courant-Friedrichs-Lewy (CFL) condition imposed
    DT   = feval( MODEL, 'dt_from_cfl', U, PRB.CFL ) ;
    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 ) ;
      Time1 = NextTime ;

      U1 = Update(MODEL, DT1, DX, Uold ) ;

      kframe              = kframe+1 ;
      UFRAMES(:,:,kframe) = U1 ;
      TFRAMES(kframe)     = Time1 ;

      NextTime = NextTime + PRB.StepTime ;

    end

    % Solution is updated according to conservative formula
    U = Update(MODEL, DT, DX, Uold ) ;

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

function Unew = Update( MODEL, DT, DX, U )
  % Compute intercell flux FLUX
  FLUX = feval( MODEL, 'nflux', DT, U ) ;
  % Solution is updated according to conservative formula
  Unew = U(:,3:end-2) + (DT/DX) .* (FLUX(:,1:end-1) - FLUX(:,2:end)) ;
end
