Mercurial > repos > public > sbplib

comparison +scheme/Staggered1DAcoustics.m @ 642:8f0ee6ba6313 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme Staggered1DAcoustics
author Martin Almquist <malmquist@stanford.edu>
date Sat, 11 Nov 2017 20:30:13 -0800
parents
children dc2918fb104d
comparison
equal deleted inserted replaced
641:070d578997f6 642:8f0ee6ba6313
1 classdef Staggered1DAcoustics < 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 end
35
36
37 methods
38 % Scheme for A*u_t + B u_x = 0,
39 % u = [p, v];
40 % A: Diagonal and A > 0,
41 % A = [a1, 0;
42 % 0, a2]
43 %
44 % B = B^T and with diagonal entries = 0.
45 % B = [0 b
46 % b 0]
47 % Here we store p on the primal grid and v on the dual
48 function obj = Staggered1DAcoustics(g, order, A, B)
49 default_arg('A',[1, 0; 0, 1]);
50 default_arg('B',[0, 1; 1, 0]);
51
52 obj.order = order;
53 obj.A = A;
54 obj.B = B;
55
56 % Grids
57 obj.m = g.size();
58 xl = g.getBoundary('l');
59 xr = g.getBoundary('r');
60 xlim = {xl, xr};
61
62 obj.grid = g;
63 obj.grid_primal = g.grids{1};
64 obj.grid_dual = g.grids{2};
65
66 % Get operators
67 ops = sbp.D1Staggered(obj.m, xlim, order);
68 obj.h = ops.h;
69
70 % Build combined operators
71 H_primal = ops.H_primal;
72 H_dual = ops.H_dual;
73 obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } );
74 obj.Hi = inv(obj.H);
75
76 D1_primal = ops.D1_primal;
77 D1_dual = ops.D1_dual;
78 D = {[], -1/A(1,1)*B(1,2)*D1_primal;...
79 -1/A(2,2)*B(2,1)*D1_dual, []};
80 obj.D = blockmatrix.toMatrix(D);
81 obj.D1_primal = D1_primal;
82 obj.D1_dual = D1_dual;
83
84 % Combined boundary operators
85 e_primal_l = ops.e_primal_l;
86 e_primal_r = ops.e_primal_r;
87 e_dual_l = ops.e_dual_l;
88 e_dual_r = ops.e_dual_r;
89 e_l = {e_primal_l, [];...
90 [] , e_dual_l};
91 e_r = {e_primal_r, [];...
92 [] , e_dual_r};
93 obj.e_l = blockmatrix.toMatrix(e_l);
94 obj.e_r = blockmatrix.toMatrix(e_r);
95
96 % Pick out first or second component of solution
97 N_primal = obj.grid_primal.N();
98 N_dual = obj.grid_dual.N();
99 obj.e_primal = [speye(N_primal, N_primal); sparse(N_dual, N_primal)];
100 obj.e_dual = [sparse(N_primal, N_dual); speye(N_dual, N_dual)];
101
102
103 end
104 % Closure functions return the operators applied to the own domain to close the boundary
105 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain.
106 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
107 % type is a string specifying the type of boundary condition if there are several.
108 % neighbour_scheme is an instance of Scheme that should be interfaced to.
109 % neighbour_boundary is a string specifying which boundary to interface to.
110 function [closure, penalty] = boundary_condition(obj, boundary, type)
111
112 default_arg('type','p');
113
114 % type = 'p' => boundary condition for p
115 % type = 'v' => boundary condition for v
116 % No other types implemented yet
117
118 % BC on the form Lu - g = 0;
119 switch type
120 case 'p'
121 L = [1, 0];
122 case 'v'
123 L = [0, 1];
124 otherwise
125 error('Boundary condition not implemented.');
126 end
127
128 % Diagonalize B
129 B = obj.B;
130 [T, Lambda] = eig(B);
131 lambda = diag(Lambda);
132
133 % Identify in- and outgoing characteristic variables
134 Iplus = lambda > 0;
135 Iminus = lambda < 0;
136
137 switch boundary
138 case 'l'
139 Iout = Iminus;
140 Iin = Iplus;
141 case 'r'
142 Iout = Iplus;
143 Iin = Iminus;
144 end
145
146 Tin = T(:,Iin);
147 Tout = T(:,Iout);
148
149 % Penalty parameters
150 A = obj.A;
151 sigma = [0; 0];
152 sigma(Iin) = lambda(Iin);
153 switch boundary
154 case 'l'
155 tau = -1*obj.e_l * inv(A) * T * sigma * inv(L*Tin);
156 closure = obj.Hi*tau*L*obj.e_l';
157
158 case 'r'
159 tau = 1*obj.e_r * inv(A) * T * sigma * inv(L*Tin);
160 closure = obj.Hi*tau*L*obj.e_r';
161
162 end
163
164 penalty = -obj.Hi*tau;
165
166 end
167
168 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
169
170 error('Staggered1DAcoustics, interface not implemented');
171
172 switch boundary
173 % Upwind coupling
174 case {'l','left'}
175 tau = -1*obj.e_l;
176 closure = obj.Hi*tau*obj.e_l';
177 penalty = -obj.Hi*tau*neighbour_scheme.e_r';
178 case {'r','right'}
179 tau = 0*obj.e_r;
180 closure = obj.Hi*tau*obj.e_r';
181 penalty = -obj.Hi*tau*neighbour_scheme.e_l';
182 end
183
184 end
185
186 function N = size(obj)
187 N = obj.m;
188 end
189
190 end
191
192 methods(Static)
193 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
194 % and bound_v of scheme schm_v.
195 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
196 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
197 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
198 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
199 end
200 end
201 end