MATLAB lessons N.3¶
A library to solve 1D hyperbolic equations with various methods
All MATLAB file for the lesson in a unique archive lesson3.zip
File ClassRoomExample.m
:open:
1close all ;
2
3LinearAdvectionModel( 'setup', 1 ) ;
4LinearAdvectionModel( 'set_test', 'classroomproblem' ) ;
5LinearAdvectionModel( 'set_numerical_flux', 'godunov' ) ;
6LinearAdvectionModel( 'set_limiter', 'godunov' ) ;
7
8PRB . CFL = 2*LinearAdvectionModel( 'get_cfl' ) ;
9PRB . FinalTime = 1.5 ;
10PRB . Ncells = 800 ; % select 100 cell inteval subdivision
11PRB . StepTime = 8 ;
12
13% call the solver
14[ UFRAMES, TFRAMES ] = Solver(PRB,@LinearAdvectionModel) ;
15
16U = UFRAMES(:,:,end) ; % extract the numerical solutions at the last computed time
17X = LinearAdvectionModel( 'get_x' ) ;
18
19handleFig = figure('Position',[10 10 1000 600],'color','w') ;
20plot(X, U,'-ob','MarkerFaceColor','b');
Linear Adbvection Model:
:open:
1% solve the linearised equations of Gas Dynamics in one space dimension
2%
3% u(t,x) + c u (t,x) = 0
4% t x
5%
6% for t > 0 and -1 <= x <= 1
7%
8function [varargout] = LinearAdvectionModel( what, varargin )
9
10 % LinearAdvectionModel_paramaters
11 % c = velocity
12 % test = which test is used
13 % nflux = which numerical flux is used
14 % bc_left = left bounday condition
15 % bc_right = right boundaruy condition
16 %
17
18 global LinearAdvectionModel_paramaters ;
19
20 switch lower(what)
21
22 case 'setup'
23 expected( what, 1, nargin, 0, nargout ) ;
24 LinearAdvectionModel_paramaters = clear_struct(LinearAdvectionModel_paramaters) ;
25 LinearAdvectionModel_paramaters . speed = varargin{1} ;
26
27 case 'set_test'
28 expected( what, 1, nargin, 0, nargout ) ; % TEST
29 ChooseTest( varargin{1} ) ;
30
31 case 'set_limiter'
32 expected( what, 1, nargin, 0, nargout ) ;
33 LinearAdvectionModel_paramaters . limiter = lower( varargin{1} ) ;
34
35 case 'set_numerical_flux'
36 expected( what, 1, nargin, 0, nargout ) ;
37 LinearAdvectionModel_paramaters . nflux = lower( varargin{1} ) ;
38
39 case 'get_dim'
40 expected( what, 0, nargin, 1, nargout ) ;
41 varargout{1} = 1 ;
42
43 case 'get_dx'
44 expected( what, 0, nargin, 1, nargout ) ;
45 varargout{1} = LinearAdvectionModel_paramaters . DX ;
46
47 case 'get_x'
48 expected( what, 0, nargin, 1, nargout ) ;
49 varargout{1} = LinearAdvectionModel_paramaters . X ;
50
51 case 'get_final_time'
52 expected( what, 0, nargin, 1, nargout ) ;
53 varargout{1} = LinearAdvectionModel_paramaters . tf ;
54
55 case 'get_cfl'
56 expected( what, 0, nargin, 1, nargout ) ;
57 varargout{1} = MaxCFL() ;
58
59 case 'names'
60 expected( what, 0, nargin, 1, nargout ) ;
61 varargout{1} = 'u' ;
62
63 case 'flux'
64 expected( what, 1, nargin, 1, nargout ) ; % U
65 varargout{1} = Flux( varargin{1} ) ;
66
67 case 'init'
68 expected( what, 1, nargin, 1, nargout ) ; % N
69 varargout{1} = Init( varargin{1} ) ;
70
71 case 'exact'
72 expected( what, 2, nargin, 2, nargout ) ; % TIME, N
73 [ varargout{1} varargout{2} ] = Exact( varargin{1}, varargin{2} ) ;
74
75 case 'muscle_nflux'
76 expected( what, 2, nargin, 1, nargout ) ; % UL, UR
77 varargout{1} = Godunov( varargin{1}, varargin{2} ) ;
78
79 case 'conservative_to_physical'
80 varargout{1} = varargin{1} ;
81
82 case 'physical_to_conservative'
83 varargout{1} = varargin{1} ;
84
85 case 'conservative_to_characteristic'
86 varargout{1} = varargin{1} ;
87
88 case 'characteristic_to_conservative'
89 varargout{1} = varargin{1} ;
90
91 case 'dt_from_cfl'
92 expected( what, 2, nargin, 1, nargout ) ; % U, CFLcoeff
93 varargout{1} = DT_from_CFL( varargin{1}, varargin{2} ) ;
94
95 case 'extrapolate_for_bc'
96 expected( what, 1, nargin, 4, nargout ) ; % U
97 [ varargout{1}, varargout{2}, varargout{3}, varargout{4} ] = Extrapolate_for_BC( varargin{1} ) ;
98
99 case 'nflux'
100 expected( what, 2, nargin, 1, nargout ) ; % DT, U
101 varargout{1} = NFlux( varargin{1}, varargin{2} ) ;
102
103 otherwise
104 error(['LinearAdvectionModel bad argument: ', what]) ;
105 end
106end
107
108function SE = clear_struct(S)
109 if isstruct(S)
110 SE = rmfield( S, fieldnames(S) ) ;
111 else
112 SE = S ;
113 end
114end
115
116function expected( what, ein, nin, eout, nout )
117 if ein ~= nin-1
118 error( ['in LinearAdvectionModel call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
119 end
120 if eout ~= nout
121 error( ['in LinearAdvectionModel call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;
122 end
123end
124
125% #######
126% # # # # # #
127% # # # # # #
128% ##### # # # ##
129% # # # # ##
130% # # # # # #
131% # ###### #### # #
132
133% compute the flux
134function F = Flux( U )
135 global LinearAdvectionModel_paramaters ;
136 speed = LinearAdvectionModel_paramaters . speed ;
137 F = speed .* U ;
138end
139
140% ###
141% # # # # #####
142% # ## # # #
143% # # # # # #
144% # # # # # #
145% # # ## # #
146% ### # # # #
147
148function U = Init(N)
149 global LinearAdvectionModel_paramaters ;
150 % Calculate mesh size DX
151 DX = 2.0/N ;
152 X = [-1+DX/2:DX:1-DX/2] ;
153
154 if isfield(LinearAdvectionModel_paramaters, 'test')
155 switch LinearAdvectionModel_paramaters . test
156 case 'classroomproblem'
157 DX = 8.0/N ;
158 X = [0+DX/2:DX:8-DX/2] ; % center of the cell
159 N1 = floor(3*N/8) ;
160 N2 = N - N1 ;
161 U = [ones(1,N1) zeros(1,N2) ] ;
162 case 'square'
163 N3 = floor(N/3) ;
164 M = N - N3 - N3 ;
165 U = [zeros(1,N3) ones(1,M) zeros(1,N3) ] ;
166 case 'smooth'
167 U = exp( -8 .* X .^ 2 ) ;
168 otherwise
169 error(['LinearAdvectionModel::Init unknown test ', LinearAdvectionModel_paramaters . test]) ;
170 end
171 else
172 error('LinearAdvectionModel::Init unknown test you must set it by set_test!') ;
173 end
174 LinearAdvectionModel_paramaters . DX = DX ;
175 LinearAdvectionModel_paramaters . X = X ;
176end
177
178
179% ##### #######
180% # # # # #### #### #### ###### # ###### #### #####
181% # # # # # # # # # # # # #
182% # ###### # # # # #### ##### # ##### #### #
183% # # # # # # # # # # # # #
184% # # # # # # # # # # # # # # # #
185% ##### # # #### #### #### ###### # ###### #### #
186%
187function ChooseTest(test)
188
189 global LinearAdvectionModel_paramaters ;
190
191 switch test
192 case 'classroomproblem'
193 tf = 3.5 ;
194 bc_left = 'free' ;
195 bc_right = 'free' ;
196 case 'square'
197 tf = 4 ;
198 bc_left = 'cyclic' ;
199 bc_right = 'cyclic' ;
200 case 'smooth'
201 tf = 4 ;
202 bc_left = 'cyclic' ;
203 bc_right = 'cyclic' ;
204 otherwise
205 error(['LinearAdvectionModel:ChooseTest test name',test] ) ;
206 end
207
208 LinearAdvectionModel_paramaters . test = test ;
209 LinearAdvectionModel_paramaters . bc_left = bc_left ;
210 LinearAdvectionModel_paramaters . bc_right = bc_right ;
211 LinearAdvectionModel_paramaters . tf = tf ;
212
213end
214
215
216
217% #######
218% # # # ## #### #####
219% # # # # # # # #
220% ##### ## # # # #
221% # ## ###### # #
222% # # # # # # # #
223% ####### # # # # #### #
224
225function [X,U] = Exact(TIME,N)
226 global LinearAdvectionModel_paramaters ;
227 speed = LinearAdvectionModel_paramaters . speed ;
228
229 % Calculate mesh size DX
230 DX = 2.0/N ;
231 X = [-1+DX/2:DX:1-DX/2] ;
232 X = X - speed*TIME ; % translate the mesh
233
234 % put the solution in [-1..1]
235 while min(X) < -1
236 IDX = find( X < -1 ) ;
237 X(IDX) = X(IDX)+2 ;
238 end
239
240 while max(X) > 1
241 IDX = find( X > 1 ) ;
242 X(IDX) = X(IDX)-2 ;
243 end
244
245 switch LinearAdvectionModel_paramaters . test
246 case 'square'
247 U = zeros(1,N) ;
248 IDX = find( X >= -1.0/3.0 & X <= 1.0/3.0 ) ;
249 U(IDX) = 1 ;
250 case 'smooth'
251 % use gauss quadrature to estimate status
252 XL = X + DX/sqrt(3)/2;
253 XR = X - DX/sqrt(3)/2 ;
254 U = (exp( -8 .* XL .^ 2 )+exp( -8 .* XR .^ 2 ))/2 ;
255 end
256 % ripristina la mesh
257 X = [-1+DX/2:DX:1-DX/2] ;
258end
259
260% ###### ####### ##### ####### #
261% # # # ###### ##### #### # # # # # #
262% # # # # # # # # ## ## # # #
263% # # # ##### # # # # # ## # # ##### #
264% # # # # ##### # # # # # # #
265% # # # # # # # # # # # # # #
266% ###### # # # # #### # # ##### # #######
267% ####### #######
268% Compute DT from CFL
269function DT = DT_from_CFL( U, CFLcoeff )
270 global LinearAdvectionModel_paramaters ;
271 speed = LinearAdvectionModel_paramaters . speed ;
272 DX = LinearAdvectionModel_paramaters . DX ;
273 DT = CFLcoeff * DX / abs(speed) ;
274end
275
276% #######
277% # # # ##### ##### ## ##### #### # ## ##### ######
278% # # # # # # # # # # # # # # # # #
279% ##### ## # # # # # # # # # # # # # #####
280% # ## # ##### ###### ##### # # # ###### # #
281% # # # # # # # # # # # # # # # #
282% ####### # # # # # # # # #### ###### # # # ######
283%
284% ###### #####
285% ###### #### ##### # # # #
286% # # # # # # # #
287% ##### # # # # ###### #
288% # # # ##### # # #
289% # # # # # # # # #
290% # #### # # ###### #####
291% #######
292
293% Extrapolate values for boundary condition
294function [ULL, UL, UR, URR] = Extrapolate_for_BC( U )
295 global LinearAdvectionModel_paramaters ;
296 L = LinearAdvectionModel_paramaters . bc_left ;
297 R = LinearAdvectionModel_paramaters . bc_right ;
298 % Left boundary
299 switch L
300 case 'free' % Apply transmissive boundary conditions
301 UL = U(:,1) ;
302 ULL = U(:,2) ;
303 case 'cyclic' % Apply cyclic boundary conditions
304 UL = U(:,end) ;
305 ULL = U(:,end-1) ;
306 otherwise % Apply reflective boundary conditions
307 UL = [ U(1,1) ; -U(2,1) ] ;
308 ULL = [ U(1,2) ; -U(2,2) ] ;
309 end
310 % Right boundary
311 switch R
312 case 'free' % Apply transmissive boundary conditions
313 UR = U(:,end) ;
314 URR = U(:,end-1) ;
315 case 'cyclic' % Apply cyclic boundary conditions
316 UR = U(:,1) ;
317 URR = U(:,2) ;
318 otherwise % Apply reflective boundary conditions
319 UR = [ U(1,end) ; -U(2,end) ] ;
320 URR = [ U(1,end-1) ; -U(2,end-1) ] ;
321 end
322end
323
324% ####### ##### ####### #
325% # # ##### ##### # # # ## # # # # #
326% # # # # # # ## ## # # # # # #
327% # # # # # # # ## # # # # # ##### #
328% # # ##### # # # # ###### # # # #
329% # # # # # # # # # # # # # #
330% ####### # # # # # # # ###### ##### # #######
331
332function CFL = MaxCFL
333 global LinearAdvectionModel_paramaters ;
334 nflux = LinearAdvectionModel_paramaters . nflux ;
335 %
336 switch nflux
337 case 'torobillett'
338 CFL = 2.0 ;
339 case { 'fromm', 'warmingbeam', 'waf', 'godunov' }
340 CFL = 1.0 ;
341 otherwise
342 CFL = StandardFlux( 'get_cfl', nflux ) ;
343 end
344end
345
346
347% # # # # # # ###### ##### # #### ## #
348% ## # # # ## ## # # # # # # # # #
349% # # # # # # ## # ##### # # # # # # #
350% # # # # # # # # ##### # # ###### #
351% # ## # # # # # # # # # # # # #
352% # # #### # # ###### # # # #### # # ######
353%
354% ###### # # # # # ###### ####
355% # # # # # # # #
356% ##### # # # ## ##### ####
357% # # # # ## # #
358% # # # # # # # # #
359% # ###### #### # # ###### ####
360
361function FLUX = NFlux( DT, U )
362
363 global LinearAdvectionModel_paramaters ;
364
365 ULL = U(:,1:end-3) ;
366 UL = U(:,2:end-2) ;
367 UR = U(:,3:end-1) ;
368 URR = U(:,4:end-0) ;
369
370 nflux = LinearAdvectionModel_paramaters . nflux ;
371 DX = LinearAdvectionModel_paramaters . DX ;
372 limiter = LinearAdvectionModel_paramaters . limiter ;
373
374 switch nflux
375 % monolitic numerical fluxes
376 case 'godunov'
377 FLUX = Godunov(UL,UR) ;
378 case 'torobillett'
379 FLUX = ToroBillett(DT,ULL,UL,UR,URR) ;
380 case 'fromm'
381 FLUX = Fromm(DT,ULL,UL,UR,URR) ;
382 case 'warmingbeam'
383 FLUX = WarmingBeam(DT,ULL,UL,UR,URR) ;
384 case 'waf'
385 FLUX = Waf(DT,U) ;
386 otherwise
387 FLUX = StandardFlux( 'nflux', nflux, @Flux, DT, DX, U, limiter ) ;
388 end
389end
390
391% _ ___ _ __ ___ _ _ _ _ _
392% |_ | |_) (_ | / \ |_) | \ |_ |_)
393% | _|_ | \ __) | \_/ | \ |_/ |_ | \
394%
395
396% #####
397% # # #### ##### # # # # #### # #
398% # # # # # # # ## # # # # #
399% # #### # # # # # # # # # # # # #
400% # # # # # # # # # # # # # # #
401% # # # # # # # # # ## # # # #
402% ##### #### ##### #### # # #### ## %
403%
404% Compute the solution of the Riemann problem
405%
406function F = Godunov( UL, UR )
407 global LinearAdvectionModel_paramaters ;
408 speed = LinearAdvectionModel_paramaters . speed ;
409 % the flux at interface i is computed
410 if speed > 0
411 F = speed .* UL ; % using the left state
412 else
413 F = speed .* UR ; % using the right state
414 end
415end
416% ####### ######
417% # #### ##### #### # # # # # ###### ##### #####
418% # # # # # # # # # # # # # # #
419% # # # # # # # ###### # # # ##### # #
420% # # # ##### # # # # # # # # # #
421% # # # # # # # # # # # # # # #
422% # #### # # #### ###### # ###### ###### ###### # #
423%
424% Compute intercell fluxes according to the Toro-Billett first order CFL-2 scheme
425% See Eq. 13.22, Chapter 13, Ref. 1
426function F = ToroBillett(DT,ULL,UL,UR,URR)
427 global LinearAdvectionModel_paramaters ;
428 speed = LinearAdvectionModel_paramaters . speed ;
429 DX = LinearAdvectionModel_paramaters . DX ;
430 CN = abs(speed*DT/DX) ;
431
432 if CN <= 1
433 % Conventional Godunov upwind flux
434 F = Godunov(UL,UR) ;
435 else
436 a = (CN-1)/CN ;
437 b = 1-a ;
438 if speed > 0
439 % Information travels from the left
440 F = a .* Flux(ULL) + b .* Flux(UL) ;
441 else
442 % Information travels from the right
443 F = a .* Flux(URR) + b .* Flux(UR) ;
444 end
445 end
446end
447
448% __ _ _ _ _ _ _ _ _ _
449% (_ |_ / / \ |\ | | \ / \ |_) | \ |_ |_)
450% __) |_ \_ \_/ | \| |_/ \_/ | \ |_/ |_ | \
451
452% #######
453% # ##### #### # # # #
454% # # # # # ## ## ## ##
455% ##### # # # # # ## # # ## #
456% # ##### # # # # # #
457% # # # # # # # # #
458% # # # #### # # # #
459%
460% Compute intercell fluxes according to the Fromm first order upwind method
461% See Eq. 13.35, Chapter 13, Ref. 1
462%
463function F = Fromm(DT,ULL,UL,UR,URR)
464 global LinearAdvectionModel_paramaters ;
465 speed = LinearAdvectionModel_paramaters . speed ;
466 DX = LinearAdvectionModel_paramaters . DX ;
467 CN = speed*DT/DX ;
468 if CN > 0
469 F = speed .* (UL + 0.25 .* (1-CN) .* (UR - ULL)) ;
470 else
471 F = speed .* (UR - 0.25 .* (1+CN) .* (URR - UR)) ;
472 end ;
473end
474
475% # # ######
476% # # # ## ##### # # # # # #### # # ###### ## # #
477% # # # # # # # ## ## # ## # # # # # # # # ## ##
478% # # # # # # # # ## # # # # # # ###### ##### # # # ## #
479% # # # ###### ##### # # # # # # # ### # # # ###### # #
480% # # # # # # # # # # # ## # # # # # # # # #
481% ## ## # # # # # # # # # #### ###### ###### # # # #
482%
483% Compute intercell fluxes according to the Warming-Beam method
484% See Eq. 13.21, Chapter 13, Ref. 1
485%
486function FLUX = WarmingBeam(DT,ULL,UL,UR,URR)
487 global LinearAdvectionModel_paramaters ;
488 speed = LinearAdvectionModel_paramaters . speed ;
489 DX = LinearAdvectionModel_paramaters . DX ;
490 CN = speed*DT/DX ;
491 if speed > 0
492 FLUX = (0.5*speed) .* ( (3-CN) .* UL - (1-CN) .* ULL ) ;
493 else
494 FLUX = (0.5*speed) .* ( (3+CN) .* UR - (1+CN) .* URR ) ;
495 end
496end
497
498% # # # #######
499% # # # # # #
500% # # # # # #
501% # # # # # #####
502% # # # ####### #
503% # # # # # #
504% ## ## # # #
505%
506% Compute the WAF flux with the exact Riemann solver,
507% to be used in explicit conservative formula, subroutine UPDATE
508
509function FLUX = Waf( DT, U )
510 global LinearAdvectionModel_paramaters ;
511 speed = LinearAdvectionModel_paramaters . speed ;
512 DX = LinearAdvectionModel_paramaters . DX ;
513 limiter = LinearAdvectionModel_paramaters . limiter ;
514 CN = speed*DT/DX ;
515
516 % Compute local jump DLOC. This is the change across the wave of speed c.
517 DLOC = U(:,3:end-1) - U(:,2:end-2) ;
518
519 % Identify the upwind direction
520 if speed > 0 % Information travels from the left
521 IUPW = -1 ;
522 else % Information travels from the right
523 IUPW = +1 ;
524 end
525
526 % Compute the upwind jump DUPW. This depends on the direction of the wave
527 % in the solution of the Riemann problem RP(UL, UR), namely the sign of c
528 DUPW = U(:,3+IUPW:end-1+IUPW) - U(:,2+IUPW:end-2+IUPW) ;
529 ONES = ones (size(DLOC)) ;
530
531 % Apply TVD condition to find WAF limiter function WAFLIM
532 WAFLIM = WafLimiter( DLOC, DUPW, CN, limiter ) ;
533 WL = 0.5 .* (ONES + WAFLIM) ;
534 WR = 0.5 .* (ONES - WAFLIM) ;
535
536 % Compute fluxes on data
537 FD = Flux( U(:,2:end-1) ) ;
538
539 % Compute WAF Flux FLUX(I) corresponding to RP(I,I+1)
540 FLUX = WL .* FD(1:end-1) + WR .* FD(2:end) ;
541
542end
543
544function WAFLIM = WafLimiter( DLOC, DUPW, CN, limiter )
545
546 TOLLIM = 1E-6 ;
547
548 % Select limiter function WAFLIM
549 ZEROS = zeros(size(DLOC)) ;
550 ONES = ones (size(DLOC)) ;
551
552 % Compute RATIO=DUPW/DLOC of upwind to local changes
553 % some workaround to avoid division by zero are added
554 RATIO = sign(DLOC).*DUPW./max(abs(DLOC),TOLLIM.*ONES) ;
555
556 switch ( limiter )
557 case 'godunov'
558 WAFLIM = sign(CN) .* ONES ;
559 return ;
560 case 'secondorder'
561 WAFLIM = CN .* ONES ;
562 return ;
563 case 'superbee'
564 B = max( min(2.*RATIO, ONES), min(RATIO, 2.*ONES) ) ;
565 case 'vanleer'
566 B = max(RATIO,ZEROS) ;
567 B = 2.*B./(ONES+B) ;
568 case 'vanalbada'
569 B = (RATIO + RATIO.^2 ) ./ (ONES + RATIO.^2) ;
570 case 'minmod'
571 B = min(RATIO, ONES) ;
572 end
573
574 % Transform to WAF limiter
575 WAFLIM = sign(CN) .* ( ONES - (ONES - abs(CN)) .* max(ZEROS, B) ) ;
576end
Slope Limiter Model:
File SlopeLimiterModel.m
:open:
1function [varargout] = SlopeLimiterModel( what, varargin )
2
3 global SlopeLimiter_parameters ;
4
5 switch ( what )
6
7 case 'setup'
8 expected( what, 1, nargin, 0, nargout ) ;
9 SlopeLimiter_parameters . limiter = lower( varargin{1} ) ;
10
11 case 'eval_slope'
12 expected( what, 1, nargin, 1, nargout ) ;
13 varargout{1} = SlopeLimiter( varargin{1} ) ;
14
15 otherwise
16 error(['SlopeLimiterModel bad argument: ' what]) ;
17 end
18end
19
20function expected( what, ein, nin, eout, nout )
21 if ein ~= nin-1
22 error( ['in SlopeLimiterModel call to ', what, ' expected ', int2str( ein ), ' input argument, found ', int2str(nin-1) ]) ;
23 end
24 if eout ~= nout
25 error( ['in SlopeLimiterModel call to ', what, ' expected ', int2str( eout ), ' output argument, found ', int2str(nout) ]) ;
26 end
27end
28
29function DELTA = SlopeLimiter( US )
30 global SlopeLimiter_parameters ;
31 limiter = SlopeLimiter_parameters . limiter ;
32 [r,c] = size(US) ;
33
34 for k=1:r
35 UL = US(r,1:end-2) ;
36 UC = US(r,2:end-1) ;
37 UR = US(r,3:end-0) ;
38
39 P = UR - UC ;
40 Q = UC - UL ;
41 EPSI = 1E-10 * ones(size(Q)) ;
42
43 % Compute RATIO=DUPW/DLOC of upwind to local changes
44 % some workaround to avoid division by zero are added
45 R = sign(Q).*P./max(abs(Q),EPSI) ;
46 S = sign(P).*Q./max(abs(P),EPSI) ;
47 DELTA(r,:) = min( Q .* StandardLimiter( limiter, R ), ...
48 P .* StandardLimiter( limiter, S ) ) ;
49 end
50end
Check Order for Linear Advection:
File LinearAdvectionCheckOrder.m
:open:
1nflux = { 'LaxFriedrichs', 'force', 'Richtmyer', ...
2 'slic', 'GodunovCentered', 'ToroBillett', ...
3 'Fromm', 'WarmingBeam', 'Waf' } ;
4grid = [ 50 100 200 400 ];
5
6for j=1:length(nflux)
7 LinearAdvectionModel( 'setup', -0.4 ) ;
8 LinearAdvectionModel( 'set_test', 'smooth' ) ;
9 LinearAdvectionModel( 'set_numerical_flux', nflux{j} ) ;
10 LinearAdvectionModel( 'set_limiter', 'superbee' ) ;
11
12 PRB . CFL = 0.9*LinearAdvectionModel( 'get_cfl' ) ;
13 PRB . FinalTime = LinearAdvectionModel( 'get_final_time' ) ;
14 PRB . StepTime = PRB . FinalTime / 100 ;
15
16 SlopeLimiterModel( 'setup', 'superbee' ) ;
17 % limiters:
18 % 'minbee', 'superbee', 'vanleer', 'vanalbada', 'minmod', 'warmingbeam', 'laxwendroff', 'fromm'
19
20 % call the solver
21 for i=1:4
22 PRB . Ncells = grid(i) ;
23 [UFRAMES,TFRAMES] = Solver(PRB,@LinearAdvectionModel) ;
24 %[UFRAMES,TFRAMES] = MUSCLSolver(PRB,@LinearAdvectionModel) ;
25 [NR1(i),NR2(i),NRi(i)] = ComputeError(UFRAMES,TFRAMES,@LinearAdvectionModel);
26 end
27 disp( sprintf('\n\n%s', nflux{j}) ) ;
28 for i=1:3
29 disp( sprintf('norm1[%d,%d] = %f',grid(i+1),grid(i),log(NR1(i)/NR1(i+1))/log(2))) ;
30 disp( sprintf('norm2[%d,%d] = %f',grid(i+1),grid(i),log(NR2(i)/NR2(i+1))/log(2))) ;
31 disp( sprintf('normi[%d,%d] = %f',grid(i+1),grid(i),log(NRi(i)/NRi(i+1))/log(2))) ;
32 end
33end
Linear Advection Example:
:open:
1test = { 'square', 'smooth' } ;
2
3nflux = {'ToroBillet'} ;
4
5nflux1 = { 'LaxFriedrichs', ...
6 'force', ...
7 'Richtmyer', ...
8 'slic', ...
9 'GodunovCentered', ...
10 'ToroBillet', ...
11 'Fromm', ...
12 'WarmingBeam', ...
13 'Waf' } ;
14
15close all ;
16
17for j=1:length(nflux)
18
19
20 ok = input( ['Start numerical flux: ' nflux{j} ' [q=quit] ' ], 's' ) ;
21 if ok == 'q'
22 break ;
23 end ;
24
25 handleFig = figure('Position',[10 10 1000 600],'color','w') ;
26
27 for i=1:length(test)
28
29 LinearAdvectionModel( 'setup', -0.4 ) ;
30 LinearAdvectionModel( 'set_test', test{i} ) ;
31 LinearAdvectionModel( 'set_numerical_flux', nflux{j} ) ;
32 LinearAdvectionModel( 'set_limiter', 'superbee' ) ;
33
34 PRB . CFL = 0.9*LinearAdvectionModel( 'get_cfl' ) ;
35 PRB . FinalTime = LinearAdvectionModel( 'get_final_time' ) ;
36 PRB . Ncells = 100 ; % select 100 cell inteval subdivision
37 PRB . StepTime = PRB . FinalTime / 100 ;
38
39 SlopeLimiterModel( 'setup', 'superbee' ) ;
40 % others limiters:
41 % 'minbee', 'superbee', 'vanleer', 'vanalbada', 'minmod', 'warmingbeam', 'laxwendroff', 'fromm'
42
43 % call the solver
44 [ UFRAMES2, TFRAMES ] = Solver(PRB,@LinearAdvectionModel) ;
45
46 PRB . CFL = 0.9 ; % change CFL for MUSCL
47 [ UFRAMES1, TFRAMES ] = MUSCLSolver(PRB, @LinearAdvectionModel, @SlopeLimiterModel) ;
48 [ Xe, Ue ] = LinearAdvectionModel( 'exact', PRB . FinalTime, 1000 ) ;
49
50 U1 = UFRAMES1(:,:,end) ;
51 U2 = UFRAMES2(:,:,end) ;
52 X = LinearAdvectionModel( 'get_x' ) ;
53
54 subplot(length(test),1,i);
55 plot(Xe, Ue,'-b', 'LineWidth',2, 'Color', 'green') ; hold on ;
56 plot(X, U1,'-ok','MarkerFaceColor','r');
57 plot(X, U2,'-ob','MarkerFaceColor','b');
58 title(['test case: ' test{i} ' solver: ' nflux{j} ' and MUSCL'], ...
59 'FontSize', 24, 'FontName', 'Arial' );
60 axis([X(1) X(end) -0.2 1.2]) ;
61 end
62end
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 case 'slic'
36 varargout{1} = Slic( Flux, DT, DX, varargin{6}, U ) ;
37 otherwise
38 error( [ 'StandardFlux nflux: ', varargin{1}, ' unknown' ] ) ;
39 end
40
41 case 'get_cfl'
42 switch numerical_flux_name
43 case { 'laxfriedrichs', 'force', 'richtmyer' }
44 varargout{1} = 1.0 ;
45 case { 'slic', 'godunovcentered' }
46 varargout{1} = 0.5 ;
47 otherwise
48 error( [ 'StandardFlux nflux: ', numerical_flux_name, ' unknown' ] ) ;
49 end
50 end
51end
52% _ _____ _ _ _ _
53% | | __ ___ __ | ___| __(_) ___ __| |_ __(_) ___| |__ ___
54% | | / _` \ \/ /____| |_ | '__| |/ _ \/ _` | '__| |/ __| '_ \/ __|
55% | |__| (_| |> <_____| _|| | | | __/ (_| | | | | (__| | | \__ \
56% |_____\__,_/_/\_\ |_| |_| |_|\___|\__,_|_| |_|\___|_| |_|___/
57%
58% Compute the Lax-Friedrichs flux
59%
60function F = LaxFriedrichs( FLUX, DT, DX, U )
61 UL = U(:,2:end-2) ;
62 UR = U(:,3:end-1) ;
63 %%
64 FL = feval( FLUX, UL) ;
65 FR = feval( FLUX, UR) ;
66 F = 0.5 .* ( FL + FR + (DX/DT) .* (UL-UR) ) ;
67end
68% _____
69% | ___|__ _ __ ___ ___
70% | |_ / _ \| '__/ __/ _ \
71% | _| (_) | | | (_| __/
72% |_| \___/|_| \___\___|
73%
74% compute an intercell flux according to the FORCE scheme.
75%
76function F = Force( FLUX, DT, DX, U )
77 UL = U(:,2:end-2) ;
78 UR = U(:,3:end-1) ;
79 %%
80 FL = feval( FLUX, UL) ;
81 FR = feval( FLUX, UR) ;
82 % Richtmyer
83 US = 0.5 .* ( UL + UR + (DT/DX) .* (FL-FR) ) ;
84 FRM = feval( FLUX, US ) ;
85 % LaxFriedrichs
86 FLF = 0.5 .* ( FL+FR + (DX/DT) .* (UL-UR) ) ;
87 F = 0.5 .* ( FRM + FLF ) ;
88end
89% ____ _ ____ _ _
90% / ___| ___ __| |_ _ _ __ _____ __/ ___|___ _ __ | |_ ___ _ __ ___ __| |
91% | | _ / _ \ / _` | | | | '_ \ / _ \ \ / / | / _ \ '_ \| __/ _ \ '__/ _ \/ _` |
92% | |_| | (_) | (_| | |_| | | | | (_) \ V /| |__| __/ | | | || __/ | | __/ (_| |
93% \____|\___/ \__,_|\__,_|_| |_|\___/ \_/ \____\___|_| |_|\__\___|_| \___|\__,_|
94%
95% compute an intercell flux according to the Godunov centred scheme (non-monotone).
96% Stability: 0 < CFL Coefficient < 0.5 * sqrt(2), about 0.7
97%
98function F = GodunovCentered( FLUX, DT, DX, U )
99 UL = U(:,2:end-2) ;
100 UR = U(:,3:end-1) ;
101 %%
102 FL = feval( FLUX, UL ) ;
103 FR = feval( FLUX, UR ) ;
104 US = 0.5 .* (UL+UR) + (DT/DX) .* ( FL - FR ) ;
105 % Compute the flux at the state US
106 F = feval( FLUX, US ) ;
107end
108% ____ _ _ _
109% | _ \(_) ___| |__ | |_ _ __ ___ _ _ ___ _ __
110% | |_) | |/ __| '_ \| __| '_ ` _ \| | | |/ _ \ '__|
111% | _ <| | (__| | | | |_| | | | | | |_| | __/ |
112% |_| \_\_|\___|_| |_|\__|_| |_| |_|\__, |\___|_|
113% |___/
114% Compute intercell fluxes according to the Richtmyer method, or two step Lax-Wendroff method
115% See Eq. 5.79, Chapter 5, Ref. 1
116%
117function F = Richtmyer( FLUX, DT, DX, U )
118 UL = U(:,2:end-2) ;
119 UR = U(:,3:end-1) ;
120 %%
121 FL = feval( FLUX, UL) ;
122 FR = feval( FLUX, UR) ;
123 US = 0.5 * ( UL + UR + (DT/DX) * (FL-FR) ) ;
124 F = feval( FLUX, US ) ;
125end
126%
127% ____ _ _
128% / ___|| (_) ___
129% \___ \| | |/ __|
130% ___) | | | (__
131% |____/|_|_|\___|
132%
133function F = Slic( FLUX, DT, DX, limiter, U )
134 %%
135 COE1 = 0.5*DT/DX ;
136 COE2 = 0.5*DX/DT ;
137
138 % Compute slope limiter functions.
139 SlopeLimiterModel('setup', limiter ) ;
140 DELTA = SlopeLimiterModel( 'eval_slope', U ) ;
141
142 % Boundary extrapolated values PIL and PIR are computed
143 PIL = U(:,2:end-1) - 0.5.* DELTA ;
144 PIR = U(:,2:end-1) + 0.5.* DELTA ;
145
146 % Compute boundary extrapolated fluxes
147 FIL = feval( FLUX, PIL ) ;
148 FIR = feval( FLUX, PIR ) ;
149
150 % Evolve boundary extrapolated values
151 % for conservedvariables in each cell i
152 DELFLUX = COE1.*(FIL - FIR) ;
153 EL = PIL + DELFLUX ;
154 ER = PIR + DELFLUX ;
155
156 % Compute intercell flux according to the FORCE method
157 CDL = ER(:,1:end-1) ;
158 CDR = EL(:,2:end) ;
159 FDL = feval( FLUX, CDL ) ;
160 FDR = feval( FLUX, CDR ) ;
161
162 % Compute Lax-Friedrichs flux
163 FLF = 0.5 .* (FDL+FDR) + COE2 .* (CDL-CDR) ;
164
165 % Compute Richtmyer state
166 URI = 0.5 .* (CDL+CDR) + COE1 .* (FDL-FDR) ;
167
168 % Compute Richtmyer flux
169 FRI = feval( FLUX, URI ) ;
170
171 % Compute FORCE intercell flux FLUX
172 F = 0.5 .*(FLF + FRI) ;
173end
Standard Limiter:
File StandardLimiter.m
:open:
1% #####
2% # # ##### ## # # ##### ## ##### #####
3% # # # # ## # # # # # # # # #
4% ##### # # # # # # # # # # # # # #
5% # # ###### # # # # # ###### ##### # #
6% # # # # # # ## # # # # # # # #
7% ##### # # # # # ##### # # # # #####
8%
9% #
10% # # # # # ##### ###### #####
11% # # ## ## # # # # #
12% # # # ## # # # ##### # #
13% # # # # # # # #####
14% # # # # # # # # #
15% ####### # # # # # ###### # #
16%
17function [varargout] = StandardLimiter( what, varargin )
18 if nargin > 1
19 R = varargin{1} ;
20 Z = zeros(size(R)) ;
21 ONE = ones(size(R)) ;
22 end
23 switch lower(what)
24 case 'charm' % [not 2nd order TVD]
25 % Zhou, G, (1995), Numerical simulations of physical discontinuities in single and multi-fluid flows
26 % for arbitrary Mach numbers, PhD Thesis, Chalmers Univ. of Tech., Goteborg, Sweden.
27 varargout{1} = R.*(3*R+1)./((R+1).^2) ;
28 case 'hcus' % [not 2nd order TVD]
29 % Waterson, N P and H Deconinck, (1995), A unified approach to the design and application of bounded
30 % higher-order convection schemes, VKI Preprint 1995-21.
31 varargout{1} = 1.5*(R+abs(R))./(R+2) ;
32 case 'hquick' % [not 2nd order TVD]
33 % Waterson, N P and H Deconinck, (1995)
34 varargout{1} = 2*(R+abs(R))./(R+3) ;
35 case 'koren'
36 % Koren, B, (1993), A robust upwind discretisation method for advection, diffusion and source terms,
37 % In: Numerical Methods for Advection-Diffusion Problems, Ed. C.B.Vreugdenhil & B.Koren, Vieweg, Braunschweig, p117.
38 varargout{1} = max(Z,min(2*R,min((1+2*R)/3,2*ONE))) ;
39 case { 'minmod', 'minbee' }
40 % Roe, P L, (1986), Characteristic-based schemes for the Euler equations, Ann. Rev. Fluid Mech., 18, p337.
41 varargout{1} = max(Z,min(ONE,R)) ;
42 case { 'mc', 'monotonizedcentraldifference' } % symmetric
43 % Van Leer, B, (1977), Towards the ultimate conservative difference scheme III.
44 % Upstream-centered finite-difference schemes for ideal compressible flow., J. Comp. Phys., 23, p263-75.
45 varargout{1} = max(Z,min((1+R)/2, 2*min(R,ONE))) ;
46 case 'osher'
47 % Chatkravathy, S R and S Osher, (1983), High resolution applications of the Osher upwind scheme for the Euler equations,
48 % AIAA Paper 83-1943, Proc. AIAA 6th Comutational Fluid Dynamics Conference, pp 363-73.
49 if nargin > 2
50 beta = varargin{2} ;
51 else
52 beta = 1.5 ;
53 end
54 varargout{1} = max(Z,min(R,beta*ONE)) ;
55 case 'ospre' % symmetric
56 % Waterson, N P and H Deconinck, (1995)
57 % A unified approach to the design and application of bounded higher-order convection schemes, VKI Preprint 1995-21.
58 R2 = R.*(R+1) ;
59 varargout{1} = 1.5 * R2./(R2+1);
60 case 'smart' % [not 2nd order TVD]
61 % Gaskell, P H and A K C Lau, (1988)
62 % Curvature-compensated convective transport: SMART, a new boundedness-preserving transport algorithm
63 % Int. J. Num. Meth. Fluids, 8, p617.
64 varargout{1} = max(Z, min(2*R,min(4*ONE, 0.25+0.75*R))) ;
65 case 'superbee' % symmetric
66 % Roe, P L, (1986), Characteristic-based schemes for the Euler equations, Ann. Rev. Fluid Mech., 18, p337.
67 varargout{1} = max(Z, max(min(ONE,2*R), min(2*ONE,R))) ;
68 case 'sweby' % symmetric
69 % Sweby, P K, (1984), High resolution schemes using flux-limiters for hyperbolic conservation laws.
70 % SIAM J. Num. Anal., 21, p995-1011.
71 if nargin > 2
72 beta = varargin{2} ;
73 else
74 beta = 1.5 ;
75 end
76 varargout{1} = max(Z,max(min(beta*R,ONE),min(R,beta*ONE))) ;
77 case 'umist'
78 % Lien, F S, and M. A. Leschziner, (1994)
79 % Upstream monotonic interpolation for scalar transport with application to complex turbulent flows
80 % Int. J. Num. Meth. Fluids, 19, p527.
81 varargout{1} = max(Z,min(2*min(R,ONE),min(0.25+0.75*R,0.75+0.25*R))) ;
82 case 'vanalbada' % symmetric
83 % Van Albada, G D, B. Van Leer and W. W. Roberts, (1982)
84 % A comparative study of computational methods in cosmic gas dynamics, Astron. Astrophysics, 108, p76.
85 varargout{1} = (R.^2+R) ./ (R.^2+1) ;
86 case 'vanalbada1' % Alternative form [not 2nd order TVD] used on high spatial order schemes
87 % Kermani, M. J., Gerber, A. G., and Stockie, J. M. (2003),
88 % Thermodynamically Based Moisture Prediction Using Roe’s Scheme, 4th Conference of Iranian AeroSpace Society
89 % Amir Kabir University of Technology, Tehran, Iran, January 27–29.
90 varargout{1} = 2*R ./ (R.^2+1) ;
91 case 'vanleer' % symmetric
92 % Van Leer, B, (1974), Towards the ultimate conservative difference scheme II.
93 % Monotonicity and conservation combined in a second order scheme. J. Comp. Phys., 14, p361-70.
94 varargout{1} = (R+abs(R)) ./ (R+1) ;
95 case 'upwind'
96 varargout{1} = Z ;
97 case 'laxwendroff'
98 varargout{1} = ONE ;
99 case 'beamwarming'
100 varargout{1} = R ;
101 case 'fromm'
102 varargout{1} = 0.5*(1+R) ;
103 case 'names'
104 varargout{1} = { 'charm', 'hcus', 'hquick', 'koren', 'minmod', 'mc', 'osher', 'ospre', 'smart', ...
105 'superbee', 'sweby', 'umist', 'vanalbada', 'vanalbada1', 'vanleer', 'upwind', ...
106 'laxwendroff', 'beamwarming', 'fromm' } ;
107 otherwise
108 error( [ 'StandardLimiter: ', what, ' unknown' ] ) ;
109 end
110end
Compute Error:
File ComputeError.m
:open:
1%
2% Compute error in various norm:
3%
4function [NR1,NR2,NRi] = ComputeError( UFRAMES, TFRAMES, MODEL )
5 [dummy,nframes] = size(TFRAMES) ;
6 Un = UFRAMES(:,:,nframes) ;
7 [M,N] = size(Un) ;
8 [Xe,Ue] = feval( MODEL, 'exact', TFRAMES(nframes), N ) ;
9 Udiff = Ue - Un ;
10 NR1 = norm(Udiff,1)/N ;
11 NR2 = norm(Udiff,2)/sqrt(N) ;
12 NRi = norm(Udiff,inf) ;
13end