Lezione del 10 maggio 2007 - Enrico Bertolazzi

:open:

%
%     A selection of 8 schemes for the linear advection equation
%
%     Name of program: LagenSolver
%
%     Purpose: to solve the linear advection equation with constant
%              coefficient by a selection of 8 schemes, namely:
%
%              The Godunov first-order upwind scheme
%              The Toro-Billett first-order upwind scheme
%              The Lax-Friedrichs scheme (first-order centred)
%              The FORCE scheme (first-order centred)
%              The Godunov first-order centred scheme
%              The Lax-Wendroff scheme (second order, oscillatory)
%              The Fromm scheme (second order, oscillatory)
%              The Warming-Beam scheme (second order, oscillatory)
%
%     Programer: Enrico Bertolazzi 
%               (based on the Fortran NUMERICA lib of prof. E. F. Toro )
%
%     Last revision: 10st May 2007
%
%     Theory is found in Chaps. 5, 7 and 13 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 transaled in MATLAB 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
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  PURPOSE:
%
%  solve the linear advection problem
%
%     u (t,x)  + a u (t,x) = 0     t > 0, -1 <= x <= 1
%      t            x
%
%     u(0,x) = u0(x)              -1 <= x <= 1
%
%  INPUT:
%     PRB A Structure describing the problem with fields
%     problem   = an integer selectiong the numerical flux used
%               1 - The Godunov first-order upwind scheme
%               2 - The Toro-Billett first-order upwind scheme
%               3 - The Lax-Friedrichs scheme (first-order centred)
%               4 - The FORCE scheme (first-order centred)
%               5 - The Godunov first-order centred scheme
%               6 - The Lax-Wendroff scheme (second order, oscillatory)
%               7 - The Fromm scheme  (second order, oscillatory)
%               8 - The Warming-Beam scheme (second order, oscillatory)
%     init      = an integer selectiong the initial data used
%               1 - square wave 
%               2 - smooth wave 
%     Ncells    = number of cell discretizing [-1,1] the domain of the problems
%     CFL       = time step advancing CFL ratio
%     a         = speed of the signal of the PDE
%     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(i,:) is the solution
%                 at time level i-th
%     TFRAMES   = vector with the time frames. TFRAMES(i) the time of the computed frame UFRAMES(i,:).
%
%  Example of use
%
%  PRB . problem   = 1   % select Godunov first-order upwind scheme
%  PRB . init      = 1   % select square wave for initial data
%  PRB . Ncells    = 100 % select 100 cell inteval subdivision
%  PRB . CFL       = 0.9 % advance with DT do not excede 0.9 * CFL
%  PRB . a         = 0.1 % speed of teh signal
%  PRB . StepTime  = 0.5 % save the frame every 0.5 second
%  PRB . FinalTime = 3   % compute the solution op to 3 second
%
%  call the solver
%  [X,UFRAMES,TFRAMES] = LagenSolver(PRB) ;
%
%  plot the solution al frame n. 30
%  plot( X, UFRAMES(30,:) );
%
function [X,UFRAMES,TFRAMES] = LagenSolver(PRB)
  N = PRB.Ncells ;
  % Initial conditions are set up
  switch( PRB . init )
    case 1
      [X,U,DX] = InitSquare(N) ;
    case 2
      [X,U,DX] = InitSmooth(N) ;
  end
  % maximum speed
  SMAX = abs(PRB.a) ;

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

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

  % Time marching procedure
  for i=1:1000

    % Boundary conditions are cyclic so we set ciclycally the vectors
    ULL = [ U(N-1), U(N), U(1:N-1) ] ;
    UL  = [ U(N), U ] ;
    UR  = [ U, U(1) ] ;
    URR = [ U(2:N), U(1), U(2) ] ;

    % Courant-Friedrichs-Lewy (CFL) condition imposed
    DT   = PRB.CFL*DX/SMAX ;
    Time = Time + DT ;

    % adjust time step if necessary to match the time grid
    if Time >= NextTime
      % reduce DT in order to match NextTime
      DT       = DT - ( Time - NextTime ) ;
      Time     = NextTime ;
      NextTime = NextTime + PRB.StepTime ;
      k        = -k ; % mark solution to be saved
    end

    % Intercell numerical fluxes are computed.
    switch( PRB . problem )
    case 0
      FLUX = Centered        (UL, UR, PRB.a, DT, DX ) ;
    case 1
      FLUX = Godunov         (UL, UR, PRB.a, DT, DX ) ;
    case 2
      FLUX = ToroBillet      (ULL, UL, UR, URR, PRB.a, DT, DX ) ;
    case 3
      FLUX = LaxFriedrichs   (UL, UR, PRB.a, DT, DX ) ;
    case 4
      FLUX = Force           (UL, UR, PRB.a, DT, DX ) ;
    case 5
      FLUX = GodunovCentered (UL, UR, PRB.a, DT, DX ) ;
    case 6
      FLUX = Richtmyer       (UL, UR, PRB.a, DT, DX ) ;
    case 7
      FLUX = Fromm           (ULL, UL, UR, URR, PRB.a, DT, DX ) ;
    case 8
      FLUX = WarmingBeam     (ULL, UL, UR, URR, PRB.a, DT, DX ) ;
    end
    % Solution is updated according to conservative formula
    U = U + (DT/DX) .* (FLUX(1:N) - FLUX(2:N+1)) ;

    % if k < 0 solution must be saved
    if k < 0 
      k            = -k+1 ;
      UFRAMES(k,:) = U ;
      TFRAMES(k)   = Time ;
    end

    if Time >= PRB.FinalTime
      break ;
    end
  end
end

%
%  Purpose: to set initial conditions for solution U
%
function [X,U,DX] = InitSmooth(N)
  % Calculate mesh size DX
  DX = 2.0/N ;
  %  Initialise arrays (index 1,2 and IDIM+3, IDIM+4 are used for "ghost cells")
  % smooth profile
  X  = [-1+DX/2:DX:1-DX/2] ;
  U  = exp( -8 .* X .^ 2 ) ;
end
%
function [X,U,DX] = InitSquare(N)
  % Calculate mesh size DX
  DX = 2.0/N ;
  N3 = floor(N/3) ;
  M  = N - N3 - N3 ;
  % Initialise arrays
  X  = [-1+DX/2:DX:1-DX/2] ;
  U  = [zeros(1,N3) ones(1,M) zeros(1,N3) ] ;
end

%
%
%  positioning of the variables
%
%     U(1)   U(2)   U(3)                             U(N)
%  +------+------+------+------+------+-           +------+
%  |      |      |      |      |      |   ...      |      |
%  F(1)  F(2)  F(3)                               F(N)  F(N+1)
%
%  U(i) is the cell average allocated to center of the cell i.
%  F(i) is the flux at the interface between cells i end cell i-1.


% Purpose: to compute intercell fluxes as left and right average: CANNOT WORK
%
function F = Centered(UL,UR,a,DT,DX)
  F = (a/2) .* (UL+UR) ;
end
% Purpose: to compute intercell fluxes according to the
%          Godunov first order upwind method
%
function F = Godunov(UL,UR,a,DT,DX)
  % the flux at interface i is computed
  if a > 0
    % using the left state
    F = a .* UL ;
  else
    % using the right state
    F = a .* UR ;
  end
end
%
% Purpose: to compute intercell fluxes according to the
%          Toro-Billett first order CFL-2 scheme
%          See Eq. 13.22, Chapter 13, Ref. 1
function [FLUX] = ToroBillet(ULL,UL,UR,URR,a,DT,DX)
  CFLNUM = a*DT/DX ;
  if ( abs(CFLNUM) <= 1 )
    % Conventional Godunov upwind flux
    FLUX = Godunov(UL,UR,a,DT,DX) ;
  else
    if a > 0
      % Information travels from the left
      cfll =  ((CFLNUM - 1)/CFLNUM) * a ;
      cfl  =  (1.0/CFLNUM) * a ;
      FLUX = cfll * ULL + cfl * UL ;
    else
      % Information travels from the right
      cfrr =  ((CFLNUM + 1)/CFLNUM) * a ;
      cfr  =  -(1.0/CFLNUM) * a ;
      FLUX = cfl * UL + cfr * UR ;
    end
  end
end

% Purpose: to compute intercell fluxes according to the
%          Lax-Friedrichs method
%          See Eq. 5.77, Chapter 5, Ref. 1
function [FLUX] = LaxFriedrichs(UL,UR,a,DT,DX)
  CFLNUM = a*DT/DX ;
  FLUX   = 0.5 .* ( a .* (UL+UR) + (DX/DT) * (UL-UR) ) ;  
end

% Purpose: to compute intercell fluxes according to the
%          Force method
%          See Eq. 7.32, Chapter 7, Ref. 1
function [FLUX] = Force(UL,UR,a,DT,DX)
  CFLNUM = a*DT/DX ;
  % Compute Lax-Friedrichs flux
  FluxLF = 0.5 .* ( a .* (UL+UR) + (DX/DT) * ( UL-UR ) ) ;
  % Compute Richtmyer flux
  FluxR  = 0.5 .* a .* ( (UL + UR) + (DT/DX) .* a .* ( UL - UR ) ) ;
  % Compute Force flux
  FLUX = 0.5*(FluxLF + FluxR) ;
end

% Purpose: to compute intercell fluxes according to the
%          Godunov first order upwind method
%          See Eq. 13.73, Chapter 13, Ref. 1
function FLUX = GodunovCentered(UL,UR,a,DT,DX)
  CFLNUM = a*DT/DX ;
  FLUX   = 0.5 .* a .* ( (1+2*CFLNUM) * UL + ( 1-2*CFLNUM) * UR ) ;
end

% Purpose: to compute intercell fluxes according to the
%          Richtmyer method, or two step Lax-Wendroff method
%          See Eq. 5.79, Chapter 5, Ref. 1
function FLUX = Richtmyer(UL,UR,a,DT,DX)
  FLUX = 0.5 .* a .* ( (UL + UR) + (DT/DX) .* a .* ( UL - UR ) ) ;
end

% Purpose: to compute intercell fluxes according to the
%          Fromm first order upwind method
%          See Eq. 13.35, Chapter 13, Ref. 1
function FLUX = Fromm(ULL,UL,UR,URR,a,DT,DX)
  CFLNUM = a*DT/DX ;
  if ( a > 0 ) 
    FLUX = a .* (UL + 0.25 .* (1-CFLNUM) .* (UR - ULL)) ;
  else
    FLUX = a .* (UR - 0.25 .* (1+CFLNUM) .* (URR - UR)) ;
  end
end

% Purpose: to compute intercell fluxes according to the
%          Warming-Beam method
%          See Eq. 13.21, Chapter 13, Ref. 1
function FLUX = WarmingBeam(ULL,UL,UR,URR,a,DT,DX)
  CFLNUM = a*DT/DX ;
  if ( a > 0 ) 
    FLUX = 0.5 .* a .* ( (3-CFLNUM) .* UL - (1-CFLNUM) .* ULL );
  else
    FLUX = 0.5 .* a .* ( (3+CFLNUM) .* UR - (1+CFLNUM) .* URR );
  end
end

:open:

PRB . problem   = 8    ; % select Godunov first-order upwind scheme
PRB . init      = 1    ; % select square wave for initial data
PRB . Ncells    = 100  ; % select 100 cell inteval subdivision
PRB . CFL       = 0.49 ; % advance with DT do not excede 0.9 * CFL
PRB . a         = 0.25 ; % speed of the signal
PRB . StepTime  = 1    ; % save the frame every 0.5 second
PRB . FinalTime = 200  ; % compute the solution op to 3 second
% call the solver
[X,UFRAMES,TFRAMES] = LagenSolver(PRB) ;
% generate movie of the solution
LagenMovie(X,UFRAMES,TFRAMES) ;

% plot the solution al frame n. 4
%[,N] = size(TFRAMES) ; 
%plot( X, UFRAMES(N,:) );

:open:

function aviobj = LagenMovie(X,UFRAMES,TFRAMES)
  [dummy,nframes] = size(TFRAMES) ;
  aviobj = avifile('example.avi')
  rect   = [100 100 612 612] ;
  fig    = figure('position',rect) ;
  set(fig,'DoubleBuffer','on');
  axis manual ;
  title('Lagen Movie');
  for i=1:nframes
    clc;
    plot(X, UFRAMES(i,:),'-ob','MarkerFaceColor','r');
    axis([-1 1 -0.5 1.5]);
    aviobj = addframe(aviobj,getframe(fig));
  end
  aviobj = close(aviobj);
end

:open:

%
%
function [X,U] = LagenExactSquare(N,TIME,SPEEDA)
  % Calculate mesh size DX
  DX = 2.0/N ;
  X  = [-1+DX/2-TIME*SPEEDA:DX:1-DX/2-TIME*SPEEDA] ;
  % put the solution in [-1..1]
  U = zeros(1,N) ;
  for i=1:N
    while ( X(i) < -1 )
      X(i) = X(i) + 2 ;
    end ;
    while ( X(i) >  1 )
      X(i) = X(i) - 2 ;
    end ;
    if X(i) >= -1/3.0 && X(i) <= 1/3.0 
      U(i) = 1 ;
    end
  end
  % smooth profile
  X = [-1+DX/2:DX:1-DX/2] ;
end

:open:

%
%  Purpose: to set initial conditions for solution U and
%           initialise other variables.
%
%  Variables:
%
%  U  Array for numerical solution
%
function [X,U] = LagenExactSmooth(N,TIME,SPEEDA)
  % Calculate mesh size DX
  DX = 2.0/N ;
  X = [-1+DX/2-TIME*SPEEDA:DX:1-DX/2-TIME*SPEEDA] ;
  % put the solution in [-1..1]
  for i=1:N
    while ( X(i) < -1 )
      X(i) = X(i) + 2 ;
    end ;
    while ( X(i) >  1 )
      X(i) = X(i) - 2 ;
    end ;
  end
  % smooth profile
  U = exp( -8 .* X .^ 2 ) ;
  X = [-1+DX/2:DX:1-DX/2] ;
end

Lezione del 10 maggio 2007¶