Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.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
1 classdef Schrodinger2d < scheme.Scheme 1 classdef Schrodinger2d < scheme.Scheme
2 2
3 % Discretizes the Laplacian with constant coefficent, 3 % Discretizes the Laplacian with constant coefficent,
4 % in the Schrödinger equation way (i.e., the discretization matrix is not necessarily 4 % in the Schrödinger equation way (i.e., the discretization matrix is not necessarily
5 % definite) 5 % definite)
6 % u_t = a*i*Laplace u 6 % u_t = a*i*Laplace u
7 % opSet should be cell array of opSets, one per dimension. This 7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction. 8 % is useful if we have periodic BC in one direction.
9 9
10 properties 10 properties
11 m % Number of points in each direction, possibly a vector 11 m % Number of points in each direction, possibly a vector
28 H, Hi % Inner products 28 H, Hi % Inner products
29 e_l, e_r 29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary 30 d1_l, d1_r % Normal derivatives at the boundary
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 36 interpolation_type % MC or AWW
37 37
38 end 38 end
165 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 165 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
166 default_arg('type','Neumann'); 166 default_arg('type','Neumann');
167 default_arg('parameter', []); 167 default_arg('parameter', []);
168 168
169 % j is the coordinate direction of the boundary 169 % j is the coordinate direction of the boundary
170 % nj: outward unit normal component. 170 % nj: outward unit normal component.
171 % nj = -1 for west, south, bottom boundaries 171 % nj = -1 for west, south, bottom boundaries
172 % nj = 1 for east, north, top boundaries 172 % nj = 1 for east, north, top boundaries
173 [j, nj] = obj.get_boundary_number(boundary); 173 [j, nj] = obj.get_boundary_number(boundary);
174 switch nj 174 switch nj
175 case 1 175 case 1
186 186
187 switch type 187 switch type
188 188
189 % Dirichlet boundary condition 189 % Dirichlet boundary condition
190 case {'D','d','dirichlet','Dirichlet'} 190 case {'D','d','dirichlet','Dirichlet'}
191 closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); 191 closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' );
192 penalty = -nj*Hi*d{j}*a*1i*H_gamma; 192 penalty = -nj*Hi*d{j}*a*1i*H_gamma;
193 193
194 % Free boundary condition 194 % Free boundary condition
195 case {'N','n','neumann','Neumann'} 195 case {'N','n','neumann','Neumann'}
196 closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); 196 closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' );
197 penalty = nj*Hi*e{j}*a*1i*H_gamma; 197 penalty = nj*Hi*e{j}*a*1i*H_gamma;
198 198
199 % Unknown boundary condition 199 % Unknown boundary condition
200 otherwise 200 otherwise
201 error('No such boundary condition: type = %s',type); 201 error('No such boundary condition: type = %s',type);
202 end 202 end
203 end 203 end
204 204
205 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 205 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
206 % u denotes the solution in the own domain 206 % u denotes the solution in the own domain
207 % v denotes the solution in the neighbour domain 207 % v denotes the solution in the neighbour domain
208 % Get neighbour boundary operator 208 % Get neighbour boundary operator
209 209
210 [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary); 210 [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary);
249 % South or west boundary 249 % South or west boundary
250 e = obj.e_l; 250 e = obj.e_l;
251 d = obj.d1_l; 251 d = obj.d1_l;
252 end 252 end
253 e = e{coord}; 253 e = e{coord};
254 d = d{coord}; 254 d = d{coord};
255 255
256 Hi = obj.Hi; 256 Hi = obj.Hi;
257 sigma = -n*1i*a/2; 257 sigma = -n*1i*a/2;
258 tau = -n*(1i*a)'/2; 258 tau = -n*(1i*a)'/2;
259 259
279 I_local2neighbour_d = interpOpSet.IF2C; 279 I_local2neighbour_d = interpOpSet.IF2C;
280 end 280 end
281 case 'AWW' 281 case 'AWW'
282 %String 'C2F' indicates that ICF2 is more accurate. 282 %String 'C2F' indicates that ICF2 is more accurate.
283 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); 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'); 284 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
285 if grid_ratio < 1 285 if grid_ratio < 1
286 % Local is coarser than neighbour 286 % Local is coarser than neighbour
287 I_neighbour2local_e = interpOpSetF2C.IF2C; 287 I_neighbour2local_e = interpOpSetF2C.IF2C;
288 I_neighbour2local_d = interpOpSetC2F.IF2C; 288 I_neighbour2local_d = interpOpSetC2F.IF2C;
289 I_local2neighbour_e = interpOpSetC2F.IC2F; 289 I_local2neighbour_e = interpOpSetC2F.IC2F;
290 I_local2neighbour_d = interpOpSetF2C.IC2F; 290 I_local2neighbour_d = interpOpSetF2C.IC2F;
291 elseif grid_ratio > 1 291 elseif grid_ratio > 1
292 % Local is finer than neighbour 292 % Local is finer than neighbour
293 I_neighbour2local_e = interpOpSetC2F.IC2F; 293 I_neighbour2local_e = interpOpSetC2F.IC2F;
294 I_neighbour2local_d = interpOpSetF2C.IC2F; 294 I_neighbour2local_d = interpOpSetF2C.IC2F;
295 I_local2neighbour_e = interpOpSetF2C.IF2C; 295 I_local2neighbour_e = interpOpSetF2C.IF2C;
296 I_local2neighbour_d = interpOpSetC2F.IF2C; 296 I_local2neighbour_d = interpOpSetC2F.IF2C;
297 end 297 end
298 otherwise 298 otherwise
299 error(['Interpolation type ' obj.interpolation_type ... 299 error(['Interpolation type ' obj.interpolation_type ...
300 ' is not available.' ]); 300 ' is not available.' ]);
301 end 301 end
302 302
303 else 303 else
304 % No interpolation required 304 % No interpolation required
305 I_neighbour2local_e = speye(m,m); 305 I_neighbour2local_e = speye(m,m);
306 I_neighbour2local_d = speye(m,m); 306 I_neighbour2local_d = speye(m,m);
307 I_local2neighbour_e = speye(m,m); 307 I_local2neighbour_e = speye(m,m);
308 I_local2neighbour_d = speye(m,m); 308 I_local2neighbour_d = speye(m,m);
309 end 309 end
310 310
311 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; 311 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' ... 312 penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ...
313 -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; 313 -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour';
314 314
315 end 315 end
316 316
317 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 317 % 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) 318 function [j, nj] = get_boundary_number(obj, boundary)
319 319