Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 756:f891758ad7a4 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/utux2d.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Jun 2018 14:30:45 -0700
parents f4595f14d696
children 459eeb99130f
comparison
equal deleted inserted replaced
755:14f0058356f2 756:f891758ad7a4
1 classdef Schrodinger2d < scheme.Scheme
2
3 % Discretizes the Laplacian with constant coefficent,
4 % in the Schrödinger equation way (i.e., the discretization matrix is not necessarily
5 % definite)
6 % u_t = a*i*Laplace u
7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction.
9
10 properties
11 m % Number of points in each direction, possibly a vector
12 h % Grid spacing
13
14 grid
15 dim
16
17 order % Order of accuracy for the approximation
18
19 % Diagonal matrix for variable coefficients
20 a % Constant coefficient
21
22 D % Total operator
23 D1 % First derivatives
24
25 % Second derivatives
26 D2
27
28 H, Hi % Inner products
29 e_l, e_r
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
33
34 H_boundary % Boundary inner products
35
36 interpolation_type % MC or AWW
37
38 end
39
40 methods
41
42 function obj = Schrodinger2d(g ,order, a, opSet, interpolation_type)
43 default_arg('interpolation_type','AWW');
44 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
45 default_arg('a',1);
46 dim = 2;
47
48 assert(isa(g, 'grid.Cartesian'))
49 if isa(a, 'function_handle')
50 a = grid.evalOn(g, a);
51 a = spdiag(a);
52 end
53
54 m = g.size();
55 m_tot = g.N();
56
57 h = g.scaling();
58 xlim = {g.x{1}(1), g.x{1}(end)};
59 ylim = {g.x{2}(1), g.x{2}(end)};
60 lim = {xlim, ylim};
61
62 % 1D operators
63 ops = cell(dim,1);
64 for i = 1:dim
65 ops{i} = opSet{i}(m(i), lim{i}, order);
66 end
67
68 I = cell(dim,1);
69 D1 = cell(dim,1);
70 D2 = cell(dim,1);
71 H = cell(dim,1);
72 Hi = cell(dim,1);
73 e_l = cell(dim,1);
74 e_r = cell(dim,1);
75 d1_l = cell(dim,1);
76 d1_r = cell(dim,1);
77
78 for i = 1:dim
79 I{i} = speye(m(i));
80 D1{i} = ops{i}.D1;
81 D2{i} = ops{i}.D2;
82 H{i} = ops{i}.H;
83 Hi{i} = ops{i}.HI;
84 e_l{i} = ops{i}.e_l;
85 e_r{i} = ops{i}.e_r;
86 d1_l{i} = ops{i}.d1_l;
87 d1_r{i} = ops{i}.d1_r;
88 end
89
90 % Constant coeff D2
91 for i = 1:dim
92 D2{i} = D2{i}(ones(m(i),1));
93 end
94
95 %====== Assemble full operators ========
96 obj.D1 = cell(dim,1);
97 obj.D2 = cell(dim,1);
98 obj.e_l = cell(dim,1);
99 obj.e_r = cell(dim,1);
100 obj.d1_l = cell(dim,1);
101 obj.d1_r = cell(dim,1);
102
103 % D1
104 obj.D1{1} = kron(D1{1},I{2});
105 obj.D1{2} = kron(I{1},D1{2});
106
107 % Boundary operators
108 obj.e_l{1} = kron(e_l{1},I{2});
109 obj.e_l{2} = kron(I{1},e_l{2});
110 obj.e_r{1} = kron(e_r{1},I{2});
111 obj.e_r{2} = kron(I{1},e_r{2});
112
113 obj.d1_l{1} = kron(d1_l{1},I{2});
114 obj.d1_l{2} = kron(I{1},d1_l{2});
115 obj.d1_r{1} = kron(d1_r{1},I{2});
116 obj.d1_r{2} = kron(I{1},d1_r{2});
117
118 % D2
119 obj.D2{1} = kron(D2{1},I{2});
120 obj.D2{2} = kron(I{1},D2{2});
121
122 % Quadratures
123 obj.H = kron(H{1},H{2});
124 obj.Hi = inv(obj.H);
125 obj.H_boundary = cell(dim,1);
126 obj.H_boundary{1} = H{2};
127 obj.H_boundary{2} = H{1};
128
129 % Differentiation matrix D (without SAT)
130 D2 = obj.D2;
131 D = sparse(m_tot,m_tot);
132 for j = 1:dim
133 D = D + a*1i*D2{j};
134 end
135 obj.D = D;
136 %=========================================%
137
138 % Misc.
139 obj.m = m;
140 obj.h = h;
141 obj.order = order;
142 obj.grid = g;
143 obj.dim = dim;
144 obj.a = a;
145 obj.e_w = obj.e_l{1};
146 obj.e_e = obj.e_r{1};
147 obj.e_s = obj.e_l{2};
148 obj.e_n = obj.e_r{2};
149 obj.d_w = obj.d1_l{1};
150 obj.d_e = obj.d1_r{1};
151 obj.d_s = obj.d1_l{2};
152 obj.d_n = obj.d1_r{2};
153 obj.interpolation_type = interpolation_type;
154
155 end
156
157
158 % Closure functions return the operators applied to the own domain to close the boundary
159 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
160 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
161 % type is a string specifying the type of boundary condition.
162 % data is a function returning the data that should be applied at the boundary.
163 % neighbour_scheme is an instance of Scheme that should be interfaced to.
164 % neighbour_boundary is a string specifying which boundary to interface to.
165 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
166 default_arg('type','Neumann');
167 default_arg('parameter', []);
168
169 % j is the coordinate direction of the boundary
170 % nj: outward unit normal component.
171 % nj = -1 for west, south, bottom boundaries
172 % nj = 1 for east, north, top boundaries
173 [j, nj] = obj.get_boundary_number(boundary);
174 switch nj
175 case 1
176 e = obj.e_r;
177 d = obj.d1_r;
178 case -1
179 e = obj.e_l;
180 d = obj.d1_l;
181 end
182
183 Hi = obj.Hi;
184 H_gamma = obj.H_boundary{j};
185 a = e{j}'*obj.a*e{j};
186
187 switch type
188
189 % Dirichlet boundary condition
190 case {'D','d','dirichlet','Dirichlet'}
191 closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' );
192 penalty = -nj*Hi*d{j}*a*1i*H_gamma;
193
194 % Free boundary condition
195 case {'N','n','neumann','Neumann'}
196 closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' );
197 penalty = nj*Hi*e{j}*a*1i*H_gamma;
198
199 % Unknown boundary condition
200 otherwise
201 error('No such boundary condition: type = %s',type);
202 end
203 end
204
205 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
206 % u denotes the solution in the own domain
207 % v denotes the solution in the neighbour domain
208 % Get neighbour boundary operator
209
210 [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary);
211 [coord, n] = get_boundary_number(obj, 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;
227 a = obj.a;
228
229 switch coord_nei
230 case 1
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
256 Hi = obj.Hi;
257 sigma = -n*1i*a/2;
258 tau = -n*(1i*a)'/2;
259
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';
312 penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ...
313 -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour';
314
315 end
316
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)
319
320 switch boundary
321 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
322 j = 1;
323 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
324 j = 2;
325 otherwise
326 error('No such boundary: boundary = %s',boundary);
327 end
328
329 switch boundary
330 case {'w','W','west','West','s','S','south','South'}
331 nj = -1;
332 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
333 nj = 1;
334 end
335 end
336
337 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
338 function [return_op] = get_boundary_operator(obj, op, boundary)
339
340 switch boundary
341 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
342 j = 1;
343 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
344 j = 2;
345 otherwise
346 error('No such boundary: boundary = %s',boundary);
347 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
369
370 function N = size(obj)
371 N = prod(obj.m);
372 end
373 end
374 end