Mercurial > repos > public > sbplib

comparison +noname/testCfl.m @ 116:ee9c03dc7f42

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
noname.testCfl: Improved relyability and clarity.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 14 Dec 2015 19:09:03 +0100
parents 9933169d2651
children 5046ff7d13b8
comparison
equal deleted inserted replaced
115:1fe783681f9f 116:ee9c03dc7f42
1 % noname.testCfl(discr, timestepper_method, T, alpha0, tol,threshold, silentFlag)
2 % Example:
3 % noname.testCfl(Discr(100,4), 'rk4', 1, [0, 1])
1 function testCfl(discr, timestepper_method, T, alpha0, tol,threshold, silentFlag) 4 function testCfl(discr, timestepper_method, T, alpha0, tol,threshold, silentFlag)
2 default_arg('tol',0.00005); 5 default_arg('tol',0.00005);
3 default_arg('threshold',1e2); 6 default_arg('threshold',1e2);
4 default_arg('silentFlag', false); 7 default_arg('silentFlag', false);
5 8
6 alpha0(2) = alpha0(1)+2*(alpha0(2)-alpha0(1)); 9 if T < alpha0(2)
10 error('Upper bound on alpha must be smaller than T');
11 end
7 12
13 testAlpha = getAlphaTester(discr, T, threshold, silentFlag, timestepper_method);
14
15 % Make sure that the upper bound is not working
16 ok = testAlpha(alpha0(2));
17 if ok % Upper bound too large!
18 error('The upper bound on alpha is stable!')
19 end
20
21 % Make sure that the lower bound is ok
22 if alpha0(1) ~= 0
23 ok = testAlpha(alpha0(1));
24 if ~ok
25 error('The lower bound on alpha is unstable!');
26 end
27 end
28
29 % Use bisection to find sharp estimate
8 while( (alpha0(2)-alpha0(1))/alpha0(1) > tol) 30 while( (alpha0(2)-alpha0(1))/alpha0(1) > tol)
9 alpha = mean(alpha0); 31 alpha = mean(alpha0);
32 fprintf('[%.3e,%.3e]: ', alpha0(1), alpha0(2));
33 ok = testAlpha(alpha);
34 if ok
35 alpha0(1) = alpha;
36 else
37 alpha0(2) = alpha;
38 end
39 end
10 40
41 fprintf('T = %-3d dof = %-4d order = %d: clf = %.4e\n',T, discr.size(), discr.order, alpha0(1));
42
43 end
44
45 function f = getAlphaTester(discr, T, threshold, silentFlag, timestepper_method)
46
47 % Returns true if cfl was ok
48 function ok = testAlpha(alpha)
11 ts = discr.getTimestepper(struct('method', timestepper_method, 'cfl', alpha)); 49 ts = discr.getTimestepper(struct('method', timestepper_method, 'cfl', alpha));
12 50
13 warning('off','all') 51 warning('off','all')
14 ts.evolve(T,true); 52 ts.evolve(T,true);
15 warning('on','all') 53 warning('on','all')
16 54
17 [v,t] = ts.getV(); 55 [v,t] = ts.getV();
18 max_val = max(v); 56 maxVal = max(v);
19 57
20 if isnan(max_val) || max_val == Inf || max_val > threshold 58 if isnan(maxVal) || maxVal == Inf || maxVal > threshold
21 alpha0(2) = alpha; 59 ok = false;
22 else 60 else
23 alpha0(1) = alpha; 61 ok = true;
24 end 62 end
25 63
26 if ~silentFlag 64 if ~silentFlag
27 fprintf('[%.3e,%.3e]: a = %.3e, max= %.2e\n',alpha0(1),alpha0(2), alpha,max_val); 65 fprintf('a = %.3e, max= %.2e\n', alpha, maxVal);
28 end 66 end
29 end 67 end
30 68
31 fprintf('T = %-3d dof = %-4d order = %d: clf = %.4e\n',T, discr.size(), discr.order, alpha0(1)); 69 f = @testAlpha;
32
33 end 70 end