Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 398:fccd746d8573 feature/SBPinTimeClarity

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor for clarity
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 02 Feb 2017 09:57:54 +0100
parents f908ce064f35
children 9ff24a14f9ef 9fd9b1bea3d2
comparison
equal deleted inserted replaced
397:4fd5bfe5d0bb 398:fccd746d8573
1 classdef SBPInTime < time.Timestepper 1 classdef SBPInTime < time.Timestepper
2 % The SBP in time method. 2 % The SBP in time method.
3 % Implemented for v_t = A*v + f(t) 3 % Implemented for v_t = A*v + f(t)
4 % k_local -- time-step 4 %
5 % Nblock -- number of points in each block 5 % Each "step" takes one block step and thus advances
6 % nodes -- points such that t_n + nodes are the points in block n. 6 % k = k_local*(blockSize-1) in time.
7 % Each "step" takes one block step and thus advances
8 % k = k_local*(Nblock-1) in time.
9 % M -- matrix used in every solve.
10 % [L,U,P,Q] = lu(M);
11 properties 7 properties
12 M 8 M % System matrix
13 L 9 L,U,P,Q % LU factorization of M
14 U
15 P
16 Q
17 A 10 A
18 Et_r 11 Et_r
19 penalty 12 penalty
20 f 13 f
21 k_local 14 k_local % step size within a block
22 k 15 k % Time size of a block k/(blockSize-1) = k_local
23 t 16 t
24 v 17 v
25 m 18 m
26 n 19 n
27 Nblock 20 blockSize % number of points in each block
28 order 21 order
29 nodes 22 nodes
30 end 23 end
31 24
32 methods 25 methods
33 function obj = SBPInTime(A, f, k, order, Nblock, t0, v0, TYPE) 26 function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize)
34 default_arg('TYPE','equidistant'); 27 default_arg('TYPE','minimal');
35 default_arg('Nblock',time.SBPInTime.smallestBlockSize(order,TYPE)); 28 default_arg('order', 8);
36 29 default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE));
30
37 obj.A = A; 31 obj.A = A;
38 obj.f = f; 32 obj.f = f;
39 obj.k_local = k; 33 obj.k_local = k/(blockSize-1);
40 obj.k = k*(Nblock-1); 34 obj.k = k;
41 obj.Nblock = Nblock; 35 obj.blockSize = blockSize;
42 obj.t = t0; 36 obj.t = t0;
43 obj.m = length(v0); 37 obj.m = length(v0);
44 obj.n = 0; 38 obj.n = 0;
45 39
46 %==== Build the time discretization matrix =====% 40 %==== Build the time discretization matrix =====%
47 switch TYPE 41 switch TYPE
48 case 'equidistant' 42 case 'equidistant'
49 ops = sbp.D2Standard(Nblock,{0,obj.k},order); 43 ops = sbp.D2Standard(blockSize,{0,obj.k},order);
50 case 'optimal' 44 case 'optimal'
51 ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order); 45 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
52 case 'minimal' 46 case 'minimal'
53 ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order,'minimal'); 47 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
54 end 48 end
55 49
56 D1 = ops.D1; 50 D1 = ops.D1;
57 HI = ops.HI; 51 HI = ops.HI;
58 e_l = ops.e_l; 52 e_l = ops.e_l;
59 e_r = ops.e_r; 53 e_r = ops.e_r;
60 obj.nodes = ops.x; 54 obj.nodes = ops.x;
61 55
62 Ix = speye(size(A)); 56 Ix = speye(size(A));
63 It = speye(Nblock,Nblock); 57 It = speye(blockSize,blockSize);
64 58
65 obj.Et_r = kron(e_r,Ix); 59 obj.Et_r = kron(e_r,Ix);
66 60
67 % Time derivative + penalty 61 % Time derivative + penalty
68 tau = 1; 62 tau = 1;
69 Mt = D1 + tau*HI*(e_l*e_l'); 63 Mt = D1 + tau*HI*(e_l*e_l');
70 64
71 % penalty to impose "data" 65 % penalty to impose "data"
72 penalty = tau*HI*e_l; 66 penalty = tau*HI*e_l;
73 obj.penalty = kron(penalty,Ix); 67 obj.penalty = kron(penalty,Ix);
74 68
75 Mx = kron(It,A); 69 Mx = kron(It,A);
76 Mt = kron(Mt,Ix); 70 Mt = kron(Mt,Ix);
77 obj.M = Mt - Mx; 71 obj.M = Mt - Mx;
78 %==============================================% 72 %==============================================%
79 73
80 % LU factorization 74 % LU factorization
81 [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M); 75 [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M);
82 76
83 % Pretend that the initial condition is the last level 77 % Pretend that the initial condition is the last level
84 % of a previous step. 78 % of a previous step.
85 obj.v = obj.Et_r * v0; 79 obj.v = obj.Et_r * v0;
86 80
87 end 81 end
88 82
89 function [v,t] = getV(obj) 83 function [v,t] = getV(obj)
90 v = obj.Et_r' * obj.v; 84 v = obj.Et_r' * obj.v;
91 t = obj.t; 85 t = obj.t;
92 end 86 end
93 87
94 function obj = step(obj) 88 function obj = step(obj)
95 obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,... 89 obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,...
96 obj.penalty, obj.f, obj.Nblock,... 90 obj.penalty, obj.f, obj.blockSize,...
97 obj.Et_r,... 91 obj.Et_r,...
98 obj.L, obj.U, obj.P, obj.Q); 92 obj.L, obj.U, obj.P, obj.Q);
99 obj.t = obj.t + obj.k; 93 obj.t = obj.t + obj.k;
100 obj.n = obj.n + obj.Nblock-1; 94 obj.n = obj.n + 1;
101 end 95 end
102 end 96 end
103 97
104
105 methods(Static) 98 methods(Static)
106
107 %
108 function [k,numberOfBlocks] = alignedTimeStep(k,Tend,Nblock)
109
110 % input k is the desired time-step
111 % Nblock is the number of points per block.
112
113 % Make sure that we reach the final time by advancing
114 % an integer number of blocks
115 kblock = (Nblock-1)*k;
116 numberOfBlocks = ceil(Tend/kblock);
117 kblock = Tend/(numberOfBlocks);
118
119 % Corrected time step
120 k = kblock/(Nblock-1);
121
122 end
123
124 function N = smallestBlockSize(order,TYPE) 99 function N = smallestBlockSize(order,TYPE)
125 default_arg('TYPE','equidistant') 100 default_arg('TYPE','equidistant')
126 101
127 switch TYPE 102 switch TYPE
128 103
129 case 'equidistant' 104 case 'equidistant'
130 switch order 105 switch order
131 case 2 106 case 2
132 N = 2; 107 N = 2;
133 case 4 108 case 4
141 case 12 116 case 12
142 N = 24; 117 N = 24;
143 otherwise 118 otherwise
144 error('Operator does not exist'); 119 error('Operator does not exist');
145 end 120 end
146 121
147 case 'optimal' 122 case 'optimal'
148 123
149 switch order 124 switch order
150 case 4 125 case 4
151 N = 8; 126 N = 8;
152 case 6 127 case 6
153 N = 12; 128 N = 12;
158 case 12 133 case 12
159 N = 24; 134 N = 24;
160 otherwise 135 otherwise
161 error('Operator does not exist'); 136 error('Operator does not exist');
162 end 137 end
163 138
164 case 'minimal' 139 case 'minimal'
165 140
166 switch order 141 switch order
167 case 4 142 case 4
168 N = 6; 143 N = 6;
169 case 6 144 case 6
170 N = 10; 145 N = 10;
175 case 12 150 case 12
176 N = 20; 151 N = 20;
177 otherwise 152 otherwise
178 error('Operator does not exist'); 153 error('Operator does not exist');
179 end 154 end
180
181 end 155 end
182
183 end 156 end
184
185 end 157 end
186
187
188
189 end 158 end