Mercurial > repos > public > sbplib

comparison +grid/Schrodinger2dCurve.m @ 490:b13d44271ead feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Schrodinger2dCurve Added
author Ylva Rydin <ylva.rydin@telia.com>
date Thu, 09 Feb 2017 11:41:21 +0100
parents
children
comparison
equal deleted inserted replaced
432:eca4ca84cf0a 490:b13d44271ead
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