Mercurial > repos > public > sbplib

comparison +noname/testCfl.m @ 118:5046ff7d13b8

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed bug and made more robust.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 15 Dec 2015 16:38:02 +0100
parents ee9c03dc7f42
children de8bcef865b0
comparison
equal deleted inserted replaced
117:8514c3d67201 118:5046ff7d13b8
4 function testCfl(discr, timestepper_method, T, alpha0, tol,threshold, silentFlag) 4 function testCfl(discr, timestepper_method, T, alpha0, tol,threshold, silentFlag)
5 default_arg('tol',0.00005); 5 default_arg('tol',0.00005);
6 default_arg('threshold',1e2); 6 default_arg('threshold',1e2);
7 default_arg('silentFlag', false); 7 default_arg('silentFlag', false);
8 8
9 if T < alpha0(2) 9 % TODO:
10 error('Upper bound on alpha must be smaller than T'); 10 % Set threshold from the initial conditions of the pde?
11 end 11 % Take a set number of steps instead of evolving to a certain time?
12 % Stop evolving when it has blown up?
12 13
13 testAlpha = getAlphaTester(discr, T, threshold, silentFlag, timestepper_method); 14 testAlpha = getAlphaTester(discr, T, threshold, silentFlag, timestepper_method);
14 15
15 % Make sure that the upper bound is not working 16 % Make sure that the upper bound is not working
16 ok = testAlpha(alpha0(2)); 17 ok = testAlpha(alpha0(2));
27 end 28 end
28 29
29 % Use bisection to find sharp estimate 30 % Use bisection to find sharp estimate
30 while( (alpha0(2)-alpha0(1))/alpha0(1) > tol) 31 while( (alpha0(2)-alpha0(1))/alpha0(1) > tol)
31 alpha = mean(alpha0); 32 alpha = mean(alpha0);
32 fprintf('[%.3e,%.3e]: ', alpha0(1), alpha0(2)); 33
33 ok = testAlpha(alpha); 34 if ~silentFlag
35 fprintf('[%.3e,%.3e]: ', alpha0(1), alpha0(2));
36 end
37
38 [ok, n_step, maxVal] = testAlpha(alpha);
39
34 if ok 40 if ok
35 alpha0(1) = alpha; 41 alpha0(1) = alpha;
36 else 42 else
37 alpha0(2) = alpha; 43 alpha0(2) = alpha;
38 end 44 end
45
46 if ~silentFlag
47 fprintf('a = %.3e, n_step=%d max = %.2e\n', alpha, n_step, maxVal);
48 end
49 end
50
51 if silentFlag
52 fprintf('Last calculated: a = %.3e, n_step=%d max = %.2e\n', alpha, n_step, maxVal);
39 end 53 end
40 54
41 fprintf('T = %-3d dof = %-4d order = %d: clf = %.4e\n',T, discr.size(), discr.order, alpha0(1)); 55 fprintf('T = %-3d dof = %-4d order = %d: clf = %.4e\n',T, discr.size(), discr.order, alpha0(1));
42 56
43 end 57 end
44 58
45 function f = getAlphaTester(discr, T, threshold, silentFlag, timestepper_method) 59 function f = getAlphaTester(discr, T, threshold, silentFlag, timestepper_method)
46 60
47 % Returns true if cfl was ok 61 % Returns true if cfl was ok
48 function ok = testAlpha(alpha) 62 function [ok, n_step, maxVal] = testAlpha(alpha)
49 ts = discr.getTimestepper(struct('method', timestepper_method, 'cfl', alpha)); 63 ts = discr.getTimestepper(struct('method', timestepper_method, 'cfl', alpha));
50 64
51 warning('off','all') 65 warning('off','all')
52 ts.evolve(T,true); 66 ts.evolve(T,true);
53 warning('on','all') 67 warning('on','all')
54 68
55 [v,t] = ts.getV(); 69 [v,t] = ts.getV();
56 maxVal = max(v); 70 maxVal = max(v);
57 71
58 if isnan(maxVal) || maxVal == Inf || maxVal > threshold 72 if isnan(maxVal) || maxVal == Inf || abs(maxVal) > threshold
59 ok = false; 73 ok = false;
60 else 74 else
61 ok = true; 75 ok = true;
62 end 76 end
63 77
64 if ~silentFlag 78 n_step = ts.n;
65 fprintf('a = %.3e, max= %.2e\n', alpha, maxVal);
66 end
67 end 79 end
68 80
69 f = @testAlpha; 81 f = @testAlpha;
70 end 82 end