Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 719:b3f8fb9cefd2 feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add interpolation to Schrödinger 2D scheme.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 11 Mar 2018 21:39:49 -0700
parents 71aa5828cbbf
children f4595f14d696
comparison
equal deleted inserted replaced
718:71aa5828cbbf 719:b3f8fb9cefd2
26 D2 26 D2
27 27
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
32 d_w, d_e, d_s, d_n
31 33
32 H_boundary % Boundary inner products 34 H_boundary % Boundary inner products
33 35
36 interpolation_type % MC or AWW
37
34 end 38 end
35 39
36 methods 40 methods
37 41
38 function obj = Schrodinger2d(g ,order, a, opSet) 42 function obj = Schrodinger2d(g ,order, a, opSet, interpolation_type)
43 default_arg('interpolation_type','AWW');
39 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); 44 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
40 default_arg('a',1); 45 default_arg('a',1);
41 dim = 2; 46 dim = 2;
42 47
43 assert(isa(g, 'grid.Cartesian')) 48 assert(isa(g, 'grid.Cartesian'))
131 obj.h = h; 136 obj.h = h;
132 obj.order = order; 137 obj.order = order;
133 obj.grid = g; 138 obj.grid = g;
134 obj.dim = dim; 139 obj.dim = dim;
135 obj.a = a; 140 obj.a = a;
141 obj.e_w = obj.e_l{1};
142 obj.e_e = obj.e_r{1};
143 obj.e_s = obj.e_l{2};
144 obj.e_n = obj.e_r{2};
145 obj.d_w = obj.d1_l{1};
146 obj.d_e = obj.d1_r{1};
147 obj.d_s = obj.d1_l{2};
148 obj.d_n = obj.d1_r{2};
149 obj.interpolation_type = interpolation_type;
136 150
137 end 151 end
138 152
139 153
140 % 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
185 end 199 end
186 200
187 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 201 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
188 % u denotes the solution in the own domain 202 % u denotes the solution in the own domain
189 % v denotes the solution in the neighbour domain 203 % v denotes the solution in the neighbour domain
190 error('Interface not implemented'); 204 % Get neighbour boundary operator
205
206 [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary);
207 [coord, n] = get_boundary_number(obj, boundary);
208 switch n_nei
209 case 1
210 % North or east boundary
211 e_neighbour = neighbour_scheme.e_r;
212 d_neighbour = neighbour_scheme.d1_r;
213 case -1
214 % South or west boundary
215 e_neighbour = neighbour_scheme.e_l;
216 d_neighbour = neighbour_scheme.d1_l;
217 end
218
219 e_neighbour = e_neighbour{coord_nei};
220 d_neighbour = d_neighbour{coord_nei};
221 H_gamma = obj.H_boundary{coord};
222 Hi = obj.Hi;
223 a = obj.a;
224
225 switch coord_nei
226 case 1
227 m_neighbour = neighbour_scheme.m(2);
228 case 2
229 m_neighbour = neighbour_scheme.m(1);
230 end
231
232 switch coord
233 case 1
234 m = obj.m(2);
235 case 2
236 m = obj.m(1);
237 end
238
239 switch n
240 case 1
241 % North or east boundary
242 e = obj.e_r;
243 d = obj.d1_r;
244 case -1
245 % South or west boundary
246 e = obj.e_l;
247 d = obj.d1_l;
248 end
249 e = e{coord};
250 d = d{coord};
251
252 Hi = obj.Hi;
253 sigma = -n*1i*a/2;
254 tau = -n*(1i*a)'/2;
255
256 grid_ratio = m/m_neighbour;
257 if grid_ratio ~= 1
258
259 [ms, index] = sort([m, m_neighbour]);
260 orders = [obj.order, neighbour_scheme.order];
261 orders = orders(index);
262
263 switch obj.interpolation_type
264 case 'MC'
265 interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2));
266 if grid_ratio < 1
267 I_neighbour2local_e = interpOpSet.IF2C;
268 I_neighbour2local_d = interpOpSet.IF2C;
269 I_local2neighbour_e = interpOpSet.IC2F;
270 I_local2neighbour_d = interpOpSet.IC2F;
271 elseif grid_ratio > 1
272 I_neighbour2local_e = interpOpSet.IC2F;
273 I_neighbour2local_d = interpOpSet.IC2F;
274 I_local2neighbour_e = interpOpSet.IF2C;
275 I_local2neighbour_d = interpOpSet.IF2C;
276 end
277 case 'AWW'
278 %String 'C2F' indicates that ICF2 is more accurate.
279 interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C');
280 interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F');
281 if grid_ratio < 1
282 % Local is coarser than neighbour
283 I_neighbour2local_e = interpOpSetF2C.IF2C;
284 I_neighbour2local_d = interpOpSetC2F.IF2C;
285 I_local2neighbour_e = interpOpSetC2F.IC2F;
286 I_local2neighbour_d = interpOpSetF2C.IC2F;
287 elseif grid_ratio > 1
288 % Local is finer than neighbour
289 I_neighbour2local_e = interpOpSetC2F.IC2F;
290 I_neighbour2local_d = interpOpSetF2C.IC2F;
291 I_local2neighbour_e = interpOpSetF2C.IF2C;
292 I_local2neighbour_d = interpOpSetC2F.IF2C;
293 end
294 otherwise
295 error(['Interpolation type ' obj.interpolation_type ...
296 ' is not available.' ]);
297 end
298
299 else
300 % No interpolation required
301 I_neighbour2local_e = speye(m,m);
302 I_neighbour2local_d = speye(m,m);
303 I_local2neighbour_e = speye(m,m);
304 I_local2neighbour_d = speye(m,m);
305 end
306
307 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d';
308 penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ...
309 -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour';
310
191 end 311 end
192 312
193 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 313 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
194 function [j, nj] = get_boundary_number(obj, boundary) 314 function [j, nj] = get_boundary_number(obj, boundary)
195 315