Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 912:1cdf5ead2a16 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor interface method for Schrodinger2d
author Martin Almquist <malmquist@stanford.edu>
date Sun, 25 Nov 2018 19:46:39 -0800
parents b9c98661ff5d
children 46f5dc61d90b
comparison
equal deleted inserted replaced
911:f7306f03f77a 912:1cdf5ead2a16
31 e_w, e_e, e_s, e_n 31 e_w, e_e, e_s, e_n
32 d_w, d_e, d_s, d_n 32 d_w, d_e, d_s, d_n
33 33
34 H_boundary % Boundary inner products 34 H_boundary % Boundary inner products
35 35
36 interpolation_type % MC or AWW
37
38 end 36 end
39 37
40 methods 38 methods
41 39
42 function obj = Schrodinger2d(g ,order, a, opSet, interpolation_type) 40 function obj = Schrodinger2d(g ,order, a, opSet)
43 default_arg('interpolation_type','AWW');
44 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); 41 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
45 default_arg('a',1); 42 default_arg('a',1);
46 dim = 2; 43 dim = 2;
47 44
48 assert(isa(g, 'grid.Cartesian')) 45 assert(isa(g, 'grid.Cartesian'))
148 obj.e_n = obj.e_r{2}; 145 obj.e_n = obj.e_r{2};
149 obj.d_w = obj.d1_l{1}; 146 obj.d_w = obj.d1_l{1};
150 obj.d_e = obj.d1_r{1}; 147 obj.d_e = obj.d1_r{1};
151 obj.d_s = obj.d1_l{2}; 148 obj.d_s = obj.d1_l{2};
152 obj.d_n = obj.d1_r{2}; 149 obj.d_n = obj.d1_r{2};
153 obj.interpolation_type = interpolation_type;
154 150
155 end 151 end
156 152
157 153
158 % Closure functions return the operators applied to the own domain to close the boundary 154 % Closure functions return the operators applied to the own domain to close the boundary
200 otherwise 196 otherwise
201 error('No such boundary condition: type = %s',type); 197 error('No such boundary condition: type = %s',type);
202 end 198 end
203 end 199 end
204 200
201 % opts Struct that specifies the interface coupling.
202 % Fields:
203 % -- interpolation: struct of interpolation operators (empty for conforming grids)
205 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts) 204 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
205 if isempty(opts)
206 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary);
207 else
208 assertType(opts, 'struct');
209 if isfield(opts, 'I_local2neighbor')
210 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
211 else
212 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts);
213 end
214 end
215 end
216
217 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
206 % u denotes the solution in the own domain 218 % u denotes the solution in the own domain
207 % v denotes the solution in the neighbour domain 219 % v denotes the solution in the neighbour domain
208 % Get neighbour boundary operator 220
209 221 % Get boundary operators
210 [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary); 222 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
211 [coord, n] = get_boundary_number(obj, boundary); 223 [e, d, H_gamma] = obj.get_boundary_ops(boundary);
212 switch n_nei
213 case 1
214 % North or east boundary
215 e_neighbour = neighbour_scheme.e_r;
216 d_neighbour = neighbour_scheme.d1_r;
217 case -1
218 % South or west boundary
219 e_neighbour = neighbour_scheme.e_l;
220 d_neighbour = neighbour_scheme.d1_l;
221 end
222
223 e_neighbour = e_neighbour{coord_nei};
224 d_neighbour = d_neighbour{coord_nei};
225 H_gamma = obj.H_boundary{coord};
226 Hi = obj.Hi; 224 Hi = obj.Hi;
227 a = obj.a; 225 a = obj.a;
228 226
229 switch coord_nei 227 % Get outward unit normal component
230 case 1 228 [~, n] = obj.get_boundary_number(boundary);
231 m_neighbour = neighbour_scheme.m(2);
232 case 2
233 m_neighbour = neighbour_scheme.m(1);
234 end
235
236 switch coord
237 case 1
238 m = obj.m(2);
239 case 2
240 m = obj.m(1);
241 end
242
243 switch n
244 case 1
245 % North or east boundary
246 e = obj.e_r;
247 d = obj.d1_r;
248 case -1
249 % South or west boundary
250 e = obj.e_l;
251 d = obj.d1_l;
252 end
253 e = e{coord};
254 d = d{coord};
255 229
256 Hi = obj.Hi; 230 Hi = obj.Hi;
257 sigma = -n*1i*a/2; 231 sigma = -n*1i*a/2;
258 tau = -n*(1i*a)'/2; 232 tau = -n*(1i*a)'/2;
259 233
260 grid_ratio = m/m_neighbour;
261 if grid_ratio ~= 1
262
263 [ms, index] = sort([m, m_neighbour]);
264 orders = [obj.order, neighbour_scheme.order];
265 orders = orders(index);
266
267 switch obj.interpolation_type
268 case 'MC'
269 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
270 if grid_ratio < 1
271 I_neighbour2local_e = interpOpSet.IF2C;
272 I_neighbour2local_d = interpOpSet.IF2C;
273 I_local2neighbour_e = interpOpSet.IC2F;
274 I_local2neighbour_d = interpOpSet.IC2F;
275 elseif grid_ratio > 1
276 I_neighbour2local_e = interpOpSet.IC2F;
277 I_neighbour2local_d = interpOpSet.IC2F;
278 I_local2neighbour_e = interpOpSet.IF2C;
279 I_local2neighbour_d = interpOpSet.IF2C;
280 end
281 case 'AWW'
282 %String 'C2F' indicates that ICF2 is more accurate.
283 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
284 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
285 if grid_ratio < 1
286 % Local is coarser than neighbour
287 I_neighbour2local_e = interpOpSetF2C.IF2C;
288 I_neighbour2local_d = interpOpSetC2F.IF2C;
289 I_local2neighbour_e = interpOpSetC2F.IC2F;
290 I_local2neighbour_d = interpOpSetF2C.IC2F;
291 elseif grid_ratio > 1
292 % Local is finer than neighbour
293 I_neighbour2local_e = interpOpSetC2F.IC2F;
294 I_neighbour2local_d = interpOpSetF2C.IC2F;
295 I_local2neighbour_e = interpOpSetF2C.IF2C;
296 I_local2neighbour_d = interpOpSetC2F.IF2C;
297 end
298 otherwise
299 error(['Interpolation type ' obj.interpolation_type ...
300 ' is not available.' ]);
301 end
302
303 else
304 % No interpolation required
305 I_neighbour2local_e = speye(m,m);
306 I_neighbour2local_d = speye(m,m);
307 I_local2neighbour_e = speye(m,m);
308 I_local2neighbour_d = speye(m,m);
309 end
310
311 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; 234 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d';
312 penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ... 235 penalty = -tau*Hi*d*H_gamma*e_neighbour' ...
313 -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; 236 -sigma*Hi*e*H_gamma*d_neighbour';
237
238 end
239
240 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,opts)
241 % u denotes the solution in the own domain
242 % v denotes the solution in the neighbour domain
243
244 % Get boundary operators
245 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
246 [e, d, H_gamma] = obj.get_boundary_ops(boundary);
247 Hi = obj.Hi;
248 a = obj.a;
249
250 % Get outward unit normal component
251 [~, n] = obj.get_boundary_number(boundary);
252
253 % Get interpolation operators from opts
254 I_u2v_good = opts.I_local2neighbor.good;
255 I_u2v_bad = opts.I_local2neighbor.bad;
256 I_v2u_good = opts.I_neighbor2local.good;
257 I_v2u_bad = opts.I_neighbor2local.bad;
258
259 sigma = -n*1i*a/2;
260 tau = -n*(1i*a)'/2;
261
262 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d';
263 penalty = -tau*Hi*d*H_gamma*I_v2u_good*e_neighbour' ...
264 -sigma*Hi*e*H_gamma*I_v2u_bad*d_neighbour';
314 265
315 end 266 end
316 267
317 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 268 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
318 function [j, nj] = get_boundary_number(obj, boundary) 269 function [j, nj] = get_boundary_number(obj, boundary)
332 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 283 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
333 nj = 1; 284 nj = 1;
334 end 285 end
335 end 286 end
336 287
337 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 288 % Returns the boundary ops and sign for the boundary specified by the string boundary.
338 function [return_op] = get_boundary_operator(obj, op, boundary) 289 % The right boundary is considered the positive boundary
290 function [e, d, H_b] = get_boundary_ops(obj, boundary)
339 291
340 switch boundary 292 switch boundary
341 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 293 case 'w'
342 j = 1; 294 e = obj.e_w;
343 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 295 d = obj.d_w;
344 j = 2; 296 H_b = obj.H_boundary{1};
297 case 'e'
298 e = obj.e_e;
299 d = obj.d_e;
300 H_b = obj.H_boundary{1};
301 case 's'
302 e = obj.e_s;
303 d = obj.d_s;
304 H_b = obj.H_boundary{2};
305 case 'n'
306 e = obj.e_n;
307 d = obj.d_n;
308 H_b = obj.H_boundary{2};
345 otherwise 309 otherwise
346 error('No such boundary: boundary = %s',boundary); 310 error('No such boundary: boundary = %s',boundary);
347 end 311 end
348
349 switch op
350 case 'e'
351 switch boundary
352 case {'w','W','west','West','s','S','south','South'}
353 return_op = obj.e_l{j};
354 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
355 return_op = obj.e_r{j};
356 end
357 case 'd'
358 switch boundary
359 case {'w','W','west','West','s','S','south','South'}
360 return_op = obj.d1_l{j};
361 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
362 return_op = obj.d1_r{j};
363 end
364 otherwise
365 error(['No such operator: operator = ' op]);
366 end
367
368 end 312 end
369 313
370 function N = size(obj) 314 function N = size(obj)
371 N = prod(obj.m); 315 N = prod(obj.m);
372 end 316 end