MATLAB lessons N.4¶
A library to solve Bughers equation with various methods
All MATLAB file for the lesson in a unique archive Burgers.zip
Burgers Model:
File BurgersModel.m
:open:
1% Model for the Burgers equations
2%
3% u(t,x) + u(t,x) * u (t,x) = 0
4% t x
5%
6% for t > 0 and -1 <= x <= 1
7%
8function [varargout] = BurgerModel( what, varargin )
9
10 % BurgerModel_paramaters
11 % test = which test is used
12 % nflux = which numerical flux is used
13 % bc_left = left bounday condition
14 % bc_right = right boundaruy condition
15
16 global BurgerModel_paramaters ;
17
18 switch lower(what)
19
20 case 'setup'
21 expected( what, 0, nargin, 0, nargout ) ;
22 BurgerModel_paramaters = clear_struct(BurgerModel_paramaters) ;
23
24 case 'set_test'
25 expected( what, 1, nargin, 0, nargout ) ; % TEST
26 ChooseTest( varargin{1} ) ;
27
28 case 'set_limiter'
29 expected( what, 1, nargin, 0, nargout ) ;
30 BurgerModel_paramaters . limiter = lower( varargin{1} ) ;
31
32 case 'set_numerical_flux'
33 expected( what, 1, nargin, 0, nargout ) ;
34 BurgerModel_paramaters . nflux = lower( varargin{1} ) ;
35
36 case 'get_dim'
37 expected( what, 0, nargin, 1, nargout ) ;
38 varargout{1} = 1 ;
39
40 case 'get_dx'
41 expected( what, 0, nargin, 1, nargout ) ;
42 varargout{1} = BurgerModel_paramaters . DX ;
43
44 case 'get_x'
45 expected( what, 0, nargin, 1, nargout ) ;
46 varargout{1} = BurgerModel_paramaters . X ;
47
48 case 'get_final_time'
49 expected( what, 0, nargin, 1, nargout ) ;
50 varargout{1} = BurgerModel_paramaters . tf ;
51
52 case 'get_cfl'
53 expected( what, 0, nargin, 1, nargout ) ;
54 varargout{1} = MaxCFL() ;
55
56 case 'names'
57 expected( what, 0, nargin, 1, nargout ) ;
58 varargout{1} = 'u' ;
59
60 case 'flux'
61 expected( what, 1, nargin, 1, nargout ) ; % U
62 varargout{1} = Flux( varargin{1} ) ;
63
64 case 'init'
65 expected( what, 1, nargin, 1, nargout ) ; % N
66 varargout{1} = Init( varargin{1} ) ;
67
68 case 'muscle_nflux'
69 expected( what, 2, nargin, 1, nargout ) ; % UL, UR
70 varargout{1} = Godunov( varargin{1}, varargin{2} ) ;
71
72 case 'conservative_to_physical'
73 varargout{1} = varargin{1} ;
74
75 case 'physical_to_conservative'
76 varargout{1} = varargin{1} ;
77
78 case 'conservative_to_characteristic'
79 varargout{1} = varargin{1} ;
80
81 case 'characteristic_to_conservative'
82 varargout{1} = varargin{1} ;
83
84 case 'dt_from_cfl'
85 expected( what, 2, nargin, 1, nargout ) ; % U, CFLcoeff
86 varargout{1} = DT_from_CFL( varargin{1}, varargin{2} ) ;
87
88 case 'extrapolate_for_bc'
89 expected( what, 1, nargin, 4, nargout ) ; % U
90 [ varargout{1}, varargout{2}, varargout{3}, varargout{4} ] = Extrapolate_for_BC( varargin{1} ) ;
91
92 case 'nflux'
93 expected( what, 2, nargin, 1, nargout ) ; % DT, U
94 varargout{1} = NFlux( varargin{1}, varargin{2} ) ;
95
96 otherwise
97 error(['BurgerModel bad argument: ', what]) ;
98 end
99end
100
101function SE = clear_struct(S)
102 if isstruct(S)
103 SE = rmfield( S, fieldnames(S) ) ;
104 else
105 SE = S ;
106 end
107end
108
109function expected( what, ein, nin, eout, nout )
110 if ein ~= nin-1
111 error( ['in BurgerModel call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
112 end
113 if eout ~= nout
114 error( ['in BurgerModel call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;
115 end
116end
117
118% #######
119% # # # # # #
120% # # # # # #
121% ##### # # # ##
122% # # # # ##
123% # # # # # #
124% # ###### #### # #
125
126% compute the flux
127function F = Flux( U )
128 F = 0.5 * U .^ 2 ;
129end
130
131% ###
132% # # # # #####
133% # ## # # #
134% # # # # # #
135% # # # # # #
136% # # ## # #
137% ### # # # #
138
139function U = Init(N)
140 global BurgerModel_paramaters ;
141 % Calculate mesh size DX
142 DX = 2.0/N ;
143 X = [-1+DX/2:DX:1-DX/2] ;
144
145 if isfield(BurgerModel_paramaters, 'test')
146 switch BurgerModel_paramaters . test
147 case 'classroomproblem'
148 DX = 8.0/N ;
149 X = [0+DX/2:DX:8-DX/2] ; % center of the cell
150 N1 = floor(3*N/8) ;
151 N2 = N - N1 ;
152 U = [ones(1,N1) zeros(1,N2) ] ;
153 case 'square'
154 N3 = floor(N/3) ;
155 M = N - N3 - N3 ;
156 U = [zeros(1,N3) ones(1,M) zeros(1,N3) ] ;
157 case 'smooth'
158 U = exp( -8 .* X .^ 2 ) ;
159 otherwise
160 error(['BurgerModel::Init unknown test ', BurgerModel_paramaters . test]) ;
161 end
162 else
163 error('BurgerModel::Init unknown test you must set it by set_test!') ;
164 end
165 BurgerModel_paramaters . DX = DX ;
166 BurgerModel_paramaters . X = X ;
167end
168
169
170% ##### #######
171% # # # # #### #### #### ###### # ###### #### #####
172% # # # # # # # # # # # # #
173% # ###### # # # # #### ##### # ##### #### #
174% # # # # # # # # # # # # #
175% # # # # # # # # # # # # # # # #
176% ##### # # #### #### #### ###### # ###### #### #
177%
178function ChooseTest(test)
179
180 global BurgerModel_paramaters ;
181
182 switch test
183 case 'classroomproblem'
184 tf = 3.5 ;
185 bc_left = 'free' ;
186 bc_right = 'free' ;
187 case 'square'
188 tf = 3 ;
189 bc_left = 'cyclic' ;
190 bc_right = 'cyclic' ;
191 case 'smooth'
192 tf = 4 ;
193 bc_left = 'cyclic' ;
194 bc_right = 'cyclic' ;
195 otherwise
196 error(['BurgerModel:ChooseTest test name',test] ) ;
197 end
198
199 BurgerModel_paramaters . test = test ;
200 BurgerModel_paramaters . bc_left = bc_left ;
201 BurgerModel_paramaters . bc_right = bc_right ;
202 BurgerModel_paramaters . tf = tf ;
203
204end
205
206% ###### ####### ##### ####### #
207% # # # ###### ##### #### # # # # # #
208% # # # # # # # # ## ## # # #
209% # # # ##### # # # # # ## # # ##### #
210% # # # # ##### # # # # # # #
211% # # # # # # # # # # # # # #
212% ###### # # # # #### # # ##### # #######
213% ####### #######
214% Compute DT from CFL
215function DT = DT_from_CFL( U, CFLcoeff )
216 global BurgerModel_paramaters ;
217 DX = BurgerModel_paramaters . DX ;
218 DT = CFLcoeff * DX / max(abs(U)) ;
219end
220
221% #######
222% # # # ##### ##### ## ##### #### # ## ##### ######
223% # # # # # # # # # # # # # # # # #
224% ##### ## # # # # # # # # # # # # # #####
225% # ## # ##### ###### ##### # # # ###### # #
226% # # # # # # # # # # # # # # # #
227% ####### # # # # # # # # #### ###### # # # ######
228%
229% ###### #####
230% ###### #### ##### # # # #
231% # # # # # # # #
232% ##### # # # # ###### #
233% # # # ##### # # #
234% # # # # # # # # #
235% # #### # # ###### #####
236% #######
237
238% Extrapolate values for boundary condition
239function [ULL, UL, UR, URR] = Extrapolate_for_BC( U )
240 global BurgerModel_paramaters ;
241 L = BurgerModel_paramaters . bc_left ;
242 R = BurgerModel_paramaters . bc_right ;
243 % Left boundary
244 switch L
245 case 'free' % Apply transmissive boundary conditions
246 UL = U(:,1) ;
247 ULL = U(:,2) ;
248 case 'cyclic' % Apply cyclic boundary conditions
249 UL = U(:,end) ;
250 ULL = U(:,end-1) ;
251 otherwise % Apply reflective boundary conditions
252 UL = [ U(1,1) ; -U(2,1) ] ;
253 ULL = [ U(1,2) ; -U(2,2) ] ;
254 end
255 % Right boundary
256 switch R
257 case 'free' % Apply transmissive boundary conditions
258 UR = U(:,end) ;
259 URR = U(:,end-1) ;
260 case 'cyclic' % Apply cyclic boundary conditions
261 UR = U(:,1) ;
262 URR = U(:,2) ;
263 otherwise % Apply reflective boundary conditions
264 UR = [ U(1,end) ; -U(2,end) ] ;
265 URR = [ U(1,end-1) ; -U(2,end-1) ] ;
266 end
267end
268
269% ####### ##### ####### #
270% # # ##### ##### # # # ## # # # # #
271% # # # # # # ## ## # # # # # #
272% # # # # # # # ## # # # # # ##### #
273% # # ##### # # # # ###### # # # #
274% # # # # # # # # # # # # # #
275% ####### # # # # # # # ###### ##### # #######
276
277function CFL = MaxCFL
278 global BurgerModel_paramaters ;
279 nflux = BurgerModel_paramaters . nflux ;
280 %
281 switch nflux
282 case 'godunov'
283 CFL = 1.0 ;
284 otherwise
285 CFL = StandardFlux( 'get_cfl', nflux ) ;
286 end
287end
288
289
290% # # # # # # ###### ##### # #### ## #
291% ## # # # ## ## # # # # # # # # #
292% # # # # # # ## # ##### # # # # # # #
293% # # # # # # # # ##### # # ###### #
294% # ## # # # # # # # # # # # # #
295% # # #### # # ###### # # # #### # # ######
296%
297% ###### # # # # # ###### ####
298% # # # # # # # #
299% ##### # # # ## ##### ####
300% # # # # ## # #
301% # # # # # # # # #
302% # ###### #### # # ###### ####
303
304function FLUX = NFlux( DT, U )
305
306 global BurgerModel_paramaters ;
307
308 ULL = U(:,1:end-3) ;
309 UL = U(:,2:end-2) ;
310 UR = U(:,3:end-1) ;
311 URR = U(:,4:end-0) ;
312
313 nflux = BurgerModel_paramaters . nflux ;
314 DX = BurgerModel_paramaters . DX ;
315
316 switch nflux
317 % monolitic numerical fluxes
318 case 'godunov'
319 FLUX = Godunov(UL,UR) ;
320 otherwise
321 FLUX = StandardFlux( 'nflux', nflux, @Flux, DT, DX, U ) ;
322 end
323end
324
325% _ ___ _ __ ___ _ _ _ _ _
326% |_ | |_) (_ | / \ |_) | \ |_ |_)
327% | _|_ | \ __) | \_/ | \ |_/ |_ | \
328%
329
330% #####
331% # # #### ##### # # # # #### # #
332% # # # # # # # ## # # # # #
333% # #### # # # # # # # # # # # # #
334% # # # # # # # # # # # # # # #
335% # # # # # # # # # ## # # # #
336% ##### #### ##### #### # # #### ## %
337%
338% Compute the solution of the Riemann problem
339%
340function F = Godunov( UL, UR )
341
342 %%USTAR = zeros(size(UL)) ;
343
344 IDX = find( UL > UR ) ;
345 % Solution is a shock wave Compute shock speed S
346 S(IDX) = 0.5*(UL(IDX) + UR(IDX)) ;
347 % Sample the state along the t-axis
348 II = IDX( find( S(IDX) >= 0 )) ;
349 USTAR(II) = UL(II) ;
350 II = IDX( find( S(IDX) < 0 )) ;
351 USTAR(II) = UR(II) ;
352
353 IDX = find( UL <= UR ) ;
354 % Solution is a rarefaction wave. There are 3 cases:
355 II = IDX(find( UL(IDX) >= 0 )) ;
356 % Right supersonic rarefaction
357 USTAR(II) = UL(II) ;
358
359 II = IDX(find( UR(IDX) <= 0 )) ;
360 % Left supersonic rarefaction
361 USTAR(II) = UR(II) ;
362
363 II = IDX(find( UL(IDX) <= 0 & UR(IDX) >= 0 )) ;
364 % Transonic rarefaction
365 USTAR(II) = 0 ;
366
367 F = 0.5 * USTAR .^ 2 ;
368
369end
Driver to test the Solver:
File BurgersExample.m
:open:
1test = { 'square', 'smooth' } ;
2
3nflux = { 'Godunov', ...
4 'force', ...
5 'LaxFriedrichs', ...
6 'Richtmyer', ...
7 'GodunovCentered' } ;
8
9close all ;
10
11for j=1:length(nflux)
12
13 ok = input( ['Start numerical flux: ' nflux{j} ' [q=quit] ' ], 's' ) ;
14 if ok == 'q'
15 break ;
16 end ;
17
18 handleFig = figure('Position',[10 10 1000 600],'color','w') ;
19
20 for i=1:length(test)
21
22 BurgersModel( 'setup' ) ;
23 BurgersModel( 'set_test', test{i} ) ;
24 BurgersModel( 'set_numerical_flux', nflux{j} ) ;
25 %%BurgersModel( 'set_limiter', 'superbee' ) ;
26
27 PRB . CFL = 0.9*BurgersModel( 'get_cfl' ) ;
28 PRB . FinalTime = BurgersModel( 'get_final_time' ) ;
29 PRB . Ncells = 100 ; % select 100 cell inteval subdivision
30 PRB . StepTime = PRB . FinalTime / 100 ;
31
32 % call the solver
33 [ UFRAMES, TFRAMES ] = Solver(PRB,@BurgersModel) ;
34
35 U = UFRAMES(:,:,end) ;
36 X = BurgersModel( 'get_x' ) ;
37
38 subplot(length(test),1,i);
39 plot(X,U,'-ob','MarkerFaceColor','b');
40 title(['test case: ' test{i} ' solver: ' nflux{j}], ...
41 'FontSize', 24, 'FontName', 'Arial' );
42 axis([X(1) X(end) -0.2 1.2]) ;
43 end
44end
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% Ncells = number of cell discretizing [-1,1] the domain of the problems
31% CFL = time step advancing CFL ratio
32% StepTime = solution is saved every StepTime
33% FinalTime = solution is computed up to FinalTime
34%
35% OUTPUT:
36% X = the cell boundary of the subdivided interval [-1,1]; X(1) = -1,..., X(N) = 1
37% UFRAMES = matrix with the frames of the computed solution.
38% UFRAMES(:,:,k) is the solution at time level k-th
39% TFRAMES = vector with the time frames. TFRAMES(i) the time of the computed frame UFRAMES(i,:).
40%
41function [ UFRAMES, TFRAMES ] = Solver(PRB,MODEL)
42
43 % Initial conditions are set up
44 U = feval( MODEL, 'init', PRB.Ncells ) ;
45 DX = feval( MODEL, 'get_dx' ) ;
46
47 kframe = 1 ;
48 Time = 0 ;
49 NextTime = PRB.StepTime ;
50
51 % save initial frame
52 UFRAMES(:,:,kframe) = U ;
53 TFRAMES(kframe) = Time ;
54
55 % Time marching procedure
56 for NITER=1:10000
57
58 % Extrapolate for Boundary conditions
59 [ ULL, UL, UR, URR ] = feval( MODEL, 'extrapolate_for_bc', U ) ;
60
61 Uold = [ ULL, UL, U, UR, URR ] ;
62
63 % Courant-Friedrichs-Lewy (CFL) condition imposed
64 DT = feval( MODEL, 'dt_from_cfl', U, PRB.CFL ) ;
65 Time = Time + DT ;
66
67 % adjust time step if necessary to match the time grid
68 while Time >= NextTime
69 % reduce DT in order to match NextTime
70 DT1 = DT - ( Time - NextTime ) ;
71 Time1 = NextTime ;
72
73 U1 = Update(MODEL, DT1, DX, Uold ) ;
74
75 kframe = kframe+1 ;
76 UFRAMES(:,:,kframe) = U1 ;
77 TFRAMES(kframe) = Time1 ;
78
79 NextTime = NextTime + PRB.StepTime ;
80
81 end
82
83 % Solution is updated according to conservative formula
84 U = Update(MODEL, DT, DX, Uold ) ;
85
86 if Time >= PRB.FinalTime ; break ; end
87 end
88 kframe = kframe+1 ;
89 UFRAMES(:,:,kframe) = U ;
90 TFRAMES(kframe) = Time ;
91end
92
93function Unew = Update( MODEL, DT, DX, U )
94 % Compute intercell flux FLUX
95 FLUX = feval( MODEL, 'nflux', DT, U ) ;
96 % Solution is updated according to conservative formula
97 Unew = U(:,3:end-2) + (DT/DX) .* (FLUX(:,1:end-1) - FLUX(:,2:end)) ;
98end
Standard Numerical Fluxes:
File StandardFlux.m
:open:
1% # # # # # # ###### ##### # #### ## #
2% ## # # # ## ## # # # # # # # # #
3% # # # # # # ## # ##### # # # # # # #
4% # # # # # # # # ##### # # ###### #
5% # ## # # # # # # # # # # # # #
6% # # #### # # ###### # # # #### # # ######
7%
8% ###### # # # # # ###### ####
9% # # # # # # # #
10% ##### # # # ## ##### ####
11% # # # # ## # #
12% # # # # # # # # #
13% # ###### #### # # ###### ####
14
15function [varargout] = StandardFlux( what, varargin )
16
17 numerical_flux_name = lower(varargin{1}) ;
18
19 switch what
20 case 'nflux'
21 Flux = varargin{2} ;
22 DT = varargin{3} ;
23 DX = varargin{4} ;
24 U = varargin{5} ;
25 %
26 switch numerical_flux_name
27 case 'laxfriedrichs'
28 varargout{1} = LaxFriedrichs( Flux, DT, DX, U ) ;
29 case 'force'
30 varargout{1} = Force( Flux, DT, DX, U ) ;
31 case 'godunovcentered'
32 varargout{1} = GodunovCentered( Flux, DT, DX, U ) ;
33 case 'richtmyer'
34 varargout{1} = Richtmyer( Flux, DT, DX, U ) ;
35 otherwise
36 error( [ 'StandardFlux nflux: ', varargin{1}, ' unknown' ] ) ;
37 end
38
39 case 'get_cfl'
40 switch numerical_flux_name
41 case { 'laxfriedrichs', 'force', 'richtmyer' }
42 varargout{1} = 1.0 ;
43 case { 'slic', 'godunovcentered' }
44 varargout{1} = 0.5 ;
45 otherwise
46 error( [ 'StandardFlux nflux: ', numerical_flux_name, ' unknown' ] ) ;
47 end
48 end
49end
50% _ _____ _ _ _ _
51% | | __ ___ __ | ___| __(_) ___ __| |_ __(_) ___| |__ ___
52% | | / _` \ \/ /____| |_ | '__| |/ _ \/ _` | '__| |/ __| '_ \/ __|
53% | |__| (_| |> <_____| _|| | | | __/ (_| | | | | (__| | | \__ \
54% |_____\__,_/_/\_\ |_| |_| |_|\___|\__,_|_| |_|\___|_| |_|___/
55%
56% Compute the Lax-Friedrichs flux
57%
58function F = LaxFriedrichs( FLUX, DT, DX, U )
59 UL = U(:,2:end-2) ;
60 UR = U(:,3:end-1) ;
61 %%
62 FL = feval( FLUX, UL) ;
63 FR = feval( FLUX, UR) ;
64 F = 0.5 .* ( FL + FR + (DX/DT) .* (UL-UR) ) ;
65end
66% _____
67% | ___|__ _ __ ___ ___
68% | |_ / _ \| '__/ __/ _ \
69% | _| (_) | | | (_| __/
70% |_| \___/|_| \___\___|
71%
72% compute an intercell flux according to the FORCE scheme.
73%
74function F = Force( FLUX, DT, DX, U )
75 UL = U(:,2:end-2) ;
76 UR = U(:,3:end-1) ;
77 %%
78 FL = feval( FLUX, UL) ;
79 FR = feval( FLUX, UR) ;
80 % Richtmyer
81 US = 0.5 .* ( UL + UR + (DT/DX) .* (FL-FR) ) ;
82 FRM = feval( FLUX, US ) ;
83 % LaxFriedrichs
84 FLF = 0.5 .* ( FL+FR + (DX/DT) .* (UL-UR) ) ;
85 F = 0.5 .* ( FRM + FLF ) ;
86end
87% ____ _ ____ _ _
88% / ___| ___ __| |_ _ _ __ _____ __/ ___|___ _ __ | |_ ___ _ __ ___ __| |
89% | | _ / _ \ / _` | | | | '_ \ / _ \ \ / / | / _ \ '_ \| __/ _ \ '__/ _ \/ _` |
90% | |_| | (_) | (_| | |_| | | | | (_) \ V /| |__| __/ | | | || __/ | | __/ (_| |
91% \____|\___/ \__,_|\__,_|_| |_|\___/ \_/ \____\___|_| |_|\__\___|_| \___|\__,_|
92%
93% compute an intercell flux according to the Godunov centred scheme (non-monotone).
94% Stability: 0 < CFL Coefficient < 0.5 * sqrt(2), about 0.7
95%
96function F = GodunovCentered( FLUX, DT, DX, U )
97 UL = U(:,2:end-2) ;
98 UR = U(:,3:end-1) ;
99 %%
100 FL = feval( FLUX, UL ) ;
101 FR = feval( FLUX, UR ) ;
102 US = 0.5 .* (UL+UR) + (DT/DX) .* ( FL - FR ) ;
103 % Compute the flux at the state US
104 F = feval( FLUX, US ) ;
105end
106% ____ _ _ _
107% | _ \(_) ___| |__ | |_ _ __ ___ _ _ ___ _ __
108% | |_) | |/ __| '_ \| __| '_ ` _ \| | | |/ _ \ '__|
109% | _ <| | (__| | | | |_| | | | | | |_| | __/ |
110% |_| \_\_|\___|_| |_|\__|_| |_| |_|\__, |\___|_|
111% |___/
112% Compute intercell fluxes according to the Richtmyer method, or two step Lax-Wendroff method
113% See Eq. 5.79, Chapter 5, Ref. 1
114%
115function F = Richtmyer( FLUX, DT, DX, U )
116 UL = U(:,2:end-2) ;
117 UR = U(:,3:end-1) ;
118 %%
119 FL = feval( FLUX, UL) ;
120 FR = feval( FLUX, UR) ;
121 US = 0.5 * ( UL + UR + (DT/DX) * (FL-FR) ) ;
122 F = feval( FLUX, US ) ;
123end