Mercurial > repos > public > sbplib

annotate +time/SBPInTimeScaled.m @ 1344:b4e5e45bd239 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Remove round off zeros from D2Nonequidistant operators
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sat, 15 Oct 2022 15:48:20 +0200
parents e95a0f2f7a8d
children 47e86b5270ad
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
746
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
1 classdef SBPInTimeScaled < time.Timestepper
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
2 % The SBP in time method.
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
3 % Implemented for A*v_t = B*v + f(t), v(0) = v0
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
4 % The resulting system of equations is
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
5 % M*u_next= K*u_prev_end + f
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
6 properties
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
7 A,B
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
8 f
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
9
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
10 k % total time step.
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
11
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
12 blockSize % number of points in each block
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
13 N % Number of components
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
14
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
15 order
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
16 nodes
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
17
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
18 Mtilde,Ktilde % System matrices
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
19 L,U,p,q % LU factorization of M
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
20 e_T
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
21
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
22 scaling
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
23 S, Sinv % Scaling matrices
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
24
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
25 % Time state
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
26 t
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
27 vtilde
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
28 n
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
29 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
30
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
31 methods
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
32 function obj = SBPInTimeScaled(A, B, f, k, t0, v0, scaling, TYPE, order, blockSize)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
33 default_arg('TYPE','gauss');
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
34 default_arg('f',[]);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
35
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
36 if(strcmp(TYPE,'gauss'))
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
37 default_arg('order',4)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
38 default_arg('blockSize',4)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
39 else
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
40 default_arg('order', 8);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
41 default_arg('blockSize',time.SBPInTimeImplicitFormulation.smallestBlockSize(order,TYPE));
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
42 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
43
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
44 obj.A = A;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
45 obj.B = B;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
46 obj.scaling = scaling;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
47
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
48 if ~isempty(f)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
49 obj.f = f;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
50 else
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
51 obj.f = @(t)sparse(length(v0),1);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
52 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
53
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
54 obj.k = k;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
55 obj.blockSize = blockSize;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
56 obj.N = length(v0);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
57
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
58 obj.n = 0;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
59 obj.t = t0;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
60
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
61 %==== Build the time discretization matrix =====%
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
62 switch TYPE
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
63 case 'equidistant'
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
64 ops = sbp.D2Standard(blockSize,{0,obj.k},order);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
65 case 'optimal'
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
66 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
67 case 'minimal'
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
68 ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal');
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
69 case 'gauss'
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
70 ops = sbp.D1Gauss(blockSize,{0,obj.k});
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
71 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
72
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
73 I = speye(size(A));
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
74 I_t = speye(blockSize,blockSize);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
75
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
76 D1 = kron(ops.D1, I);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
77 HI = kron(ops.HI, I);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
78 e_0 = kron(ops.e_l, I);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
79 e_T = kron(ops.e_r, I);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
80 obj.nodes = ops.x;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
81
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
82 % Convert to form M*w = K*v0 + f(t)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
83 tau = kron(I_t, A) * e_0;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
84 M = kron(I_t, A)*D1 + HI*tau*e_0' - kron(I_t, B);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
85
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
86 K = HI*tau;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
87
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
88 obj.S = kron(I_t, spdiag(scaling));
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
89 obj.Sinv = kron(I_t, spdiag(1./scaling));
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
90
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
91 obj.Mtilde = obj.Sinv*M*obj.S;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
92 obj.Ktilde = obj.Sinv*K*spdiag(scaling);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
93 obj.e_T = e_T;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
94
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
95
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
96 % LU factorization
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
97 [obj.L,obj.U,obj.p,obj.q] = lu(obj.Mtilde, 'vector');
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
98
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
99 obj.vtilde = (1./obj.scaling).*v0;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
100 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
101
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
102 function [v,t] = getV(obj)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
103 v = obj.scaling.*obj.vtilde;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
104 t = obj.t;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
105 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
106
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
107 function obj = step(obj)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
108 forcing = zeros(obj.blockSize*obj.N,1);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
109
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
110 for i = 1:obj.blockSize
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
111 forcing((1 + (i-1)*obj.N):(i*obj.N)) = obj.f(obj.t + obj.nodes(i));
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
112 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
113
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
114 RHS = obj.Sinv*forcing + obj.Ktilde*obj.vtilde;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
115
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
116 y = obj.L\RHS(obj.p);
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
117 z = obj.U\y;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
118
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
119 w = zeros(size(z));
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
120 w(obj.q) = z;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
121
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
122 obj.vtilde = obj.e_T'*w;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
123
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
124 obj.t = obj.t + obj.k;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
125 obj.n = obj.n + 1;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
126 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
127 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
128
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
129 methods(Static)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
130 function N = smallestBlockSize(order,TYPE)
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
131 default_arg('TYPE','gauss')
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
132
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
133 switch TYPE
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
134 case 'gauss'
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
135 N = 4;
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
136 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
137 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
138 end
e95a0f2f7a8d Add file that was forgotten.
Jonatan Werpers <jonatan@werpers.com>
parents:
diff changeset
139 end