Mercurial > repos > public > sbplib

comparison +noname/calculateSolution.m @ 167:15baeb35f74e feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in changes from default.
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 23 Feb 2016 13:25:43 +0100
parents 2b133d833668
children bd99ea1fc733
comparison
equal deleted inserted replaced
166:7cb97c1988d9 167:15baeb35f74e
1 % Calculates the solution of discretization for a given set of ms ts and orders. 1 % Calculates the solution of discretization for a given set of ms ts and orders.
2 % discrHand -- function handle to a Discretization constructor 2 % discrHand -- function handle to a Discretization constructor
3 % method -- time stepping method
4 % m -- grid parameter 3 % m -- grid parameter
5 % order -- order of accuracy of the approximtion 4 % order -- order of accuracy of the approximtion
6 % T -- time to calculate solution for 5 % T -- time to calculate solution for
6 % tsOpt -- options for the time stepper creation.
7 % input paramters m, t, order may all be vectors. 7 % input paramters m, t, order may all be vectors.
8 function [] = calculateSolution(filename, discrHand, method, m, T_in, order, force_flag) 8 function [] = calculateSolution(filename, name, discrHand, m, T_in, order, tsOpt, force_flag)
9 default_arg('force_flag',false); 9 default_arg('force_flag',false);
10 default_arg('tsOpt', []);
10 11
11 if exist(filename,'file') && ~force_flag 12 if exist(filename,'file') && ~force_flag
12 fprintf('File ''%s'' already exist.',filename); 13 fprintf('File ''%s'' already exist.',filename);
13 do_append = yesnoQuestion('Do you want to append to it?'); 14 do_append = yesnoQuestion('Do you want to append to it?');
14 if ~do_append 15 if ~do_append
17 end 18 end
18 end 19 end
19 20
20 sf = SolutionFile(filename); 21 sf = SolutionFile(filename);
21 22
22
23
24 orderWidth = findFieldWidth('%d',order); 23 orderWidth = findFieldWidth('%d',order);
25 mWidth = findFieldWidth('%d',m); 24 mWidth = findFieldWidth('%d',m);
26 TWidth = findFieldWidth('%d',T_in); 25 TWidth = findFieldWidth('%d',T_in);
27 26
28 for i = 1:length(order) 27 for i = 1:length(order)
29 for j = 1:length(m) 28 for j = 1:length(m)
30 T = sort(T_in); % Make sure times are sorted 29 T = sort(T_in); % Make sure times are sorted
31 30
32 discr = discrHand(m(j),order(i)); 31 discr = discrHand(m(j),order(i));
33 k_max = discr.getTimestep(method); 32 k_max = discr.getTimestep(tsOpt);
34 33
35 % Do we want to to save the initial conditions? 34 % Do we want to to save the initial conditions?
36 if T(1) == 0 35 if T(1) == 0
37 snapshot = discr.getTimeSnapshot(0); 36 snapshot = discr.getTimeSnapshot(0);
38 saveToFile(sf, method, order(i), m(j),T(1), snapshot, NaN, NaN, discr); 37 saveToFile(sf, name, order(i), m(j),T(1), snapshot, NaN, NaN, discr);
39 T(1) = []; 38 T(1) = [];
40 end 39 end
41 40
42 % Find out if times to be calulated are integer multiples of the smallest one. 41 % Find out if times to be calulated are integer multiples of the smallest one.
43 time_multiples = T/T(1); 42 time_multiples = T/T(1);
45 is_int_multiples = all(time_multiples == int64(time_multiples)); 44 is_int_multiples = all(time_multiples == int64(time_multiples));
46 45
47 if is_int_multiples 46 if is_int_multiples
48 fprintf('Calculating time series in increments\n'); 47 fprintf('Calculating time series in increments\n');
49 else 48 else
50 fprintf('Restarting for each time in timeseries\n'); 49 fprintf('RESTARTING for each time in timeseries\n');
50 fprintf('If this is not what you want try giving T in integer multiples.\n');
51 end 51 end
52 52
53 % T now contains all the times we need to step to, 53 % T now contains all the times we need to step to,
54 % if T contained 0 it has now been removed. 54 % if T contained 0 it has now been removed.
55 55
56 if is_int_multiples 56 if is_int_multiples
57 % Times are integer multiples, we can save time 57 % Times are integer multiples, we can save time
58 [k,N] = alignedTimestep(k_max,T(1)); 58 [k,N] = alignedTimestep(k_max,T(1));
59 ts = discr.getTimestepper(method,k); 59 tsOpt.k = k;
60 ts = discr.getTimestepper(tsOpt);
60 runtime = 0; 61 runtime = 0;
61 for l = 1:length(T) 62 for l = 1:length(T)
62 end_step = N * time_multiples(l); 63 end_step = N * time_multiples(l);
63 fprintf('[order = %-*d, m = %-*d, T = %-*d]: ',orderWidth,order(i),mWidth,m(j),TWidth,T(l)); 64 fprintf('[order = %-*d, m = %-*d, T = %-*d]: ',orderWidth,order(i),mWidth,m(j),TWidth,T(l));
64 clock_start = tic(); 65 clock_start = tic();
65 ts.stepN(end_step-ts.n,true); 66 ts.stepN(end_step-ts.n,true);
66 runtime = runtime + toc(clock_start); 67 runtime = runtime + toc(clock_start);
67 snapshot = discr.getTimeSnapshot(ts); 68 snapshot = discr.getTimeSnapshot(ts);
68 saveToFile(sf, method, order(i), m(j),T(l), snapshot, runtime, k, discr); 69 saveToFile(sf, name, order(i), m(j),T(l), snapshot, runtime, k, discr);
69 fprintf('Done! (%.3fs)\n',runtime); 70 fprintf('Done! (%.3fs)\n',runtime);
70 end 71 end
71 else 72 else
72 % Times are not interger multiples, we have to start from 0 every time. 73 % Times are not interger multiples, we have to start from 0 every time.
73 for l = 1:length(T) 74 for l = 1:length(T)
74 [k,N] = alignedTimestep(k_max,T(l)); 75 [k,N] = alignedTimestep(k_max,T(l));
75 ts = discr.getTimestepper(method,k); 76 tsOpt.k = k;
77 ts = discr.getTimestepper(tsOpt);
76 fprintf('[order = %-*d, m = %-*d, T = %-*d]: ',orderWidth,order(i),mWidth,m(j),TWidth,T(l)); 78 fprintf('[order = %-*d, m = %-*d, T = %-*d]: ',orderWidth,order(i),mWidth,m(j),TWidth,T(l));
77 clock_start = tic(); 79 clock_start = tic();
78 [v,t] = ts.stepN(N-ts.n,true); 80 [v,t] = ts.stepN(N-ts.n,true);
79 runtime = toc(clock_start); 81 runtime = toc(clock_start);
80 snapshot = discr.getTimeSnapshot(ts); 82 snapshot = discr.getTimeSnapshot(ts);
81 saveToFile(sf, method, order(i), m(j),T(l), snapshot, runtime, k, discr); 83 saveToFile(sf, name, order(i), m(j),T(l), snapshot, runtime, k, discr);
82 fprintf('Done! (%.3fs)\n',runtime); 84 fprintf('Done! (%.3fs)\n',runtime);
83 end 85 end
84 86
85 end 87 end
86 sf.stupidSave(); 88 sf.stupidSave();
87 end 89 end
88 end 90 end
89 end 91 end
90 92
91 93
92 function saveToFile(sf, method, order, m, T, snapshot, runtime, k, discr) 94 function saveToFile(sf, name, order, m, T, snapshot, runtime, k, discr)
93 key.method = method; 95 key.name = name;
94 key.order = order; 96 key.order = order;
95 key.m = m; 97 key.m = m;
96 key.T = T; 98 key.T = T;
97 99
98 entry.repr = snapshot; 100 entry.repr = snapshot;
99 entry.runtime = runtime; 101 entry.runtime = runtime;
100 entry.k = k; 102 entry.k = k;
101 entry.discr = discr; 103 entry.discr = discr;