Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcousticsVariable.m @ 652:be941bb0a11a feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add staggered 1D variable coefficient. Convergence study working.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 04 Dec 2017 11:11:03 -0800
parents
children 2351a7690e8a
comparison
equal deleted inserted replaced
651:4ee7d15bd8e6 652:be941bb0a11a
1 classdef Staggered1DAcousticsVariable < scheme.Scheme
2 properties
3 m % Number of points of primal grid in each direction, possibly a vector
4 h % Grid spacing
5
6 % Grids
7 grid % Total grid object
8 grid_primal
9 grid_dual
10
11 order % Order accuracy for the approximation
12
13 H % Combined norm
14 Hi % Inverse
15 D % Semi-discrete approximation matrix
16
17 D1_primal
18 D1_dual
19
20 % Pick out left or right boundary
21 e_l
22 e_r
23
24 % Initial data
25 v0
26
27 % Pick out primal or dual component
28 e_primal
29 e_dual
30
31 % System matrices
32 A
33 B
34
35 %
36 A_l, A_r
37 B_l, B_r
38 end
39
40
41 methods
42 % Scheme for A*u_t + B u_x = 0,
43 % u = [p, v];
44 % A: Diagonal and A > 0,
45 % A = [a1, 0;
46 % 0, a2]
47 %
48 % B = B^T and with diagonal entries = 0.
49 % B = [0 b
50 % b 0]
51 % Here we store p on the primal grid and v on the dual
52 function obj = Staggered1DAcousticsVariable(g, order, A, B)
53 default_arg('B',{@(x)0*x, @(x)0*x+1; @(x)0*x+1, @(x)0*x});
54 default_arg('A',{@(x)0*x+1, @(x)0*x; @(x)0*x, @(x)0*x+1});
55
56 obj.order = order;
57 obj.A = A;
58 obj.B = B;
59
60 % Grids
61 obj.m = g.size();
62 xl = g.getBoundary('l');
63 xr = g.getBoundary('r');
64 xlim = {xl, xr};
65
66 % Boundary matrices
67 obj.A_l = [A{1,1}(xl), A{1,2}(xl);....
68 A{2,1}(xl), A{2,2}(xl)];
69 obj.A_r = [A{1,1}(xr), A{1,2}(xr);....
70 A{2,1}(xr), A{2,2}(xr)];
71 obj.B_l = [B{1,1}(xl), B{1,2}(xl);....
72 B{2,1}(xl), B{2,2}(xl)];
73 obj.B_r = [B{1,1}(xr), B{1,2}(xr);....
74 B{2,1}(xr), B{2,2}(xr)];
75
76 obj.grid = g;
77 obj.grid_primal = g.grids{1};
78 obj.grid_dual = g.grids{2};
79 x_primal = obj.grid_primal.points();
80 x_dual = obj.grid_dual.points();
81
82 % Get operators
83 ops = sbp.D1StaggeredUpwind(obj.m, xlim, order);
84 obj.h = ops.h;
85
86 % Build combined operators
87 H_primal = ops.H_primal;
88 H_dual = ops.H_dual;
89 obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } );
90 obj.Hi = inv(obj.H);
91
92 D1_primal = ops.D1_primal;
93 D1_dual = ops.D1_dual;
94 A11_B12 = spdiag(-1./A{1,1}(x_primal).*B{1,2}(x_primal), 0);
95 A22_B21 = spdiag(-1./A{2,2}(x_dual).*B{2,1}(x_dual), 0);
96 D = {[], A11_B12*D1_primal;...
97 A22_B21*D1_dual, []};
98 obj.D = blockmatrix.toMatrix(D);
99 obj.D1_primal = D1_primal;
100 obj.D1_dual = D1_dual;
101
102 % Combined boundary operators
103 e_primal_l = ops.e_primal_l;
104 e_primal_r = ops.e_primal_r;
105 e_dual_l = ops.e_dual_l;
106 e_dual_r = ops.e_dual_r;
107 e_l = {e_primal_l, [];...
108 [] , e_dual_l};
109 e_r = {e_primal_r, [];...
110 [] , e_dual_r};
111 obj.e_l = blockmatrix.toMatrix(e_l);
112 obj.e_r = blockmatrix.toMatrix(e_r);
113
114 % Pick out first or second component of solution
115 N_primal = obj.grid_primal.N();
116 N_dual = obj.grid_dual.N();
117 obj.e_primal = [speye(N_primal, N_primal); sparse(N_dual, N_primal)];
118 obj.e_dual = [sparse(N_primal, N_dual); speye(N_dual, N_dual)];
119
120
121 end
122 % Closure functions return the operators applied to the own domain to close the boundary
123 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
124 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
125 % type is a string specifying the type of boundary condition if there are several.
126 % neighbour_scheme is an instance of Scheme that should be interfaced to.
127 % neighbour_boundary is a string specifying which boundary to interface to.
128 function [closure, penalty] = boundary_condition(obj, boundary, type)
129
130 default_arg('type','p');
131
132 % type = 'p' => boundary condition for p
133 % type = 'v' => boundary condition for v
134 % No other types implemented yet
135
136 % BC on the form Lu - g = 0;
137
138 % Get boundary matrices
139 switch boundary
140 case 'l'
141 A = full(obj.A_l);
142 B = full(obj.B_l);
143 case 'r'
144 A = full(obj.A_r);
145 B = full(obj.B_r);
146 end
147
148 % Diagonalize B
149 [T, Lambda] = eig(B);
150 lambda = diag(Lambda);
151
152 % Identify in- and outgoing characteristic variables
153 Iplus = lambda > 0;
154 Iminus = lambda < 0;
155
156 switch boundary
157 case 'l'
158 Iout = Iminus;
159 Iin = Iplus;
160 case 'r'
161 Iout = Iplus;
162 Iin = Iminus;
163 end
164
165 Tin = T(:,Iin);
166 Tout = T(:,Iout);
167
168 switch type
169 case 'p'
170 L = [1, 0];
171 case 'v'
172 L = [0, 1];
173 case 'characteristic'
174 L = Tin';
175 otherwise
176 error('Boundary condition not implemented.');
177 end
178
179 % Penalty parameters
180 sigma = [0; 0];
181 sigma(Iin) = lambda(Iin);
182 switch boundary
183 case 'l'
184 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin);
185 closure = obj.Hi*tau*L*obj.e_l';
186
187 case 'r'
188 tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin);
189 closure = obj.Hi*tau*L*obj.e_r';
190
191 end
192
193 penalty = -obj.Hi*tau;
194
195 end
196
197 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
198
199 error('Staggered1DAcoustics, interface not implemented');
200
201 switch boundary
202 % Upwind coupling
203 case {'l','left'}
204 tau = -1*obj.e_l;
205 closure = obj.Hi*tau*obj.e_l';
206 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
207 case {'r','right'}
208 tau = 0*obj.e_r;
209 closure = obj.Hi*tau*obj.e_r';
210 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
211 end
212
213 end
214
215 function N = size(obj)
216 N = obj.m;
217 end
218
219 end
220
221 methods(Static)
222 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
223 % and bound_v of scheme schm_v.
224 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
225 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
226 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
227 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
228 end
229 end
230 end