|
490
|
1 classdef Schrodinger2dCurve < scheme.Scheme
|
|
|
2 properties
|
|
|
3 m % Number of points in each direction, possibly a vector
|
|
|
4 h % Grid spacing
|
|
|
5
|
|
|
6 grid
|
|
|
7
|
|
|
8 order % Order accuracy for the approximation
|
|
|
9
|
|
|
10 D % non-stabalized scheme operator
|
|
|
11 M % Derivative norm
|
|
|
12 H % Discrete norm
|
|
|
13 Hi
|
|
|
14 H_u, H_v % Norms in the x and y directions
|
|
|
15 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
|
|
|
16 Hi_u, Hi_v
|
|
|
17 Hiu, Hiv
|
|
|
18 D1_v D1_u
|
|
|
19 D2_v D2_u
|
|
|
20 Du Dv
|
|
|
21
|
|
|
22
|
|
|
23 e_w, e_e, e_s, e_n
|
|
|
24 du_w, dv_w
|
|
|
25 du_e, dv_e
|
|
|
26 du_s, dv_s
|
|
|
27 du_n, dv_n
|
|
|
28 g_1
|
|
|
29 g_2
|
|
|
30
|
|
|
31 p,p_tau
|
|
|
32 end
|
|
|
33
|
|
|
34 methods
|
|
|
35 function obj = Schrodinger2dCurve(g ,order, opSet, p,p_tau)
|
|
|
36 default_arg('opSet',@sbp.D2Variable);
|
|
|
37 default_arg('c', 1);
|
|
|
38
|
|
|
39 assert(isa(g, 'grid.Curvilinear'))
|
|
|
40
|
|
|
41 obj.p=p;
|
|
|
42 obj.p_tau=p_tau;
|
|
|
43 obj.c=1;
|
|
|
44
|
|
|
45 m = g.size();
|
|
|
46 m_u = m(1);
|
|
|
47 m_v = m(2);
|
|
|
48 m_tot = g.N();
|
|
|
49
|
|
|
50 h = g.scaling();
|
|
|
51 h_u = h(1);
|
|
|
52 h_v = h(2);
|
|
|
53
|
|
|
54 % Operators
|
|
|
55 ops_u = opSet(m_u, {0, 1}, order);
|
|
|
56 ops_v = opSet(m_v, {0, 1}, order);
|
|
|
57
|
|
|
58 I_u = speye(m_u);
|
|
|
59 I_v = speye(m_v);
|
|
|
60
|
|
|
61 D1_u = ops_u.D1;
|
|
|
62
|
|
|
63 H_u = ops_u.H;
|
|
|
64 Hi_u = ops_u.HI;
|
|
|
65 e_l_u = ops_u.e_l;
|
|
|
66 e_r_u = ops_u.e_r;
|
|
|
67 d1_l_u = ops_u.d1_l;
|
|
|
68 d1_r_u = ops_u.d1_r;
|
|
|
69
|
|
|
70 obj.D1_v = ops_v.D1;
|
|
|
71 obj.D2_v = ops_v.D2;
|
|
|
72 H_v = ops_v.H;
|
|
|
73 Hi_v = ops_v.HI;
|
|
|
74 e_l_v = ops_v.e_l;
|
|
|
75 e_r_v = ops_v.e_r;
|
|
|
76 d1_l_v = ops_v.d1_l;
|
|
|
77 d1_r_v = ops_v.d1_r;
|
|
|
78
|
|
|
79 obj.Du = kr(D1_u,I_v);
|
|
|
80 obj.Dv = kr(I_u,D1_v);
|
|
|
81
|
|
|
82 obj.H = kr(H_u,H_v);
|
|
|
83 obj.Hi = kr(Hi_u,Hi_v);
|
|
|
84 obj.Hu = kr(H_u,I_v);
|
|
|
85 obj.Hv = kr(I_u,H_v);
|
|
|
86 obj.Hiu = kr(Hi_u,I_v);
|
|
|
87 obj.Hiv = kr(I_u,Hi_v);
|
|
|
88
|
|
|
89 obj.e_w = kr(e_l_u,I_v);
|
|
|
90 obj.e_e = kr(e_r_u,I_v);
|
|
|
91 obj.e_s = kr(I_u,e_l_v);
|
|
|
92 obj.e_n = kr(I_u,e_r_v);
|
|
|
93 obj.du_w = kr(d1_l_u,I_v);
|
|
|
94 obj.dv_w = (obj.e_w'*Dv)';
|
|
|
95 obj.du_e = kr(d1_r_u,I_v);
|
|
|
96 obj.dv_e = (obj.e_e'*Dv)';
|
|
|
97 obj.du_s = (obj.e_s'*Du)';
|
|
|
98 obj.dv_s = kr(I_u,d1_l_v);
|
|
|
99 obj.du_n = (obj.e_n'*Du)';
|
|
|
100 obj.dv_n = kr(I_u,d1_r_v);
|
|
|
101
|
|
|
102 % obj.x_u = x_u;
|
|
|
103 % obj.x_v = x_v;
|
|
|
104 % obj.y_u = y_u;
|
|
|
105 % obj.y_v = y_v;
|
|
|
106
|
|
|
107 obj.m = m;
|
|
|
108 obj.h = [h_u h_v];
|
|
|
109 obj.order = order;
|
|
|
110 obj.grid = g;
|
|
|
111
|
|
|
112
|
|
|
113 end
|
|
|
114
|
|
|
115
|
|
|
116 function [D ]= d_fun(obj,t)
|
|
|
117 % Metric derivatives
|
|
|
118 ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3,obj.p.p4);
|
|
|
119 ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3,obj.p_tau.p4);
|
|
|
120
|
|
|
121 coords = parametrization.ti.map();
|
|
|
122 coords_tau = parametrization.ti_tau.map();
|
|
|
123 x = coords(:,1);
|
|
|
124 y = coords(:,2);
|
|
|
125
|
|
|
126 x_tau = coords_tau(:,1);
|
|
|
127 y_tau = coords_tau(:,2);
|
|
|
128
|
|
|
129 x_u = obj.Du*x;
|
|
|
130 x_v = obj.Dv*x;
|
|
|
131 y_u = obj.Du*y;
|
|
|
132 y_v = obj.Dv*y;
|
|
|
133
|
|
|
134 J = x_u.*y_v - x_v.*y_u;
|
|
|
135 a11 = 1./J.* (x_v.^2 + y_v.^2);
|
|
|
136 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
|
|
|
137 a22 = 1./J .* (x_u.^2 + y_u.^2);
|
|
|
138 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
|
|
|
139
|
|
|
140 % Assemble full operators
|
|
|
141 L_12 = spdiags(a12, 0, m_tot, m_tot);
|
|
|
142 Duv = obj.Du*L_12*obj.Dv;
|
|
|
143 Dvu = obj.Dv*L_12*obj.Du;
|
|
|
144
|
|
|
145 Duu = sparse(m_tot);
|
|
|
146 Dvv = sparse(m_tot);
|
|
|
147 ind = grid.funcToMatrix(g, 1:m_tot);
|
|
|
148
|
|
|
149 for i = 1:m_v
|
|
|
150 D = D2_u(a11(ind(:,i)));
|
|
|
151 p = ind(:,i);
|
|
|
152 Duu(p,p) = D;
|
|
|
153 end
|
|
|
154
|
|
|
155 for i = 1:m_u
|
|
|
156 D = D2_v(a22(ind(i,:)));
|
|
|
157 p = ind(i,:);
|
|
|
158 Dvv(p,p) = D;
|
|
|
159 end
|
|
|
160
|
|
|
161
|
|
|
162 J = spdiags(J, 0, m_tot, m_tot);
|
|
|
163 Ji = spdiags(1./J, 0, m_tot, m_tot);
|
|
|
164 obj.g_1 = x_tau.*y_v-y_tau.*x_v;
|
|
|
165 obj.g_2 = -x_tau.*y_u + y_tau.*x_u;
|
|
|
166
|
|
|
167 %Add the flux splitting
|
|
|
168 D = Ji*(obj.g_1*obj.Du + obj.g_2*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv));
|
|
|
169
|
|
|
170 % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
|
|
|
171 % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
|
|
|
172
|
|
|
173 end
|
|
|
174
|
|
|
175 % Closure functions return the opertors applied to the own doamin to close the boundary
|
|
|
176 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
|
|
|
177 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
|
|
|
178 % type is a string specifying the type of boundary condition if there are several.
|
|
|
179 % data is a function returning the data that should be applied at the boundary.
|
|
|
180 % neighbour_scheme is an instance of Scheme that should be interfaced to.
|
|
|
181 % neighbour_boundary is a string specifying which boundary to interface to.
|
|
|
182 function [closure, penalty] = boundary_condition(obj, boundary, parameter)
|
|
|
183 default_arg('type','neumann');
|
|
|
184 default_arg('parameter', []);
|
|
|
185
|
|
|
186 % v denotes the solution in the neighbour domain
|
|
|
187 tuning = 1.2;
|
|
|
188 % tuning = 20.2;
|
|
|
189 [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
|
|
|
190
|
|
|
191 a_n = spdiag(coeff_n);
|
|
|
192 a_t = spdiag(coeff_t);
|
|
|
193
|
|
|
194 F = (s * a_n * d_n' + s * a_t*d_t')';
|
|
|
195
|
|
|
196 u = obj;
|
|
|
197
|
|
|
198 b1 = gamm*u.lambda./u.a11.^2;
|
|
|
199 b2 = gamm*u.lambda./u.a22.^2;
|
|
|
200
|
|
|
201 tau1 = -1./b1 - 1./b2;
|
|
|
202 tau1 = tuning * spdiag(tau1);
|
|
|
203 sig1 = 1;
|
|
|
204
|
|
|
205 a = e'*g;
|
|
|
206 tau2 = (-1*s*a - abs(a))/4;
|
|
|
207
|
|
|
208 penalty_parameter_1 = halfnorm_inv_n*(tau1 + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e;
|
|
|
209 penalty_parameter_2 = halfnorm_inv_n(tau2)*e;
|
|
|
210
|
|
|
211 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e';
|
|
|
212 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1;
|
|
|
213
|
|
|
214 end
|
|
|
215
|
|
|
216 function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor,g] = get_boundary_ops(obj, boundary)
|
|
|
217
|
|
|
218 % gridMatrix = zeros(obj.m(2),obj.m(1));
|
|
|
219 % gridMatrix(:) = 1:numel(gridMatrix);
|
|
|
220
|
|
|
221 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
|
|
|
222
|
|
|
223 switch boundary
|
|
|
224 case 'w'
|
|
|
225 e = obj.e_w;
|
|
|
226 d_n = obj.du_w;
|
|
|
227 d_t = obj.dv_w;
|
|
|
228 s = -1;
|
|
|
229
|
|
|
230 I = ind(1,:);
|
|
|
231 coeff_n = obj.a11(I);
|
|
|
232 coeff_t = obj.a12(I);
|
|
|
233 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
|
|
|
234 case 'e'
|
|
|
235 e = obj.e_e;
|
|
|
236 d_n = obj.du_e;
|
|
|
237 d_t = obj.dv_e;
|
|
|
238 s = 1;
|
|
|
239
|
|
|
240 I = ind(end,:);
|
|
|
241 coeff_n = obj.a11(I);
|
|
|
242 coeff_t = obj.a12(I);
|
|
|
243 scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
|
|
|
244 case 's'
|
|
|
245 e = obj.e_s;
|
|
|
246 d_n = obj.dv_s;
|
|
|
247 d_t = obj.du_s;
|
|
|
248 s = -1;
|
|
|
249
|
|
|
250 I = ind(:,1)';
|
|
|
251 coeff_n = obj.a22(I);
|
|
|
252 coeff_t = obj.a12(I);
|
|
|
253 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
|
|
|
254 case 'n'
|
|
|
255 e = obj.e_n;
|
|
|
256 d_n = obj.dv_n;
|
|
|
257 d_t = obj.du_n;
|
|
|
258 s = 1;
|
|
|
259
|
|
|
260 I = ind(:,end)';
|
|
|
261 coeff_n = obj.a22(I);
|
|
|
262 coeff_t = obj.a12(I);
|
|
|
263 scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
|
|
|
264 otherwise
|
|
|
265 error('No such boundary: boundary = %s',boundary);
|
|
|
266 end
|
|
|
267
|
|
|
268 switch boundary
|
|
|
269 case {'w','e'}
|
|
|
270 halfnorm_inv_n = obj.Hiu;
|
|
|
271 halfnorm_inv_t = obj.Hiv;
|
|
|
272 halfnorm_t = obj.Hv;
|
|
|
273 gamm = obj.gamm_u;
|
|
|
274 g=obj.g_1;
|
|
|
275 case {'s','n'}
|
|
|
276 halfnorm_inv_n = obj.Hiv;
|
|
|
277 halfnorm_inv_t = obj.Hiu;
|
|
|
278 halfnorm_t = obj.Hu;
|
|
|
279 gamm = obj.gamm_v;
|
|
|
280 g=obj.g_2;
|
|
|
281 end
|
|
|
282 end
|
|
|
283
|
|
|
284 function N = size(obj)
|
|
|
285 N = prod(obj.m);
|
|
|
286 end
|
|
|
287
|
|
|
288
|
|
|
289 end
|
|
|
290 end |
