Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 401:9ff24a14f9ef feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge clarity changes for SBPInTime.
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 02 Feb 2017 10:07:49 +0100
parents 14f2be4fe9c1 fccd746d8573
children 38173ea263ed
comparison
equal deleted inserted replaced
400:14f2be4fe9c1 401:9ff24a14f9ef
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
6 % nodes -- points such that t_n + nodes are the points in block n.
7 % Each "step" takes one block step and thus advances 5 % Each "step" takes one block step and thus advances
8 % k = k_local*(Nblock-1) in time. 6 % k = k_local*(blockSize-1) in time.
9 % M -- matrix used in every solve.
10 % [L,U,P,Q] = lu(M);
11 properties 7 properties
12 M % System matrix 8 M % System matrix
13 L,U,P % LU factorization of M 9 L,U,P,Q % LU factorization of M
14 Q
15 A 10 A
16 Et_r 11 Et_r
17 penalty 12 penalty
18 f 13 f
19 k_local 14 k_local % step size within a block
20 k 15 k % Time size of a block k/(blockSize-1) = k_local
21 t 16 t
22 v 17 v
23 m 18 m
24 n 19 n
25 Nblock 20 blockSize % number of points in each block
26 order 21 order
27 nodes 22 nodes
28 end 23 end
29 24
30 methods 25 methods
31 function obj = SBPInTime(A, f, k, order, Nblock, t0, v0, TYPE) 26 function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize)
32 default_arg('TYPE','equidistant'); 27 default_arg('TYPE','minimal');
33 default_arg('Nblock',time.SBPInTime.smallestBlockSize(order,TYPE)); 28 default_arg('order', 8);
29 default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE));
34 30
35 obj.A = A; 31 obj.A = A;
36 obj.f = f; 32 obj.f = f;
37 obj.k_local = k; 33 obj.k_local = k/(blockSize-1);
38 obj.k = k*(Nblock-1); 34 obj.k = k;
39 obj.Nblock = Nblock; 35 obj.blockSize = blockSize;
40 obj.t = t0; 36 obj.t = t0;
41 obj.m = length(v0); 37 obj.m = length(v0);
42 obj.n = 0; 38 obj.n = 0;
43 39
44 %==== Build the time discretization matrix =====% 40 %==== Build the time discretization matrix =====%
45 switch TYPE 41 switch TYPE
46 case 'equidistant' 42 case 'equidistant'
47 ops = sbp.D2Standard(Nblock,{0,obj.k},order); 43 ops = sbp.D2Standard(blockSize,{0,obj.k},order);
48 case 'optimal' 44 case 'optimal'
49 ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order); 45 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
50 case 'minimal' 46 case 'minimal'
51 ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order,'minimal'); 47 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
52 end 48 end
53 49
54 D1 = ops.D1; 50 D1 = ops.D1;
55 HI = ops.HI; 51 HI = ops.HI;
56 e_l = ops.e_l; 52 e_l = ops.e_l;
57 e_r = ops.e_r; 53 e_r = ops.e_r;
58 obj.nodes = ops.x; 54 obj.nodes = ops.x;
59 55
60 Ix = speye(size(A)); 56 Ix = speye(size(A));
61 It = speye(Nblock,Nblock); 57 It = speye(blockSize,blockSize);
62 58
63 obj.Et_r = kron(e_r,Ix); 59 obj.Et_r = kron(e_r,Ix);
64 60
65 % Time derivative + penalty 61 % Time derivative + penalty
66 tau = 1; 62 tau = 1;
89 t = obj.t; 85 t = obj.t;
90 end 86 end
91 87
92 function obj = step(obj) 88 function obj = step(obj)
93 obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,... 89 obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,...
94 obj.penalty, obj.f, obj.Nblock,... 90 obj.penalty, obj.f, obj.blockSize,...
95 obj.Et_r,... 91 obj.Et_r,...
96 obj.L, obj.U, obj.P, obj.Q); 92 obj.L, obj.U, obj.P, obj.Q);
97 obj.t = obj.t + obj.k; 93 obj.t = obj.t + obj.k;
98 obj.n = obj.n + obj.Nblock-1; 94 obj.n = obj.n + 1;
99 end 95 end
100 end 96 end
101 97
102 98
103 methods(Static) 99 methods(Static)
104
105 %
106 function [k,numberOfBlocks] = alignedTimeStep(k,Tend,Nblock)
107
108 % input k is the desired time-step
109 % Nblock is the number of points per block.
110
111 % Make sure that we reach the final time by advancing
112 % an integer number of blocks
113 kblock = (Nblock-1)*k;
114 numberOfBlocks = ceil(Tend/kblock);
115 kblock = Tend/(numberOfBlocks);
116
117 % Corrected time step
118 k = kblock/(Nblock-1);
119
120 end
121
122 function N = smallestBlockSize(order,TYPE) 100 function N = smallestBlockSize(order,TYPE)
123 default_arg('TYPE','equidistant') 101 default_arg('TYPE','equidistant')
124 102
125 switch TYPE 103 switch TYPE
126 104
173 case 12 151 case 12
174 N = 20; 152 N = 20;
175 otherwise 153 otherwise
176 error('Operator does not exist'); 154 error('Operator does not exist');
177 end 155 end
178
179 end 156 end
180
181 end 157 end
182
183 end 158 end
184
185
186
187 end 159 end