Mercurial > repos > public > sbplib

annotate +time/CdiffImplicit.m @ 978:1a30dbe99c7c

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor CdiffImplicit to take input arguments in the right order
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 07 Jan 2019 16:15:49 +0100
parents ea8fefc326b2
children 47e86b5270ad
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef CdiffImplicit < time.Timestepper
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 properties
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
3 A, B, C
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 AA, BB, CC
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
5 G
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 k
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 t
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 v, v_prev
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9 n
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11 % LU factorization
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 L,U,p,q
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 methods
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 % Solves
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
17 % A*v_tt + B*v_t + C*v = G(t)
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
18 % v(t0) = v0
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
19 % v_t(t0) = v0t
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
20 % starting at time t0 with timestep
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
21 % Using
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
22 % A*Dp*Dm*v_n + B*D0*v_n + C*I0*v_n = G(t_n)
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
23 function obj = CdiffImplicit(A, B, C, G, v0, v0t, k, t0)
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
24 m = length(v0);
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
25 default_arg('A', speye(m));
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
26 default_arg('B', sparse(m,m));
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
27 default_arg('G', @(t) sparse(m,1));
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 default_arg('t0', 0);
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30 obj.A = A;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 obj.B = B;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 obj.C = C;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 obj.G = G;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
35 % Rewrite as AA*v_(n+1) + BB*v_n + CC*v_(n-1) = G(t_n)
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
36 AA = A/k^2 + B/(2*k) + C/2;
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
37 BB = -2*A/k^2;
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
38 CC = A/k^2 - B/(2*k) + C/2;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 obj.AA = AA;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 obj.BB = BB;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 obj.CC = CC;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
44 v_prev = v0;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 I = speye(m);
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
46 v = v0 + k*v0t;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 if ~issparse(A) || ~issparse(B) || ~issparse(C)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 error('LU factorization with full pivoting only works for sparse matrices.')
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 [L,U,p,q] = lu(AA,'vector');
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 obj.L = L;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 obj.U = U;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 obj.p = p;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57 obj.q = q;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 obj.k = k;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60 obj.t = t0+k;
385
b361a04894e3 Fix bug in implicit cdiff. Add stepTo method for timesteppers.
Jonatan Werpers <jonatan@werpers.com>
parents: 383
diff changeset
61 obj.n = 1;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 obj.v = v;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 obj.v_prev = v_prev;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 function [v,t] = getV(obj)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67 v = obj.v;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 t = obj.t;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 function [vt,t] = getVt(obj)
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
72 vt = (obj.v-obj.v_prev)/obj.k; % Could be improved using u_tt = f(u))
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 t = obj.t;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75
976
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
76 % Calculate the conserved energy (Dm*v_n)^2_A + Im*v_n^2_B
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
77 function E = getEnergy(obj)
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
78 v = obj.v;
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
79 vp = obj.v_prev;
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
80 vt = (obj.v - obj.v_prev)/obj.k;
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
81
978
1a30dbe99c7c Refactor CdiffImplicit to take input arguments in the right order
Jonatan Werpers <jonatan@werpers.com>
parents: 976
diff changeset
82 E = vt'*obj.A*vt + 1/2*(v'*obj.C*v + vp'*obj.C*vp);
976
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
83 end
ea8fefc326b2 Add getEnergy methods to CdiffImplicit
Jonatan Werpers <jonatan@werpers.com>
parents: 956
diff changeset
84
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85 function obj = step(obj)
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 b = obj.G(obj.t) - obj.BB*obj.v - obj.CC*obj.v_prev;
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
87 obj.v_prev = obj.v;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
89 % % Backslash
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
90 % obj.v = obj.AA\b;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
92 % LU with column pivot
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
93 y = obj.L\b(obj.p);
383
151b08366d63 Fix bug in CdiffImplicit.
Jonatan Werpers <jonatan@werpers.com>
parents: 380
diff changeset
94 z = obj.U\y;
151b08366d63 Fix bug in CdiffImplicit.
Jonatan Werpers <jonatan@werpers.com>
parents: 380
diff changeset
95 obj.v(obj.q) = z;
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96
380
280ae4dbf93b CdiffImplicit: Fixed a bug. Switched to LU.
Jonatan Werpers <jonatan@werpers.com>
parents: 379
diff changeset
97 % Update time
379
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98 obj.t = obj.t + obj.k;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 obj.n = obj.n + 1;
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101 end
ca73ee0623e5 Added an implicit central time stepping scheme.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 end