MATLAB lessons N.5¶
A solver for a 2 equation hyperbolic linear system
All MATLAB file for the lesson in a unique archive 2x2.zip
Linear2x2 Model:
File Linear2x2Model.m
:open:
1%
2% solve the linearised Hyperbolic system in one space dimension
3%
4% / \ / \ / \
5% | q1 | | 0 -a | | 0 |
6% D_t | | + D_x | | = | |
7% | q2 | | -a 0 | | 0 |
8% \ / \ / \ /
9%
10% for t > 0 and a <= x <= b
11%
12% P = parameter
13%
14% what
15% 'setup'
16% 'flux' compute flux
17%
18function [varargout] = Linear2x2Model( what, varargin )
19
20 global Linear2x2_paramaters ;
21
22 switch lower( what )
23
24 case 'setup'
25 expected( what, 1, nargin, 0, nargout ) ; % a
26 Linear2x2_paramaters = clear_struct(Linear2x2_paramaters) ;
27 Linear2x2_paramaters . a = varargin{1} ;
28
29 case 'get_dim'
30 expected( what, 0, nargin, 1, nargout ) ;
31 varargout{1} = 2 ;
32
33 case 'get_names'
34 expected( what, 0, nargin, 1, nargout ) ;
35 varargout{1} = { 'q1', 'q2' } ;
36
37 case 'flux'
38 expected( what, 1, nargin, 1, nargout ) ; % U
39 varargout{1} = Flux( varargin{1} ) ;
40
41 case 'riemann'
42 expected( what, 4, nargin, 1, nargout ) ; % UL, UR, T, X
43 UL = varargin{1} ;
44 UR = varargin{2} ;
45 T = varargin{3} ;
46 X = varargin{4} ;
47 varargout{1} = RiemannExact( UL, UR, T, X ) ;
48
49 case 'max_speed'
50 expected( what, 1, nargin, 1, nargout ) ; % U
51 varargout{1} = maxSpeed( varargin{1} ) ;
52
53 case 'extrapolate'
54 expected( what, 3, nargin, 1, nargout ) ; % U, left_bc, right_bc
55 U = varargin{1} ;
56 left_bc = varargin{2} ;
57 right_bc = varargin{3} ;
58 varargout{1} = Extrapolate( U, left_bc, right_bc ) ;
59
60 otherwise
61 error(['Linear2x2Model bad argument: ', what]) ;
62 end
63end
64
65function SE = clear_struct(S)
66 if isstruct(S)
67 SE = rmfield( S, fieldnames(S) ) ;
68 else
69 SE = S ;
70 end
71end
72
73function expected( what, ein, nin, eout, nout )
74 if ein ~= nin-1
75 error( ['in Linear2x2 call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
76 end
77 if eout ~= nout
78 error( ['in Linear2x2 call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;
79 end
80end
81
82% #######
83% # # # # # #
84% # # # # # #
85% ##### # # # ##
86% # # # # ##
87% # # # # # #
88% # ###### #### # #
89% compute the flux
90function F = Flux( US )
91 global Linear2x2_paramaters ;
92 a = Linear2x2_paramaters . a ;
93 F = [ 0, -a ; -a, 0 ] * US ;
94end
95
96% ######
97% # # # ###### # # ## # # # #
98% # # # # ## ## # # ## # ## #
99% ###### # ##### # ## # # # # # # # # #
100% # # # # # # ###### # # # # # #
101% # # # # # # # # # ## # ##
102% # # # ###### # # # # # # # #
103%
104function QS = RiemannStar( QL, QR )
105 global Linear2x2_paramaters ;
106 a = Linear2x2_paramaters . a ;
107 A = (QR + QL) ./ 2 ;
108 B = (QR - QL) ./ 2 ;
109 QS(1,:) = A(1,:) + B(2,:) ;
110 QS(2,:) = A(2,:) + B(1,:) ;
111end
112
113function US = RiemannExact( UL, UR, T, X )
114 global Linear2x2_paramaters ;
115 a = Linear2x2_paramaters . a ;
116
117 l1 = -abs(a) ;
118 l2 = +abs(a) ;
119
120 IDX = find( X - l1 * T >= 0 & X - l2 * T <= 0 ) ;
121 US(:,IDX) = RiemannStar( UL, UR ) * ones( size(IDX) ) ;
122
123 IDX = find( X - l1 * T < 0 ) ;
124 US(:,IDX) = UL * ones( size(IDX) ) ;
125
126 IDX = find( X - l2 * T > 0 ) ;
127 US(:,IDX) = UR * ones( size(IDX) ) ;
128end
129
130% #####
131% # # ## # # # # ##### ###### ###### #####
132% ## ## # # # # # # # # # # #
133% # ## # # # ## ##### # # ##### ##### # #
134% # # ###### ## # ##### # # # #
135% # # # # # # # # # # # # #
136% # # # # # # ##### # ###### ###### #####
137%
138function MS = maxSpeed( U )
139 global Linear2x2_paramaters ;
140 a = Linear2x2_paramaters . a ;
141 MS = abs(a) ;
142end
143
144% #######
145% # # # ##### ##### ## ##### #### # ## ##### ######
146% # # # # # # # # # # # # # # # # #
147% ##### ## # # # # # # # # # # # # # #####
148% # ## # ##### ###### ##### # # # ###### # #
149% # # # # # # # # # # # # # # # #
150% ####### # # # # # # # # #### ###### # # # ######
151%
152function Uextra = Extrapolate( U, L, R )
153
154 global Linear2x2_paramaters ;
155 a = Linear2x2_paramaters . a ;
156
157 % Left boundary
158 switch lower(L)
159 case 'free' % Apply transmissive boundary conditions
160 UL = U(:,1) ;
161 ULL = U(:,2) ;
162 case 'cyclic' % Apply cyclic boundary conditions
163 UL = U(:,end) ;
164 ULL = U(:,end-1) ;
165 otherwise
166 error( ['in Linear2x2 call to Extrapolate L=', L, ' unknown boundary condition']) ;
167 end
168 % Right boundary
169 switch lower(R)
170 case 'free' % Apply transmissive boundary conditions
171 UR = U(:,end) ;
172 URR = U(:,end-1) ;
173 case 'cyclic' % Apply cyclic boundary conditions
174 UR = U(:,1) ;
175 URR = U(:,2) ;
176 otherwise
177 error( ['in Linear2x2 call to Extrapolate R=', R, ' unknown boundary condition']) ;
178 end
179
180 % build vector of extrapolated values
181 Uextra = [ ULL, UL, U, UR, URR ] ;
182end
Linear2x2 Example:
File Linear2x2Example.m
:open:
1test = { 'sod' } ;
2nflux = { 'LaxFriedrichs' } ;
3
4% setup the model
5Linear2x2Model( 'setup', 1 ) ;
6
7% get names of variable
8names = Linear2x2Model( 'get_names' ) ;
9
10% Initial left and right state
11UL = [0.5 ; 2 ];
12UR = [1 ; 0.5 ];
13
14PRB . Uinit = [ UL * ones(1,50) UR * ones(1,50) ] ; % initial condition
15PRB . Xnodes = [ -1:2/100:1 ] ; % boundary cells nodes
16PRB . FinalTime = 0.1 ; % final time
17PRB . StepTime = 0.01 ; % save a frame every StepTime
18PRB . CFL = 0.9 ;
19PRB . left_bc = 'free' ; % left boundary condition type
20PRB . right_bc = 'free' ; % right boundary condition type
21
22% Exact solution of the riemann problem sampled of fine mesh
23Xe = [ -1:2/1000:1 ] ;
24Ue = Linear2x2Model( 'riemann', UL, UR, PRB . FinalTime, Xe ) ;
25
26for j=1:length(nflux)
27
28 ok = input( ['Start flux: ' nflux{j} ' [q=quit] ' ], 's' ) ;
29 if ok == 'q' ; break ; end ;
30
31 handleFig = figure('Position',[10+j*25 10+j*25 1000 800],'color','w') ;
32
33 % call the solver
34 [ UFRAMES, TFRAMES ] = Solver(PRB,@Linear2x2Model,@LaxFriedrichs) ;
35
36 % evalute cell center
37 X = (PRB . Xnodes(1:end-1)+PRB . Xnodes(2:end))/2 ;
38
39 % extract last frame
40 U = UFRAMES(:,:,end) ;
41
42 for i=1:2
43 subplot(2,1,i);
44 plot(Xe, Ue(i,:), '-b', 'LineWidth', 2, 'Color', 'green') ; hold on ;
45 plot(X, U(i,:), '-ok', 'MarkerFaceColor','r');
46 title([names{i} ' solver: ' nflux{j}], 'FontSize', 24, 'FontName', 'Arial' );
47 end
48end
Solver:
File Solver.m
:open:
1%
2% Numerical scheme for the Conservation Law
3%
4% Name of program: Solver.m
5%
6% Purpose: ........
7%
8% Programer: Enrico Bertolazzi
9% (based on the Fortran NUMERICA lib of prof. E. F. Toro )
10%
11% Last revision: 11 April 2008
12%
13% Theory is found in Chaps. 5, 13 and 14 of Reference 1
14% and in original references therein
15%
16% 1. E. F. Toro, "Riemann Solvers and Numerical Methods for Fluid Dynamics"
17% Springer-Verlag, 1997, Second Edition, 1999
18%
19% This program is inspired from the NUMERICA library:
20%
21% NUMERICA
22% A Library of Source Codes for Teaching,
23% Research and Applications,
24% by E. F. Toro
25% Published by NUMERITEK LTD,
26% Website: www.numeritek.coM
27%
28% INPUT:
29% PRB A Structure describing the problem with fields
30% Uinit = Initial solution
31% Xnodes = Boundary cell nodes
32% StepTime = solution is saved every StepTime
33% FinalTime = solution is computed up to FinalTime
34% CFL = advancing time step CFL
35%
36% MODEL A function with variable arguments.
37% First argument simulate method name in Object Oriented Paradigm
38% MODEL must support the following call:
39% Uextr = MODEL( 'extrapolate', U, left_bc, right_bc )
40% Extrapolate 4 cell of states(2 on left 2 on right) for boudary condition
41% S = MODEL( 'max_speed', U )
42% Evaluate maximum speed signal of the flux for CFL computation
43% F = MODEL( 'flux', U )
44% Ealuate model flux
45% OUTPUT:
46% UFRAMES = matrix with the frames of the computed solution.
47% UFRAMES(:,:,k) is the solution at time level k-th
48% TFRAMES = vector with the time frames. TFRAMES(i) the time of the computed frame UFRAMES(i,:).
49%
50function [ UFRAMES, TFRAMES ] = Solver(PRB,MODEL,NFLUX)
51
52 % Initial conditions are set up
53 U = PRB . Uinit ;
54
55 % Compute cell size
56 DX = PRB . Xnodes(2) - PRB . Xnodes(1) ; % equal spaced nodes!!
57
58 % save initial frame
59 UFRAMES(:,:,1) = U ;
60 TFRAMES(1) = 0 ;
61
62 kframe = 1 ; % number of saved framed
63 Time = 0 ; % initial time
64 NextTime = PRB.StepTime ; % Time to save next frame
65
66 % Time marching procedure
67 for NITER=1:10000
68
69 % Extrapolate for Boundary conditions
70 Uextra = feval( MODEL, 'extrapolate', U, PRB.left_bc, PRB.right_bc ) ;
71
72 % Courant-Friedrichs-Lewy (CFL) condition imposed
73 DT = DX * PRB.CFL / feval( MODEL, 'max_speed', U ) ;
74
75 % update time
76 Time = Time + DT ;
77
78 % adjust time step if necessary to match the time grid
79 while Time >= NextTime
80 % reduce DT in order to match NextTime
81 DT1 = DT - ( Time - NextTime ) ;
82
83 % Compute intercell flux FLUX
84 FLUX = feval( NFLUX, MODEL, DT1, DX, Uextra ) ;
85
86 kframe = kframe+1 ;
87 % Solution is updated according to conservative formula and saved in the frame
88 UFRAMES(:,:,kframe) = U + (DT/DX) * (FLUX(:,1:end-1) - FLUX(:,2:end)) ;
89 TFRAMES(kframe) = NextTime ;
90
91 % update next time frame
92 NextTime = NextTime + PRB.StepTime ;
93 end
94
95 % Compute intercell flux FLUX
96 FLUX = feval( NFLUX, MODEL, DT, DX, Uextra ) ;
97 % Solution is updated according to conservative formula
98 U = U + (DT/DX) * ( FLUX(:,1:end-1) - FLUX(:,2:end) ) ;
99
100 if Time >= PRB.FinalTime ; break ; end
101 end
102
103 % save last frame
104 kframe = kframe+1 ;
105 UFRAMES(:,:,kframe) = U ;
106 TFRAMES(kframe) = Time ;
107end
LaxFriedrichs Numerical Fluxes:
File LaxFriedrichs.m
:open:
1% _ _____ _ _ _ _
2% | | __ ___ __ | ___| __(_) ___ __| |_ __(_) ___| |__ ___
3% | | / _` \ \/ /____| |_ | '__| |/ _ \/ _` | '__| |/ __| '_ \/ __|
4% | |__| (_| |> <_____| _|| | | | __/ (_| | | | | (__| | | \__ \
5% |_____\__,_/_/\_\ |_| |_| |_|\___|\__,_|_| |_|\___|_| |_|___/
6%
7% Compute the Lax-Friedrichs flux
8%
9function Flux = LaxFriedrichs( MODEL, DT, DX, US )
10 FS = feval( MODEL, 'flux', US ) ;
11 %
12 UL = US(:,2:end-2) ;
13 UR = US(:,3:end-1) ;
14 FL = FS(:,2:end-2) ;
15 FR = FS(:,3:end-1) ;
16 %
17 Flux = 0.5 .* ( FL + FR + (DT/DX) * (UL-UR) ) ;
18end