Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 905:459eeb99130f feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Include type as (optional) input parameter in the interface method of all schemes.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 22 Nov 2018 22:03:44 -0800
parents f4595f14d696
children b9c98661ff5d
comparison
equal deleted inserted replaced
904:14b093a344eb 905:459eeb99130f
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 grid % Grid 5 grid % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
7 v0 % Initial data 7 v0 % Initial data
8 8
9 a % Wave speed a = [a1, a2]; 9 a % Wave speed a = [a1, a2];
10 % Can either be a constant vector or a cell array of function handles. 10 % Can either be a constant vector or a cell array of function handles.
11 11
12 H % Discrete norm 12 H % Discrete norm
13 H_x, H_y % Norms in the x and y directions 13 H_x, H_y % Norms in the x and y directions
14 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms 14 Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms
15 15
16 % Derivatives 16 % Derivatives
17 Dx, Dy 17 Dx, Dy
18 18
19 % Boundary operators 19 % Boundary operators
20 e_w, e_e, e_s, e_n 20 e_w, e_e, e_s, e_n
21 21
22 D % Total discrete operator 22 D % Total discrete operator
23 23
24 % String, type of interface coupling 24 % String, type of interface coupling
25 % Default: 'upwind' 25 % Default: 'upwind'
26 % Other: 'centered' 26 % Other: 'centered'
27 coupling_type 27 coupling_type
28 28
29 % String, type of interpolation operators 29 % String, type of interpolation operators
30 % Default: 'AWW' (Almquist Wang Werpers) 30 % Default: 'AWW' (Almquist Wang Werpers)
31 % Other: 'MC' (Mattsson Carpenter) 31 % Other: 'MC' (Mattsson Carpenter)
32 interpolation_type 32 interpolation_type
33 33
34 34
35 % Cell array, damping on upwstream and downstream sides. 35 % Cell array, damping on upwstream and downstream sides.
36 interpolation_damping 36 interpolation_damping
37 37
38 end 38 end
39 39
40 40
41 methods 41 methods
42 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping) 42 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping)
43 43
44 default_arg('interpolation_damping',{0,0}); 44 default_arg('interpolation_damping',{0,0});
45 default_arg('interpolation_type','AWW'); 45 default_arg('interpolation_type','AWW');
46 default_arg('coupling_type','upwind'); 46 default_arg('coupling_type','upwind');
47 default_arg('a',1/sqrt(2)*[1, 1]); 47 default_arg('a',1/sqrt(2)*[1, 1]);
48 default_arg('opSet',@sbp.D2Standard); 48 default_arg('opSet',@sbp.D2Standard);
49 49
50 assert(isa(g, 'grid.Cartesian')) 50 assert(isa(g, 'grid.Cartesian'))
51 if iscell(a) 51 if iscell(a)
52 a1 = grid.evalOn(g, a{1}); 52 a1 = grid.evalOn(g, a{1});
53 a2 = grid.evalOn(g, a{2}); 53 a2 = grid.evalOn(g, a{2});
54 a = {spdiag(a1), spdiag(a2)}; 54 a = {spdiag(a1), spdiag(a2)};
55 else 55 else
56 a = {a(1), a(2)}; 56 a = {a(1), a(2)};
57 end 57 end
58 58
59 m = g.size(); 59 m = g.size();
60 m_x = m(1); 60 m_x = m(1);
61 m_y = m(2); 61 m_y = m(2);
62 m_tot = g.N(); 62 m_tot = g.N();
63 63
68 % Operator sets 68 % Operator sets
69 ops_x = opSet(m_x, xlim, order); 69 ops_x = opSet(m_x, xlim, order);
70 ops_y = opSet(m_y, ylim, order); 70 ops_y = opSet(m_y, ylim, order);
71 Ix = speye(m_x); 71 Ix = speye(m_x);
72 Iy = speye(m_y); 72 Iy = speye(m_y);
73 73
74 % Norms 74 % Norms
75 Hx = ops_x.H; 75 Hx = ops_x.H;
76 Hy = ops_y.H; 76 Hy = ops_y.H;
77 Hxi = ops_x.HI; 77 Hxi = ops_x.HI;
78 Hyi = ops_y.HI; 78 Hyi = ops_y.HI;
79 79
80 obj.H_x = Hx; 80 obj.H_x = Hx;
81 obj.H_y = Hy; 81 obj.H_y = Hy;
82 obj.H = kron(Hx,Hy); 82 obj.H = kron(Hx,Hy);
83 obj.Hi = kron(Hxi,Hyi); 83 obj.Hi = kron(Hxi,Hyi);
84 obj.Hx = kron(Hx,Iy); 84 obj.Hx = kron(Hx,Iy);
85 obj.Hy = kron(Ix,Hy); 85 obj.Hy = kron(Ix,Hy);
86 obj.Hxi = kron(Hxi,Iy); 86 obj.Hxi = kron(Hxi,Iy);
87 obj.Hyi = kron(Ix,Hyi); 87 obj.Hyi = kron(Ix,Hyi);
88 88
89 % Derivatives 89 % Derivatives
90 Dx = ops_x.D1; 90 Dx = ops_x.D1;
91 Dy = ops_y.D1; 91 Dy = ops_y.D1;
92 obj.Dx = kron(Dx,Iy); 92 obj.Dx = kron(Dx,Iy);
93 obj.Dy = kron(Ix,Dy); 93 obj.Dy = kron(Ix,Dy);
94 94
95 % Boundary operators 95 % Boundary operators
96 obj.e_w = kr(ops_x.e_l, Iy); 96 obj.e_w = kr(ops_x.e_l, Iy);
97 obj.e_e = kr(ops_x.e_r, Iy); 97 obj.e_e = kr(ops_x.e_r, Iy);
98 obj.e_s = kr(Ix, ops_y.e_l); 98 obj.e_s = kr(Ix, ops_y.e_l);
99 obj.e_n = kr(Ix, ops_y.e_r); 99 obj.e_n = kr(Ix, ops_y.e_r);
115 % data is a function returning the data that should be applied at the boundary. 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. 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. 117 % neighbour_boundary is a string specifying which boundary to interface to.
118 function [closure, penalty] = boundary_condition(obj,boundary,type) 118 function [closure, penalty] = boundary_condition(obj,boundary,type)
119 default_arg('type','dirichlet'); 119 default_arg('type','dirichlet');
120 120
121 sigma = -1; % Scalar penalty parameter 121 sigma = -1; % Scalar penalty parameter
122 switch boundary 122 switch boundary
123 case {'w','W','west','West'} 123 case {'w','W','west','West'}
124 tau = sigma*obj.a{1}*obj.e_w*obj.H_y; 124 tau = sigma*obj.a{1}*obj.e_w*obj.H_y;
125 closure = obj.Hi*tau*obj.e_w'; 125 closure = obj.Hi*tau*obj.e_w';
126 126
127 case {'s','S','south','South'} 127 case {'s','S','south','South'}
128 tau = sigma*obj.a{2}*obj.e_s*obj.H_x; 128 tau = sigma*obj.a{2}*obj.e_s*obj.H_x;
129 closure = obj.Hi*tau*obj.e_s'; 129 closure = obj.Hi*tau*obj.e_s';
130 end 130 end
131 penalty = -obj.Hi*tau; 131 penalty = -obj.Hi*tau;
132 132
133 end 133 end
134 134
135 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 135 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
136 136
137 % Get neighbour boundary operator 137 % Get neighbour boundary operator
138 switch neighbour_boundary 138 switch neighbour_boundary
139 case {'e','E','east','East'} 139 case {'e','E','east','East'}
140 e_neighbour = neighbour_scheme.e_e; 140 e_neighbour = neighbour_scheme.e_e;
141 m_neighbour = neighbour_scheme.m(2); 141 m_neighbour = neighbour_scheme.m(2);
147 m_neighbour = neighbour_scheme.m(1); 147 m_neighbour = neighbour_scheme.m(1);
148 case {'s','S','south','South'} 148 case {'s','S','south','South'}
149 e_neighbour = neighbour_scheme.e_s; 149 e_neighbour = neighbour_scheme.e_s;
150 m_neighbour = neighbour_scheme.m(1); 150 m_neighbour = neighbour_scheme.m(1);
151 end 151 end
152 152
153 switch obj.coupling_type 153 switch obj.coupling_type
154 154
155 % Upwind coupling (energy dissipation) 155 % Upwind coupling (energy dissipation)
156 case 'upwind' 156 case 'upwind'
157 sigma_ds = -1; %"Downstream" penalty 157 sigma_ds = -1; %"Downstream" penalty
158 sigma_us = 0; %"Upstream" penalty 158 sigma_us = 0; %"Upstream" penalty
159 159
167 end 167 end
168 168
169 % Check grid ratio for interpolation 169 % Check grid ratio for interpolation
170 switch boundary 170 switch boundary
171 case {'w','W','west','West','e','E','east','East'} 171 case {'w','W','west','West','e','E','east','East'}
172 m = obj.m(2); 172 m = obj.m(2);
173 case {'s','S','south','South','n','N','north','North'} 173 case {'s','S','south','South','n','N','north','North'}
174 m = obj.m(1); 174 m = obj.m(1);
175 end 175 end
176 grid_ratio = m/m_neighbour; 176 grid_ratio = m/m_neighbour;
177 if grid_ratio ~= 1 177 if grid_ratio ~= 1
195 I_local2neighbour_ds = interpOpSet.IF2C; 195 I_local2neighbour_ds = interpOpSet.IF2C;
196 end 196 end
197 case 'AWW' 197 case 'AWW'
198 %String 'C2F' indicates that ICF2 is more accurate. 198 %String 'C2F' indicates that ICF2 is more accurate.
199 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); 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'); 200 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
201 if grid_ratio < 1 201 if grid_ratio < 1
202 % Local is coarser than neighbour 202 % Local is coarser than neighbour
203 I_neighbour2local_us = interpOpSetC2F.IF2C; 203 I_neighbour2local_us = interpOpSetC2F.IF2C;
204 I_neighbour2local_ds = interpOpSetF2C.IF2C; 204 I_neighbour2local_ds = interpOpSetF2C.IF2C;
205 I_local2neighbour_us = interpOpSetC2F.IC2F; 205 I_local2neighbour_us = interpOpSetC2F.IC2F;
206 I_local2neighbour_ds = interpOpSetF2C.IC2F; 206 I_local2neighbour_ds = interpOpSetF2C.IC2F;
207 elseif grid_ratio > 1 207 elseif grid_ratio > 1
208 % Local is finer than neighbour 208 % Local is finer than neighbour
209 I_neighbour2local_us = interpOpSetF2C.IC2F; 209 I_neighbour2local_us = interpOpSetF2C.IC2F;
210 I_neighbour2local_ds = interpOpSetC2F.IC2F; 210 I_neighbour2local_ds = interpOpSetC2F.IC2F;
211 I_local2neighbour_us = interpOpSetF2C.IF2C; 211 I_local2neighbour_us = interpOpSetF2C.IF2C;
212 I_local2neighbour_ds = interpOpSetC2F.IF2C; 212 I_local2neighbour_ds = interpOpSetC2F.IF2C;
213 end 213 end
214 otherwise 214 otherwise
215 error(['Interpolation type ' obj.interpolation_type ... 215 error(['Interpolation type ' obj.interpolation_type ...
216 ' is not available.' ]); 216 ' is not available.' ]);
217 end 217 end
218 218
219 else 219 else
220 % No interpolation required 220 % No interpolation required
221 I_neighbour2local_us = speye(m,m); 221 I_neighbour2local_us = speye(m,m);
222 I_neighbour2local_ds = speye(m,m); 222 I_neighbour2local_ds = speye(m,m);
223 end 223 end
224 224
225 int_damp_us = obj.interpolation_damping{1}; 225 int_damp_us = obj.interpolation_damping{1};
226 int_damp_ds = obj.interpolation_damping{2}; 226 int_damp_ds = obj.interpolation_damping{2};
227 227
228 I = speye(m,m); 228 I = speye(m,m);
229 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us; 229 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
236 closure = obj.Hi*tau*obj.e_w'; 236 closure = obj.Hi*tau*obj.e_w';
237 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 237 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
238 238
239 beta = int_damp_ds*obj.a{1}... 239 beta = int_damp_ds*obj.a{1}...
240 *obj.e_w*obj.H_y; 240 *obj.e_w*obj.H_y;
241 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; 241 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w';
242 case {'e','E','east','East'} 242 case {'e','E','east','East'}
243 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; 243 tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y;
244 closure = obj.Hi*tau*obj.e_e'; 244 closure = obj.Hi*tau*obj.e_e';
245 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 245 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
246 246
247 beta = int_damp_us*obj.a{1}... 247 beta = int_damp_us*obj.a{1}...
248 *obj.e_e*obj.H_y; 248 *obj.e_e*obj.H_y;
249 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; 249 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e';
250 case {'s','S','south','South'} 250 case {'s','S','south','South'}
251 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; 251 tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x;
252 closure = obj.Hi*tau*obj.e_s'; 252 closure = obj.Hi*tau*obj.e_s';
253 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 253 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
254 254
255 beta = int_damp_ds*obj.a{2}... 255 beta = int_damp_ds*obj.a{2}...
256 *obj.e_s*obj.H_x; 256 *obj.e_s*obj.H_x;
257 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; 257 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s';
260 closure = obj.Hi*tau*obj.e_n'; 260 closure = obj.Hi*tau*obj.e_n';
261 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 261 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
262 262
263 beta = int_damp_us*obj.a{2}... 263 beta = int_damp_us*obj.a{2}...
264 *obj.e_n*obj.H_x; 264 *obj.e_n*obj.H_x;
265 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; 265 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n';
266 end 266 end
267 267
268 268
269 end 269 end
270 270
271 function N = size(obj) 271 function N = size(obj)
272 N = obj.m; 272 N = obj.m;
273 end 273 end
274 274
275 end 275 end