Mercurial > repos > public > sbplib

comparison +time/SBPInTime.m @ 423:a2cb0d4f4a02 feature/grids

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