Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 756:f891758ad7a4 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Jun 2018 14:30:45 -0700
parents f4595f14d696
children 459eeb99130f
comparison
equal deleted inserted replaced
755:14f0058356f2 756:f891758ad7a4
1 classdef Utux2D < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5 grid % Grid
6 order % Order accuracy for the approximation
7 v0 % Initial data
8
9 a % Wave speed a = [a1, a2];
10 % Can either be a constant vector or a cell array of function handles.
11
12 H % Discrete norm
13 H_x, H_y % Norms in the x and y directions
14 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
15
16 % Derivatives
17 Dx, Dy
18
19 % Boundary operators
20 e_w, e_e, e_s, e_n
21
22 D % Total discrete operator
23
24 % String, type of interface coupling
25 % Default: 'upwind'
26 % Other: 'centered'
27 coupling_type
28
29 % String, type of interpolation operators
30 % Default: 'AWW' (Almquist Wang Werpers)
31 % Other: 'MC' (Mattsson Carpenter)
32 interpolation_type
33
34
35 % Cell array, damping on upwstream and downstream sides.
36 interpolation_damping
37
38 end
39
40
41 methods
42 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping)
43
44 default_arg('interpolation_damping',{0,0});
45 default_arg('interpolation_type','AWW');
46 default_arg('coupling_type','upwind');
47 default_arg('a',1/sqrt(2)*[1, 1]);
48 default_arg('opSet',@sbp.D2Standard);
49
50 assert(isa(g, 'grid.Cartesian'))
51 if iscell(a)
52 a1 = grid.evalOn(g, a{1});
53 a2 = grid.evalOn(g, a{2});
54 a = {spdiag(a1), spdiag(a2)};
55 else
56 a = {a(1), a(2)};
57 end
58
59 m = g.size();
60 m_x = m(1);
61 m_y = m(2);
62 m_tot = g.N();
63
64 xlim = {g.x{1}(1), g.x{1}(end)};
65 ylim = {g.x{2}(1), g.x{2}(end)};
66 obj.grid = g;
67
68 % Operator sets
69 ops_x = opSet(m_x, xlim, order);
70 ops_y = opSet(m_y, ylim, order);
71 Ix = speye(m_x);
72 Iy = speye(m_y);
73
74 % Norms
75 Hx = ops_x.H;
76 Hy = ops_y.H;
77 Hxi = ops_x.HI;
78 Hyi = ops_y.HI;
79
80 obj.H_x = Hx;
81 obj.H_y = Hy;
82 obj.H = kron(Hx,Hy);
83 obj.Hi = kron(Hxi,Hyi);
84 obj.Hx = kron(Hx,Iy);
85 obj.Hy = kron(Ix,Hy);
86 obj.Hxi = kron(Hxi,Iy);
87 obj.Hyi = kron(Ix,Hyi);
88
89 % Derivatives
90 Dx = ops_x.D1;
91 Dy = ops_y.D1;
92 obj.Dx = kron(Dx,Iy);
93 obj.Dy = kron(Ix,Dy);
94
95 % Boundary operators
96 obj.e_w = kr(ops_x.e_l, Iy);
97 obj.e_e = kr(ops_x.e_r, Iy);
98 obj.e_s = kr(Ix, ops_y.e_l);
99 obj.e_n = kr(Ix, ops_y.e_r);
100
101 obj.m = m;
102 obj.h = [ops_x.h ops_y.h];
103 obj.order = order;
104 obj.a = a;
105 obj.coupling_type = coupling_type;
106 obj.interpolation_type = interpolation_type;
107 obj.interpolation_damping = interpolation_damping;
108 obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy);
109
110 end
111 % Closure functions return the opertors applied to the own domain to close the boundary
112 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
113 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
114 % type is a string specifying the type of boundary condition if there are several.
115 % data is a function returning the data that should be applied at the boundary.
116 % neighbour_scheme is an instance of Scheme that should be interfaced to.
117 % neighbour_boundary is a string specifying which boundary to interface to.
118 function [closure, penalty] = boundary_condition(obj,boundary,type)
119 default_arg('type','dirichlet');
120
121 sigma = -1; % Scalar penalty parameter
122 switch boundary
123 case {'w','W','west','West'}
124 tau = sigma*obj.a{1}*obj.e_w*obj.H_y;
125 closure = obj.Hi*tau*obj.e_w';
126
127 case {'s','S','south','South'}
128 tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
129 closure = obj.Hi*tau*obj.e_s';
130 end
131 penalty = -obj.Hi*tau;
132
133 end
134
135 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
136
137 % Get neighbour boundary operator
138 switch neighbour_boundary
139 case {'e','E','east','East'}
140 e_neighbour = neighbour_scheme.e_e;
141 m_neighbour = neighbour_scheme.m(2);
142 case {'w','W','west','West'}
143 e_neighbour = neighbour_scheme.e_w;
144 m_neighbour = neighbour_scheme.m(2);
145 case {'n','N','north','North'}
146 e_neighbour = neighbour_scheme.e_n;
147 m_neighbour = neighbour_scheme.m(1);
148 case {'s','S','south','South'}
149 e_neighbour = neighbour_scheme.e_s;
150 m_neighbour = neighbour_scheme.m(1);
151 end
152
153 switch obj.coupling_type
154
155 % Upwind coupling (energy dissipation)
156 case 'upwind'
157 sigma_ds = -1; %"Downstream" penalty
158 sigma_us = 0; %"Upstream" penalty
159
160 % Energy-preserving coupling (no energy dissipation)
161 case 'centered'
162 sigma_ds = -1/2; %"Downstream" penalty
163 sigma_us = 1/2; %"Upstream" penalty
164
165 otherwise
166 error(['Interface coupling type ' coupling_type ' is not available.'])
167 end
168
169 % Check grid ratio for interpolation
170 switch boundary
171 case {'w','W','west','West','e','E','east','East'}
172 m = obj.m(2);
173 case {'s','S','south','South','n','N','north','North'}
174 m = obj.m(1);
175 end
176 grid_ratio = m/m_neighbour;
177 if grid_ratio ~= 1
178
179 [ms, index] = sort([m, m_neighbour]);
180 orders = [obj.order, neighbour_scheme.order];
181 orders = orders(index);
182
183 switch obj.interpolation_type
184 case 'MC'
185 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
186 if grid_ratio < 1
187 I_neighbour2local_us = interpOpSet.IF2C;
188 I_neighbour2local_ds = interpOpSet.IF2C;
189 I_local2neighbour_us = interpOpSet.IC2F;
190 I_local2neighbour_ds = interpOpSet.IC2F;
191 elseif grid_ratio > 1
192 I_neighbour2local_us = interpOpSet.IC2F;
193 I_neighbour2local_ds = interpOpSet.IC2F;
194 I_local2neighbour_us = interpOpSet.IF2C;
195 I_local2neighbour_ds = interpOpSet.IF2C;
196 end
197 case 'AWW'
198 %String 'C2F' indicates that ICF2 is more accurate.
199 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
200 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
201 if grid_ratio < 1
202 % Local is coarser than neighbour
203 I_neighbour2local_us = interpOpSetC2F.IF2C;
204 I_neighbour2local_ds = interpOpSetF2C.IF2C;
205 I_local2neighbour_us = interpOpSetC2F.IC2F;
206 I_local2neighbour_ds = interpOpSetF2C.IC2F;
207 elseif grid_ratio > 1
208 % Local is finer than neighbour
209 I_neighbour2local_us = interpOpSetF2C.IC2F;
210 I_neighbour2local_ds = interpOpSetC2F.IC2F;
211 I_local2neighbour_us = interpOpSetF2C.IF2C;
212 I_local2neighbour_ds = interpOpSetC2F.IF2C;
213 end
214 otherwise
215 error(['Interpolation type ' obj.interpolation_type ...
216 ' is not available.' ]);
217 end
218
219 else
220 % No interpolation required
221 I_neighbour2local_us = speye(m,m);
222 I_neighbour2local_ds = speye(m,m);
223 end
224
225 int_damp_us = obj.interpolation_damping{1};
226 int_damp_ds = obj.interpolation_damping{2};
227
228 I = speye(m,m);
229 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
230 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
231
232
233 switch boundary
234 case {'w','W','west','West'}
235 tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y;
236 closure = obj.Hi*tau*obj.e_w';
237 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
238
239 beta = int_damp_ds*obj.a{1}...
240 *obj.e_w*obj.H_y;
241 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w';
242 case {'e','E','east','East'}
243 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
244 closure = obj.Hi*tau*obj.e_e';
245 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
246
247 beta = int_damp_us*obj.a{1}...
248 *obj.e_e*obj.H_y;
249 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e';
250 case {'s','S','south','South'}
251 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
252 closure = obj.Hi*tau*obj.e_s';
253 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
254
255 beta = int_damp_ds*obj.a{2}...
256 *obj.e_s*obj.H_x;
257 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s';
258 case {'n','N','north','North'}
259 tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x;
260 closure = obj.Hi*tau*obj.e_n';
261 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
262
263 beta = int_damp_us*obj.a{2}...
264 *obj.e_n*obj.H_x;
265 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n';
266 end
267
268
269 end
270
271 function N = size(obj)
272 N = obj.m;
273 end
274
275 end
276
277 methods(Static)
278 % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u
279 % and bound_v of scheme schm_v.
280 % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l')
281 function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v)
282 [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v);
283 [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
284 end
285 end
286 end