Mercurial > repos > public > sbplib

annotate +time/SBPInTimeScaled.m @ 1198:2924b3a9b921 feature/d2_compatible

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add OpSet for fully compatible D2Variable, created from regular D2Variable by replacing d1 by first row of D1. Formal reduction by one order of accuracy at the boundary point.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 16 Aug 2019 14:30:28 -0700
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