Mercurial > repos > public > sbplib

comparison +time/CdiffNonlin.m @ 1113:47e86b5270ad feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change name of property k to dt in time.Timestepper
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 Apr 2019 22:40:55 +0200
parents f5e14e5986b5
children
comparison
equal deleted inserted replaced
1112:835c8fa456ec 1113:47e86b5270ad
1 classdef CdiffNonlin < time.Timestepper 1 classdef CdiffNonlin < time.Timestepper
2 properties 2 properties
3 D 3 D
4 E 4 E
5 S 5 S
6 k 6 dt
7 t 7 t
8 v 8 v
9 v_prev 9 v_prev
10 n 10 n
11 end 11 end
12 12
13 13
14 methods 14 methods
15 function obj = CdiffNonlin(D, E, S, k, t0,n0, v, v_prev) 15 function obj = CdiffNonlin(D, E, S, dt, t0,n0, v, v_prev)
16 m = size(D(v),1); 16 m = size(D(v),1);
17 default_arg('E',0); 17 default_arg('E',0);
18 default_arg('S',0); 18 default_arg('S',0);
19 19
20 if isnumeric(S) 20 if isnumeric(S)
31 % default_arg('S',sparse(m,1)); 31 % default_arg('S',sparse(m,1));
32 32
33 obj.D = D; 33 obj.D = D;
34 obj.E = E; 34 obj.E = E;
35 obj.S = S; 35 obj.S = S;
36 obj.k = k; 36 obj.dt = dt;
37 obj.t = t0; 37 obj.t = t0;
38 obj.n = n0; 38 obj.n = n0;
39 obj.v = v; 39 obj.v = v;
40 obj.v_prev = v_prev; 40 obj.v_prev = v_prev;
41 end 41 end
44 v = obj.v; 44 v = obj.v;
45 t = obj.t; 45 t = obj.t;
46 end 46 end
47 47
48 function [vt,t] = getVt(obj) 48 function [vt,t] = getVt(obj)
49 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) 49 vt = (obj.v-obj.v_prev)/obj.dt; % Could be improved using u_tt = f(u))
50 t = obj.t; 50 t = obj.t;
51 end 51 end
52 52
53 function obj = step(obj) 53 function obj = step(obj)
54 D = obj.D(obj.v); 54 D = obj.D(obj.v);
63 j = union(rows,cols); 63 j = union(rows,cols);
64 i = setdiff(1:m,j); 64 i = setdiff(1:m,j);
65 65
66 66
67 %% Calculate matrices need for the timestep 67 %% Calculate matrices need for the timestep
68 % Before optimization: A = 1/k^2 * I - 1/(2*k)*E; 68 % Before optimization: A = 1/dt^2 * I - 1/(2*dt)*E;
69 k = obj.k; 69 dt = obj.dt;
70 70
71 Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j); 71 Aj = 1/dt^2 * I(j,j) - 1/(2*dt)*E(j,j);
72 B = 2/k^2 * I + D; 72 B = 2/dt^2 * I + D;
73 C = -1/k^2 * I - 1/(2*k)*E; 73 C = -1/dt^2 * I - 1/(2*dt)*E;
74 74
75 %% Take the timestep 75 %% Take the timestep
76 v = obj.v; 76 v = obj.v;
77 v_prev = obj.v_prev; 77 v_prev = obj.v_prev;
78 78
79 % Want to solve the seq A*v_next = b where 79 % Want to solve the seq A*v_next = b where
80 b = (B*v + C*v_prev + S); 80 b = (B*v + C*v_prev + S);
81 81
82 % Before optimization: obj.v = A\b; 82 % Before optimization: obj.v = A\b;
83 83
84 obj.v(i) = k^2*b(i); 84 obj.v(i) = dt^2*b(i);
85 obj.v(j) = Aj\b(j); 85 obj.v(j) = Aj\b(j);
86 86
87 obj.v_prev = v; 87 obj.v_prev = v;
88 88
89 %% Update state of the timestepper 89 %% Update state of the timestepper
90 obj.t = obj.t + obj.k; 90 obj.t = obj.t + obj.dt;
91 obj.n = obj.n + 1; 91 obj.n = obj.n + 1;
92 end 92 end
93 end 93 end
94 end 94 end