Mercurial > repos > public > sbplib

comparison +time/CdiffNonlin.m @ 4:8e14b5a577a6

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Attempt att optimization of CdiffNonlin.
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 18 Sep 2015 15:12:44 +0200
parents 5ae4f23d9130
children b18d3d201a71
comparison
equal deleted inserted replaced
3:5205251db8c3 4:8e14b5a577a6
11 end 11 end
12 12
13 13
14 methods 14 methods
15 function obj = CdiffNonlin(D, E, S, k, t0, v, v_prev) 15 function obj = CdiffNonlin(D, E, S, k, t0, v, v_prev)
16 default_arg('S',0); 16 m = size(D(v),1);
17 default_arg('E',0); 17 default_arg('E',@(v)sparse(m,m));
18 18 default_arg('S',@(v,t)sparse(m,1));
19 if isnumeric(S) && S == 0
20 S = @(v)0;
21 end
22
23 if isnumeric(E) && E == 0
24 E = @(v)0;
25 end
26 19
27 20
28 % m = size(D,1); 21 % m = size(D,1);
29 % default_arg('E',sparse(m,m)); 22 % default_arg('E',sparse(m,m));
30 % default_arg('S',sparse(m,1)); 23 % default_arg('S',sparse(m,1));
47 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) 40 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
48 t = obj.t; 41 t = obj.t;
49 end 42 end
50 43
51 function obj = step(obj) 44 function obj = step(obj)
52 [obj.v, obj.v_prev] = time.cdiff.cdiff(obj.v, obj.v_prev, obj.k, obj.D(obj.v), obj.E(obj.v), obj.S(obj.v)); 45 D = obj.D(obj.v);
46 E = obj.E(obj.v);
47 S = obj.S(obj.v);
48
49 m = size(D,1);
50 I = speye(m);
51
52 %% Calculate for which indices we need to solve system of equations
53 [rows,cols] = find(E);
54 j = union(rows,cols);
55 i = setdiff(1:m,j);
56
57
58 %% Calculate matrices need for the timestep
59 % Before optimization: A = 1/k^2 * I - 1/(2*k)*E;
60 k = obj.k;
61 Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j);
62 B = 2/k^2 * I + D;
63 C = -1/k^2 * I - 1/(2*k)*E;
64
65 %% Take the timestep
66 v = obj.v;
67
68 % Before optimization: obj.v = A\(B*v + C*v_prev + S);
69 obj.v(i) = k^2*(B(i,i)*v(i) -1/k^2*obj.v_prev(i) + S(i));
70 obj.v(j) = Aj\(B(j,j)*v(j) + C(j,j)*obj.v_prev(j) + S(j));
71 obj.v_prev = v;
72
73 %% Update state of the timestepper
53 obj.t = obj.t + obj.k; 74 obj.t = obj.t + obj.k;
54 obj.n = obj.n + 1; 75 obj.n = obj.n + 1;
55 end 76 end
56 end 77 end
57 end 78 end