Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 400:14f2be4fe9c1 feature/beams

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