Mercurial > repos > public > sbplib

comparison +time/SBPInTimeImplicitFormulation.m @ 886:8894e9c49e40 feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default for latest changes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 15 Nov 2018 16:36:21 -0800
parents 5df7f99206b2
children 47e86b5270ad
comparison
equal deleted inserted replaced
816:b5e5b195da1e 886:8894e9c49e40
1 classdef SBPInTimeImplicitFormulation < time.Timestepper
2 % The SBP in time method.
3 % Implemented for A*v_t = B*v + f(t), v(0) = v0
4 properties
5 A,B
6 f
7
8 k % total time step.
9
10 blockSize % number of points in each block
11 N % Number of components
12
13 order
14 nodes
15
16 M,K % System matrices
17 L,U,p,q % LU factorization of M
18 e_T
19
20 % Time state
21 t
22 v
23 n
24 end
25
26 methods
27 function obj = SBPInTimeImplicitFormulation(A, B, f, k, t0, v0, TYPE, order, blockSize)
28
29 default_arg('TYPE','gauss');
30 default_arg('f',[]);
31
32 if(strcmp(TYPE,'gauss'))
33 default_arg('order',4)
34 default_arg('blockSize',4)
35 else
36 default_arg('order', 8);
37 default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE));
38 end
39
40 obj.A = A;
41 obj.B = B;
42
43 if ~isempty(f)
44 obj.f = f;
45 else
46 obj.f = @(t)sparse(length(v0),1);
47 end
48
49 obj.k = k;
50 obj.blockSize = blockSize;
51 obj.N = length(v0);
52
53 obj.n = 0;
54 obj.t = t0;
55
56 %==== Build the time discretization matrix =====%
57 switch TYPE
58 case 'equidistant'
59 ops = sbp.D2Standard(blockSize,{0,obj.k},order);
60 case 'optimal'
61 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
62 case 'minimal'
63 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
64 case 'gauss'
65 ops = sbp.D1Gauss(blockSize,{0,obj.k});
66 end
67
68 I = speye(size(A));
69 I_t = speye(blockSize,blockSize);
70
71 D1 = kron(ops.D1, I);
72 HI = kron(ops.HI, I);
73 e_0 = kron(ops.e_l, I);
74 e_T = kron(ops.e_r, I);
75 obj.nodes = ops.x;
76
77 % Convert to form M*w = K*v0 + f(t)
78 tau = kron(I_t, A) * e_0;
79 M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B);
80
81 K = HI*tau;
82
83 obj.M = M;
84 obj.K = K;
85 obj.e_T = e_T;
86
87 % LU factorization
88 [obj.L,obj.U,obj.p,obj.q] = lu(obj.M, 'vector');
89
90 obj.v = v0;
91 end
92
93 function [v,t] = getV(obj)
94 v = obj.v;
95 t = obj.t;
96 end
97
98 function obj = step(obj)
99 RHS = zeros(obj.blockSize*obj.N,1);
100
101 for i = 1:obj.blockSize
102 RHS((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i));
103 end
104
105 RHS = RHS + obj.K*obj.v;
106
107 y = obj.L\RHS(obj.p);
108 z = obj.U\y;
109
110 w = zeros(size(z));
111 w(obj.q) = z;
112
113 obj.v = obj.e_T'*w;
114
115 obj.t = obj.t + obj.k;
116 obj.n = obj.n + 1;
117 end
118 end
119
120 methods(Static)
121 function N = smallestBlockSize(order,TYPE)
122 default_arg('TYPE','gauss')
123
124 switch TYPE
125 case 'gauss'
126 N = 4;
127 end
128 end
129 end
130 end