Mercurial > repos > public > sbplib

annotate +time/CdiffImplicit.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
parents 1a30dbe99c7c
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