Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 365:f908ce064f35

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added SBP in time timestepper.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 24 Jan 2017 14:01:18 +0100
parents
children fccd746d8573 14f2be4fe9c1
comparison
equal deleted inserted replaced
363:8c8fd649bae3 365:f908ce064f35
1 classdef SBPInTime < time.Timestepper
2 % The SBP in time method.
3 % Implemented for v_t = A*v + f(t)
4 % k_local -- time-step
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
8 % k = k_local*(Nblock-1) in time.
9 % M -- matrix used in every solve.
10 % [L,U,P,Q] = lu(M);
11 properties
12 M
13 L
14 U
15 P
16 Q
17 A
18 Et_r
19 penalty
20 f
21 k_local
22 k
23 t
24 v
25 m
26 n
27 Nblock
28 order
29 nodes
30 end
31
32 methods
33 function obj = SBPInTime(A, f, k, order, Nblock, t0, v0, TYPE)
34 default_arg('TYPE','equidistant');
35 default_arg('Nblock',time.SBPInTime.smallestBlockSize(order,TYPE));
36
37 obj.A = A;
38 obj.f = f;
39 obj.k_local = k;
40 obj.k = k*(Nblock-1);
41 obj.Nblock = Nblock;
42 obj.t = t0;
43 obj.m = length(v0);
44 obj.n = 0;
45
46 %==== Build the time discretization matrix =====%
47 switch TYPE
48 case 'equidistant'
49 ops = sbp.D2Standard(Nblock,{0,obj.k},order);
50 case 'optimal'
51 ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order);
52 case 'minimal'
53 ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order,'minimal');
54 end
55
56 D1 = ops.D1;
57 HI = ops.HI;
58 e_l = ops.e_l;
59 e_r = ops.e_r;
60 obj.nodes = ops.x;
61
62 Ix = speye(size(A));
63 It = speye(Nblock,Nblock);
64
65 obj.Et_r = kron(e_r,Ix);
66
67 % Time derivative + penalty
68 tau = 1;
69 Mt = D1 + tau*HI*(e_l*e_l');
70
71 % penalty to impose "data"
72 penalty = tau*HI*e_l;
73 obj.penalty = kron(penalty,Ix);
74
75 Mx = kron(It,A);
76 Mt = kron(Mt,Ix);
77 obj.M = Mt - Mx;
78 %==============================================%
79
80 % LU factorization
81 [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M);
82
83 % Pretend that the initial condition is the last level
84 % of a previous step.
85 obj.v = obj.Et_r * v0;
86
87 end
88
89 function [v,t] = getV(obj)
90 v = obj.Et_r' * obj.v;
91 t = obj.t;
92 end
93
94 function obj = step(obj)
95 obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,...
96 obj.penalty, obj.f, obj.Nblock,...
97 obj.Et_r,...
98 obj.L, obj.U, obj.P, obj.Q);
99 obj.t = obj.t + obj.k;
100 obj.n = obj.n + obj.Nblock-1;
101 end
102 end
103
104
105 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)
125 default_arg('TYPE','equidistant')
126
127 switch TYPE
128
129 case 'equidistant'
130 switch order
131 case 2
132 N = 2;
133 case 4
134 N = 8;
135 case 6
136 N = 12;
137 case 8
138 N = 16;
139 case 10
140 N = 20;
141 case 12
142 N = 24;
143 otherwise
144 error('Operator does not exist');
145 end
146
147 case 'optimal'
148
149 switch order
150 case 4
151 N = 8;
152 case 6
153 N = 12;
154 case 8
155 N = 16;
156 case 10
157 N = 20;
158 case 12
159 N = 24;
160 otherwise
161 error('Operator does not exist');
162 end
163
164 case 'minimal'
165
166 switch order
167 case 4
168 N = 6;
169 case 6
170 N = 10;
171 case 8
172 N = 12;
173 case 10
174 N = 16;
175 case 12
176 N = 20;
177 otherwise
178 error('Operator does not exist');
179 end
180
181 end
182
183 end
184
185 end
186
187
188
189 end