Mercurial > repos > public > sbplib

comparison +noname/calculateSolution.m @ 17:30ae48efc7ae

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