Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2d.m @ 718:71aa5828cbbf feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Schrödinger scheme for 2d single block. Will develop to multiblock with interpolation.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 03 Mar 2018 16:18:33 -0800
parents
children b3f8fb9cefd2
comparison
equal deleted inserted replaced
717:8e4274ee6dd8 718:71aa5828cbbf
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
32 H_boundary % Boundary inner products
33
34 end
35
36 methods
37
38 function obj = Schrodinger2d(g ,order, a, opSet)
39 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
40 default_arg('a',1);
41 dim = 2;
42
43 assert(isa(g, 'grid.Cartesian'))
44
45 m = g.size();
46 m_tot = g.N();
47
48 h = g.scaling();
49 xlim = {g.x{1}(1), g.x{1}(end)};
50 ylim = {g.x{2}(1), g.x{2}(end)};
51 lim = {xlim, ylim};
52
53 % 1D operators
54 ops = cell(dim,1);
55 for i = 1:dim
56 ops{i} = opSet{i}(m(i), lim{i}, order);
57 end
58
59 I = cell(dim,1);
60 D1 = cell(dim,1);
61 D2 = cell(dim,1);
62 H = cell(dim,1);
63 Hi = cell(dim,1);
64 e_l = cell(dim,1);
65 e_r = cell(dim,1);
66 d1_l = cell(dim,1);
67 d1_r = cell(dim,1);
68
69 for i = 1:dim
70 I{i} = speye(m(i));
71 D1{i} = ops{i}.D1;
72 D2{i} = ops{i}.D2;
73 H{i} = ops{i}.H;
74 Hi{i} = ops{i}.HI;
75 e_l{i} = ops{i}.e_l;
76 e_r{i} = ops{i}.e_r;
77 d1_l{i} = ops{i}.d1_l;
78 d1_r{i} = ops{i}.d1_r;
79 end
80
81 % Constant coeff D2
82 for i = 1:dim
83 D2{i} = D2{i}(ones(m(i),1));
84 end
85
86 %====== Assemble full operators ========
87 obj.D1 = cell(dim,1);
88 obj.D2 = cell(dim,1);
89 obj.e_l = cell(dim,1);
90 obj.e_r = cell(dim,1);
91 obj.d1_l = cell(dim,1);
92 obj.d1_r = cell(dim,1);
93
94 % D1
95 obj.D1{1} = kron(D1{1},I{2});
96 obj.D1{2} = kron(I{1},D1{2});
97
98 % Boundary operators
99 obj.e_l{1} = kron(e_l{1},I{2});
100 obj.e_l{2} = kron(I{1},e_l{2});
101 obj.e_r{1} = kron(e_r{1},I{2});
102 obj.e_r{2} = kron(I{1},e_r{2});
103
104 obj.d1_l{1} = kron(d1_l{1},I{2});
105 obj.d1_l{2} = kron(I{1},d1_l{2});
106 obj.d1_r{1} = kron(d1_r{1},I{2});
107 obj.d1_r{2} = kron(I{1},d1_r{2});
108
109 % D2
110 obj.D2{1} = kron(D2{1},I{2});
111 obj.D2{2} = kron(I{1},D2{2});
112
113 % Quadratures
114 obj.H = kron(H{1},H{2});
115 obj.Hi = inv(obj.H);
116 obj.H_boundary = cell(dim,1);
117 obj.H_boundary{1} = H{2};
118 obj.H_boundary{2} = H{1};
119
120 % Differentiation matrix D (without SAT)
121 D2 = obj.D2;
122 D = sparse(m_tot,m_tot);
123 for j = 1:dim
124 D = D + a*1i*D2{j};
125 end
126 obj.D = D;
127 %=========================================%
128
129 % Misc.
130 obj.m = m;
131 obj.h = h;
132 obj.order = order;
133 obj.grid = g;
134 obj.dim = dim;
135 obj.a = a;
136
137 end
138
139
140 % Closure functions return the operators applied to the own domain to close the boundary
141 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
142 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
143 % type is a string specifying the type of boundary condition.
144 % data is a function returning the data that should be applied at the boundary.
145 % neighbour_scheme is an instance of Scheme that should be interfaced to.
146 % neighbour_boundary is a string specifying which boundary to interface to.
147 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
148 default_arg('type','Neumann');
149 default_arg('parameter', []);
150
151 % j is the coordinate direction of the boundary
152 % nj: outward unit normal component.
153 % nj = -1 for west, south, bottom boundaries
154 % nj = 1 for east, north, top boundaries
155 [j, nj] = obj.get_boundary_number(boundary);
156 switch nj
157 case 1
158 e = obj.e_r;
159 d = obj.d1_r;
160 case -1
161 e = obj.e_l;
162 d = obj.d1_l;
163 end
164
165 Hi = obj.Hi;
166 H_gamma = obj.H_boundary{j};
167 a = obj.a;
168
169 switch type
170
171 % Dirichlet boundary condition
172 case {'D','d','dirichlet','Dirichlet'}
173 closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' );
174 penalty = -nj*Hi*d{j}*a*1i*H_gamma;
175
176 % Free boundary condition
177 case {'N','n','neumann','Neumann'}
178 closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' );
179 penalty = nj*Hi*e{j}*a*1i*H_gamma;
180
181 % Unknown boundary condition
182 otherwise
183 error('No such boundary condition: type = %s',type);
184 end
185 end
186
187 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
188 % u denotes the solution in the own domain
189 % v denotes the solution in the neighbour domain
190 error('Interface not implemented');
191 end
192
193 % 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)
195
196 switch boundary
197 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
198 j = 1;
199 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
200 j = 2;
201 otherwise
202 error('No such boundary: boundary = %s',boundary);
203 end
204
205 switch boundary
206 case {'w','W','west','West','s','S','south','South'}
207 nj = -1;
208 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
209 nj = 1;
210 end
211 end
212
213 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
214 function [return_op] = get_boundary_operator(obj, op, boundary)
215
216 switch boundary
217 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
218 j = 1;
219 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
220 j = 2;
221 otherwise
222 error('No such boundary: boundary = %s',boundary);
223 end
224
225 switch op
226 case 'e'
227 switch boundary
228 case {'w','W','west','West','s','S','south','South'}
229 return_op = obj.e_l{j};
230 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
231 return_op = obj.e_r{j};
232 end
233 case 'd'
234 switch boundary
235 case {'w','W','west','West','s','S','south','South'}
236 return_op = obj.d1_l{j};
237 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
238 return_op = obj.d1_r{j};
239 end
240 otherwise
241 error(['No such operator: operator = ' op]);
242 end
243
244 end
245
246 function N = size(obj)
247 N = prod(obj.m);
248 end
249 end
250 end