Mercurial > repos > public > sbplib

annotate +time/CdiffNonlin.m @ 1142:cff49fba3cc8 rv-interpolation

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Closing branch
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Mon, 05 Aug 2019 10:49:21 +0200
parents d1f9dd55a2b0
children b5e5b195da1e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 D
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 E
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 S
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 k
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 t
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 v
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 v_prev
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 n
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 obj.D = D;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 obj.E = E;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35 obj.S = S;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
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
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 obj.v = v;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 obj.v_prev = v_prev;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43 function [v,t] = getV(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 v = obj.v;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 function [vt,t] = getVt(obj)
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 t = obj.t;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53 function obj = step(obj)
4
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
54 D = obj.D(obj.v);
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
55 E = obj.E(obj.v);
30
cd2e28c5ecd2 Corrected parameter treatment in nonlin timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents: 13
diff changeset
56 S = obj.S(obj.v,obj.t);
4
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
57
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
58 m = size(D,1);
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
59 I = speye(m);
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
60
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
61 %% 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
62 [rows,cols] = find(E);
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
63 j = union(rows,cols);
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
64 i = setdiff(1:m,j);
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
65
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 matrices need for the timestep
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
68 % 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
69 k = obj.k;
31
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
70
4
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
71 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
72 B = 2/k^2 * I + D;
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
73 C = -1/k^2 * I - 1/(2*k)*E;
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
74
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
75 %% Take the timestep
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
76 v = obj.v;
31
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
77 v_prev = obj.v_prev;
4
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
78
31
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
79 % 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
80 b = (B*v + C*v_prev + S);
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
81
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
82 % Before optimization: obj.v = A\b;
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
83
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
84 obj.v(i) = k^2*b(i);
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
85 obj.v(j) = Aj\b(j);
d1f9dd55a2b0 Fixed bug in step() for CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 30
diff changeset
86
4
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
87 obj.v_prev = v;
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
88
8e14b5a577a6 Attempt att optimization of CdiffNonlin.
Jonatan Werpers <jonatan@werpers.com>
parents: 1
diff changeset
89 %% Update state of the timestepper
0
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90 obj.t = obj.t + obj.k;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 obj.n = obj.n + 1;
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 end
48b6fb693025 Initial commit.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94 end