Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 1033:037f203b9bf5 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with branch feature/advectioRV to utilize the +rv package
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:44:12 +0100
parents 3dd7f87c9a1b
children 78db023a7fe3
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
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 end
37
38 methods
39
40 function obj = Schrodinger2d(g ,order, a, opSet)
41 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
42 default_arg('a',1);
43 dim = 2;
44
45 assertType(g, 'grid.Cartesian');
46 if isa(a, 'function_handle')
47 a = grid.evalOn(g, a);
48 a = spdiag(a);
49 end
50
51 m = g.size();
52 m_tot = g.N();
53
54 h = g.scaling();
55 xlim = {g.x{1}(1), g.x{1}(end)};
56 ylim = {g.x{2}(1), g.x{2}(end)};
57 lim = {xlim, ylim};
58
59 % 1D operators
60 ops = cell(dim,1);
61 for i = 1:dim
62 ops{i} = opSet{i}(m(i), lim{i}, order);
63 end
64
65 I = cell(dim,1);
66 D1 = cell(dim,1);
67 D2 = cell(dim,1);
68 H = cell(dim,1);
69 Hi = cell(dim,1);
70 e_l = cell(dim,1);
71 e_r = cell(dim,1);
72 d1_l = cell(dim,1);
73 d1_r = cell(dim,1);
74
75 for i = 1:dim
76 I{i} = speye(m(i));
77 D1{i} = ops{i}.D1;
78 D2{i} = ops{i}.D2;
79 H{i} = ops{i}.H;
80 Hi{i} = ops{i}.HI;
81 e_l{i} = ops{i}.e_l;
82 e_r{i} = ops{i}.e_r;
83 d1_l{i} = ops{i}.d1_l;
84 d1_r{i} = ops{i}.d1_r;
85 end
86
87 % Constant coeff D2
88 for i = 1:dim
89 D2{i} = D2{i}(ones(m(i),1));
90 end
91
92 %====== Assemble full operators ========
93 obj.D1 = cell(dim,1);
94 obj.D2 = cell(dim,1);
95 obj.e_l = cell(dim,1);
96 obj.e_r = cell(dim,1);
97 obj.d1_l = cell(dim,1);
98 obj.d1_r = cell(dim,1);
99
100 % D1
101 obj.D1{1} = kron(D1{1},I{2});
102 obj.D1{2} = kron(I{1},D1{2});
103
104 % Boundary operators
105 obj.e_l{1} = kron(e_l{1},I{2});
106 obj.e_l{2} = kron(I{1},e_l{2});
107 obj.e_r{1} = kron(e_r{1},I{2});
108 obj.e_r{2} = kron(I{1},e_r{2});
109
110 obj.d1_l{1} = kron(d1_l{1},I{2});
111 obj.d1_l{2} = kron(I{1},d1_l{2});
112 obj.d1_r{1} = kron(d1_r{1},I{2});
113 obj.d1_r{2} = kron(I{1},d1_r{2});
114
115 % D2
116 obj.D2{1} = kron(D2{1},I{2});
117 obj.D2{2} = kron(I{1},D2{2});
118
119 % Quadratures
120 obj.H = kron(H{1},H{2});
121 obj.Hi = inv(obj.H);
122 obj.H_boundary = cell(dim,1);
123 obj.H_boundary{1} = H{2};
124 obj.H_boundary{2} = H{1};
125
126 % Differentiation matrix D (without SAT)
127 D2 = obj.D2;
128 D = sparse(m_tot,m_tot);
129 for j = 1:dim
130 D = D + a*1i*D2{j};
131 end
132 obj.D = D;
133 %=========================================%
134
135 % Misc.
136 obj.m = m;
137 obj.h = h;
138 obj.order = order;
139 obj.grid = g;
140 obj.dim = dim;
141 obj.a = a;
142 obj.e_w = obj.e_l{1};
143 obj.e_e = obj.e_r{1};
144 obj.e_s = obj.e_l{2};
145 obj.e_n = obj.e_r{2};
146 obj.d_w = obj.d1_l{1};
147 obj.d_e = obj.d1_r{1};
148 obj.d_s = obj.d1_l{2};
149 obj.d_n = obj.d1_r{2};
150
151 end
152
153
154 % Closure functions return the operators applied to the own domain to close the boundary
155 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
156 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
157 % type is a string specifying the type of boundary condition.
158 % data is a function returning the data that should be applied at the boundary.
159 % neighbour_scheme is an instance of Scheme that should be interfaced to.
160 % neighbour_boundary is a string specifying which boundary to interface to.
161 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
162 default_arg('type','Neumann');
163 default_arg('parameter', []);
164
165 % j is the coordinate direction of the boundary
166 % nj: outward unit normal component.
167 % nj = -1 for west, south, bottom boundaries
168 % nj = 1 for east, north, top boundaries
169 [j, nj] = obj.get_boundary_number(boundary);
170 switch nj
171 case 1
172 e = obj.e_r;
173 d = obj.d1_r;
174 case -1
175 e = obj.e_l;
176 d = obj.d1_l;
177 end
178
179 Hi = obj.Hi;
180 H_gamma = obj.H_boundary{j};
181 a = e{j}'*obj.a*e{j};
182
183 switch type
184
185 % Dirichlet boundary condition
186 case {'D','d','dirichlet','Dirichlet'}
187 closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' );
188 penalty = -nj*Hi*d{j}*a*1i*H_gamma;
189
190 % Free boundary condition
191 case {'N','n','neumann','Neumann'}
192 closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' );
193 penalty = nj*Hi*e{j}*a*1i*H_gamma;
194
195 % Unknown boundary condition
196 otherwise
197 error('No such boundary condition: type = %s',type);
198 end
199 end
200
201 % type Struct that specifies the interface coupling.
202 % Fields:
203 % -- interpolation: type of interpolation, default 'none'
204 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
205
206 defaultType.interpolation = 'none';
207 default_struct('type', defaultType);
208
209 switch type.interpolation
210 case {'none', ''}
211 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
212 case {'op','OP'}
213 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
214 otherwise
215 error('Unknown type of interpolation: %s ', type.interpolation);
216 end
217 end
218
219 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
220 % u denotes the solution in the own domain
221 % v denotes the solution in the neighbour domain
222
223 % Get boundary operators
224 [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
225 [e, d, H_gamma] = obj.get_boundary_ops(boundary);
226 Hi = obj.Hi;
227 a = obj.a;
228
229 % Get outward unit normal component
230 [~, n] = obj.get_boundary_number(boundary);
231
232 Hi = obj.Hi;
233 sigma = -n*1i*a/2;
234 tau = -n*(1i*a)'/2;
235
236 closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d';
237 penalty = -tau*Hi*d*H_gamma*e_neighbour' ...
238 -sigma*Hi*e*H_gamma*d_neighbour';
239
240 end
241
242 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
243
244 % User can request special interpolation operators by specifying type.interpOpSet
245 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
246 interpOpSet = type.interpOpSet;
247
248 % u denotes the solution in the own domain
249 % v denotes the solution in the neighbour domain
250 [e_v, d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
251 [e_u, d_u, H_gamma] = obj.get_boundary_ops(boundary);
252 Hi = obj.Hi;
253 a = obj.a;
254
255 % Get outward unit normal component
256 [~, n] = obj.get_boundary_number(boundary);
257
258 % Find the number of grid points along the interface
259 m_u = size(e_u, 2);
260 m_v = size(e_v, 2);
261
262 % Build interpolation operators
263 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
264 Iu2v = intOps.Iu2v;
265 Iv2u = intOps.Iv2u;
266
267 sigma = -n*1i*a/2;
268 tau = -n*(1i*a)'/2;
269
270 closure = tau*Hi*d_u*H_gamma*e_u' + sigma*Hi*e_u*H_gamma*d_u';
271 penalty = -tau*Hi*d_u*H_gamma*Iv2u.good*e_v' ...
272 -sigma*Hi*e_u*H_gamma*Iv2u.bad*d_v';
273
274 end
275
276 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
277 function [j, nj] = get_boundary_number(obj, boundary)
278
279 switch boundary
280 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
281 j = 1;
282 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
283 j = 2;
284 otherwise
285 error('No such boundary: boundary = %s',boundary);
286 end
287
288 switch boundary
289 case {'w','W','west','West','s','S','south','South'}
290 nj = -1;
291 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
292 nj = 1;
293 end
294 end
295
296 % Returns the boundary ops and sign for the boundary specified by the string boundary.
297 % The right boundary is considered the positive boundary
298 function [e, d, H_b] = get_boundary_ops(obj, boundary)
299
300 switch boundary
301 case 'w'
302 e = obj.e_w;
303 d = obj.d_w;
304 H_b = obj.H_boundary{1};
305 case 'e'
306 e = obj.e_e;
307 d = obj.d_e;
308 H_b = obj.H_boundary{1};
309 case 's'
310 e = obj.e_s;
311 d = obj.d_s;
312 H_b = obj.H_boundary{2};
313 case 'n'
314 e = obj.e_n;
315 d = obj.d_n;
316 H_b = obj.H_boundary{2};
317 otherwise
318 error('No such boundary: boundary = %s',boundary);
319 end
320 end
321
322 function N = size(obj)
323 N = prod(obj.m);
324 end
325 end
326 end