Mercurial > repos > public > sbplib

comparison +scheme/Utux2D.m @ 666:2d85f17a8aec feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add possibility for damping terms at interpolation interface.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 16 Oct 2017 21:56:12 -0700
parents b7b3c11fab4d
children f4595f14d696
comparison
equal deleted inserted replaced
610:b7b3c11fab4d 666:2d85f17a8aec
29 % Default: 'AWW' (Almquist Wang Werpers) 29 % Default: 'AWW' (Almquist Wang Werpers)
30 % Other: 'MC' (Mattsson Carpenter) 30 % Other: 'MC' (Mattsson Carpenter)
31 interpolation_type 31 interpolation_type
32 32
33 33
34 % Cell array, damping on upwstream and downstream sides.
35 interpolation_damping
36
34 end 37 end
35 38
36 39
37 methods 40 methods
38 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type) 41 function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping)
39 42
43 default_arg('interpolation_damping',{0,0});
40 default_arg('interpolation_type','AWW'); 44 default_arg('interpolation_type','AWW');
41 default_arg('coupling_type','upwind'); 45 default_arg('coupling_type','upwind');
42 default_arg('a',1/sqrt(2)*[1, 1]); 46 default_arg('a',1/sqrt(2)*[1, 1]);
43 default_arg('opSet',@sbp.D2Standard); 47 default_arg('opSet',@sbp.D2Standard);
44 assert(isa(g, 'grid.Cartesian')) 48 assert(isa(g, 'grid.Cartesian'))
89 obj.h = [ops_x.h ops_y.h]; 93 obj.h = [ops_x.h ops_y.h];
90 obj.order = order; 94 obj.order = order;
91 obj.a = a; 95 obj.a = a;
92 obj.coupling_type = coupling_type; 96 obj.coupling_type = coupling_type;
93 obj.interpolation_type = interpolation_type; 97 obj.interpolation_type = interpolation_type;
98 obj.interpolation_damping = interpolation_damping;
94 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); 99 obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy);
95 100
96 end 101 end
97 % Closure functions return the opertors applied to the own domain to close the boundary 102 % Closure functions return the opertors applied to the own domain to close the boundary
98 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 103 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
170 case 'MC' 175 case 'MC'
171 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); 176 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
172 if grid_ratio < 1 177 if grid_ratio < 1
173 I_neighbour2local_us = interpOpSet.IF2C; 178 I_neighbour2local_us = interpOpSet.IF2C;
174 I_neighbour2local_ds = interpOpSet.IF2C; 179 I_neighbour2local_ds = interpOpSet.IF2C;
180 I_local2neighbour_us = interpOpSet.IC2F;
181 I_local2neighbour_ds = interpOpSet.IC2F;
175 elseif grid_ratio > 1 182 elseif grid_ratio > 1
176 I_neighbour2local_us = interpOpSet.IC2F; 183 I_neighbour2local_us = interpOpSet.IC2F;
177 I_neighbour2local_ds = interpOpSet.IC2F; 184 I_neighbour2local_ds = interpOpSet.IC2F;
185 I_local2neighbour_us = interpOpSet.IF2C;
186 I_local2neighbour_ds = interpOpSet.IF2C;
178 end 187 end
179 case 'AWW' 188 case 'AWW'
180 %String 'C2F' indicates that ICF2 is more accurate. 189 %String 'C2F' indicates that ICF2 is more accurate.
181 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); 190 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
182 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); 191 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
183 if grid_ratio < 1 192 if grid_ratio < 1
184 % Local is coarser than neighbour 193 % Local is coarser than neighbour
185 I_neighbour2local_us = interpOpSetC2F.IF2C; 194 I_neighbour2local_us = interpOpSetC2F.IF2C;
186 I_neighbour2local_ds = interpOpSetF2C.IF2C; 195 I_neighbour2local_ds = interpOpSetF2C.IF2C;
196 I_local2neighbour_us = interpOpSetC2F.IC2F;
197 I_local2neighbour_ds = interpOpSetF2C.IC2F;
187 elseif grid_ratio > 1 198 elseif grid_ratio > 1
188 % Local is finer than neighbour 199 % Local is finer than neighbour
189 I_neighbour2local_us = interpOpSetF2C.IC2F; 200 I_neighbour2local_us = interpOpSetF2C.IC2F;
190 I_neighbour2local_ds = interpOpSetC2F.IC2F; 201 I_neighbour2local_ds = interpOpSetC2F.IC2F;
202 I_local2neighbour_us = interpOpSetF2C.IF2C;
203 I_local2neighbour_ds = interpOpSetC2F.IF2C;
191 end 204 end
192 otherwise 205 otherwise
193 error(['Interpolation type ' obj.interpolation_type ... 206 error(['Interpolation type ' obj.interpolation_type ...
194 ' is not available.' ]); 207 ' is not available.' ]);
195 end 208 end
198 % No interpolation required 211 % No interpolation required
199 I_neighbour2local_us = speye(m,m); 212 I_neighbour2local_us = speye(m,m);
200 I_neighbour2local_ds = speye(m,m); 213 I_neighbour2local_ds = speye(m,m);
201 end 214 end
202 215
216 int_damp_us = obj.interpolation_damping{1};
217 int_damp_ds = obj.interpolation_damping{2};
218
219 I = speye(m,m);
220 I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us;
221 I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds;
222
223
203 switch boundary 224 switch boundary
204 case {'w','W','west','West'} 225 case {'w','W','west','West'}
205 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; 226 tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y;
206 closure = obj.Hi*tau*obj.e_w'; 227 closure = obj.Hi*tau*obj.e_w';
207 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 228 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
229
230 beta = int_damp_ds*obj.a(1)...
231 *obj.e_w*obj.H_y;
232 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w';
208 case {'e','E','east','East'} 233 case {'e','E','east','East'}
209 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; 234 tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y;
210 closure = obj.Hi*tau*obj.e_e'; 235 closure = obj.Hi*tau*obj.e_e';
211 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 236 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
237
238 beta = int_damp_us*obj.a(1)...
239 *obj.e_e*obj.H_y;
240 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e';
212 case {'s','S','south','South'} 241 case {'s','S','south','South'}
213 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; 242 tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x;
214 closure = obj.Hi*tau*obj.e_s'; 243 closure = obj.Hi*tau*obj.e_s';
215 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; 244 penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour';
245
246 beta = int_damp_ds*obj.a(2)...
247 *obj.e_s*obj.H_x;
248 closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s';
216 case {'n','N','north','North'} 249 case {'n','N','north','North'}
217 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; 250 tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x;
218 closure = obj.Hi*tau*obj.e_n'; 251 closure = obj.Hi*tau*obj.e_n';
219 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; 252 penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour';
253
254 beta = int_damp_us*obj.a(2)...
255 *obj.e_n*obj.H_x;
256 closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n';
220 end 257 end
221 258
222 259
223 end 260 end
224 261