Mercurial > repos > public > sbplib
annotate +time/CdiffNonlin.m @ 886:8894e9c49e40 feature/timesteppers
Merge with default for latest changes
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Thu, 15 Nov 2018 16:36:21 -0800 |
| parents | b5e5b195da1e |
| children | f5e14e5986b5 |
| rev | line source |
|---|---|
|
1
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
1 classdef CdiffNonlin < time.Timestepper |
| 0 | 2 properties |
| 3 D | |
| 4 E | |
| 5 S | |
| 6 k | |
| 7 t | |
| 8 v | |
| 9 v_prev | |
| 10 n | |
| 11 end | |
| 12 | |
| 13 | |
| 14 methods | |
|
13
b18d3d201a71
Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
4
diff
changeset
|
15 function obj = CdiffNonlin(D, E, S, k, t0,n0, v, v_prev) |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
16 m = size(D(v),1); |
|
30
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
17 default_arg('E',0); |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
18 default_arg('S',0); |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
19 |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
20 if isnumeric(S) |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
21 S = @(v,t)S; |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
22 end |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
23 |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
24 if isnumeric(E) |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
25 E = @(v)E; |
|
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
26 end |
| 0 | 27 |
|
1
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
28 |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
29 % m = size(D,1); |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
30 % default_arg('E',sparse(m,m)); |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
31 % default_arg('S',sparse(m,1)); |
|
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
32 |
| 0 | 33 obj.D = D; |
| 34 obj.E = E; | |
| 35 obj.S = S; | |
| 36 obj.k = k; | |
|
1
5ae4f23d9130
Added CdiffNonlin timestepper. Probably fixed a bug with Cdiff. Added default arguments to Rk4SecondOrderNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
0
diff
changeset
|
37 obj.t = t0; |
|
13
b18d3d201a71
Fixed initialization of step counter in timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
4
diff
changeset
|
38 obj.n = n0; |
| 0 | 39 obj.v = v; |
| 40 obj.v_prev = v_prev; | |
| 41 end | |
| 42 | |
| 43 function [v,t] = getV(obj) | |
| 44 v = obj.v; | |
| 45 t = obj.t; | |
| 46 end | |
| 47 | |
| 48 function [vt,t] = getVt(obj) | |
| 49 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u)) | |
| 50 t = obj.t; | |
| 51 end | |
| 52 | |
|
816
b5e5b195da1e
Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
31
diff
changeset
|
53 function state = getState(obj) |
|
b5e5b195da1e
Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
31
diff
changeset
|
54 [v, t] = obj.getV(); |
|
b5e5b195da1e
Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
31
diff
changeset
|
55 [vt] = obj.getVt(); |
|
b5e5b195da1e
Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
31
diff
changeset
|
56 state = struct('v', v, 'vt', vt, 't', t, 'k', obj.k); |
|
b5e5b195da1e
Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
31
diff
changeset
|
57 end |
|
b5e5b195da1e
Add getState to timesteppers, returning the relevant state of the timestepper
Vidar Stiernström <vidar.stiernstrom@it.uu.se>
parents:
31
diff
changeset
|
58 |
| 0 | 59 function obj = step(obj) |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
60 D = obj.D(obj.v); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
61 E = obj.E(obj.v); |
|
30
cd2e28c5ecd2
Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents:
13
diff
changeset
|
62 S = obj.S(obj.v,obj.t); |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
63 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
64 m = size(D,1); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
65 I = speye(m); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
66 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
67 %% Calculate for which indices we need to solve system of equations |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
68 [rows,cols] = find(E); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
69 j = union(rows,cols); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
70 i = setdiff(1:m,j); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
71 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
72 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
73 %% Calculate matrices need for the timestep |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
74 % Before optimization: A = 1/k^2 * I - 1/(2*k)*E; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
75 k = obj.k; |
|
31
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
76 |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
77 Aj = 1/k^2 * I(j,j) - 1/(2*k)*E(j,j); |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
78 B = 2/k^2 * I + D; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
79 C = -1/k^2 * I - 1/(2*k)*E; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
80 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
81 %% Take the timestep |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
82 v = obj.v; |
|
31
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
83 v_prev = obj.v_prev; |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
84 |
|
31
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
85 % Want to solve the seq A*v_next = b where |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
86 b = (B*v + C*v_prev + S); |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
87 |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
88 % Before optimization: obj.v = A\b; |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
89 |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
90 obj.v(i) = k^2*b(i); |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
91 obj.v(j) = Aj\b(j); |
|
d1f9dd55a2b0
Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
30
diff
changeset
|
92 |
|
4
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
93 obj.v_prev = v; |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
94 |
|
8e14b5a577a6
Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents:
1
diff
changeset
|
95 %% Update state of the timestepper |
| 0 | 96 obj.t = obj.t + obj.k; |
| 97 obj.n = obj.n + 1; | |
| 98 end | |
| 99 end | |
| 100 end |