Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 497:4905446f165e feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added 2D interface to shrodinger
author Ylva Rydin <ylva.rydin@telia.com>
date Sat, 25 Feb 2017 12:44:01 +0100
parents 6b8297f66c91
children 324c927d8b1d
comparison
equal deleted inserted replaced
496:437fad4a47e1 497:4905446f165e
18 Hiu, Hiv 18 Hiu, Hiv
19 D1_v, D1_u 19 D1_v, D1_u
20 D2_v, D2_u 20 D2_v, D2_u
21 Du, Dv 21 Du, Dv
22 x,y 22 x,y
23 23 b1, b2
24 b1_u,b2_v
25 DU, DV, DUU, DUV, DVU, DVV
26
24 e_w, e_e, e_s, e_n 27 e_w, e_e, e_s, e_n
25 du_w, dv_w 28 du_w, dv_w
26 du_e, dv_e 29 du_e, dv_e
27 du_s, dv_s 30 du_s, dv_s
28 du_n, dv_n 31 du_n, dv_n
29 g_1, g_2 32 g_1, g_2
30 c 33 c
34 ind
35 t_up
36
31 a11, a12, a22 37 a11, a12, a22
32 m_tot, m_u, m_v 38 m_tot, m_u, m_v
33 p,p_tau 39 p
34 Ji 40 Ji
35 end 41 end
36 42
37 methods 43 methods
38 function obj = Schrodinger2dCurve(g ,order, opSet,p,p_tau) 44 function obj = Schrodinger2dCurve(g ,order, opSet,p)
39 default_arg('opSet',@sbp.D2Variable); 45 default_arg('opSet',@sbp.D2Variable);
40 default_arg('c', 1); 46 default_arg('c', 1);
41 47
42 obj.p=p; 48 obj.p=p;
43 obj.p_tau=p_tau;
44 obj.c=1; 49 obj.c=1;
45 50
46 m = g.size(); 51 m = g.size();
47 obj.m_u = m(1); 52 obj.m_u = m(1);
48 obj.m_v = m(2); 53 obj.m_v = m(2);
49 obj.m_tot = g.N(); 54 obj.m_tot = g.N();
50 55 obj.grid=g;
56
51 h = g.scaling(); 57 h = g.scaling();
52 h_u = h(1); 58 h_u = h(1);
53 h_v = h(2); 59 h_v = h(2);
54 60
55 % Operators 61 % Operators
78 d1_l_v = ops_v.d1_l; 84 d1_l_v = ops_v.d1_l;
79 d1_r_v = ops_v.d1_r; 85 d1_r_v = ops_v.d1_r;
80 86
81 obj.Du = kr(obj.D1_u,I_v); 87 obj.Du = kr(obj.D1_u,I_v);
82 obj.Dv = kr(I_u,obj.D1_v); 88 obj.Dv = kr(I_u,obj.D1_v);
83 89
84 obj.H = kr(H_u,H_v); 90 obj.H = kr(H_u,H_v);
85 obj.Hi = kr(Hi_u,Hi_v); 91 obj.Hi = kr(Hi_u,Hi_v);
86 obj.Hu = kr(H_u,I_v); 92 obj.Hu = kr(H_u,I_v);
87 obj.Hv = kr(I_u,H_v); 93 obj.Hv = kr(I_u,H_v);
88 obj.Hiu = kr(Hi_u,I_v); 94 obj.Hiu = kr(Hi_u,I_v);
89 obj.Hiv = kr(I_u,Hi_v); 95 obj.Hiv = kr(I_u,Hi_v);
90 96
91 obj.e_w = kr(e_l_u,I_v); 97 obj.e_w = kr(e_l_u,I_v);
92 obj.e_e = kr(e_r_u,I_v); 98 obj.e_e = kr(e_r_u,I_v);
93 obj.e_s = kr(I_u,e_l_v); 99 obj.e_s = kr(I_u,e_l_v);
94 obj.e_n = kr(I_u,e_r_v); 100 obj.e_n = kr(I_u,e_r_v);
95 obj.du_w = kr(d1_l_u,I_v); 101 obj.du_w = kr(d1_l_u,I_v);
98 obj.dv_e = (obj.e_e'*obj.Dv)'; 104 obj.dv_e = (obj.e_e'*obj.Dv)';
99 obj.du_s = (obj.e_s'*obj.Du)'; 105 obj.du_s = (obj.e_s'*obj.Du)';
100 obj.dv_s = kr(I_u,d1_l_v); 106 obj.dv_s = kr(I_u,d1_l_v);
101 obj.du_n = (obj.e_n'*obj.Du)'; 107 obj.du_n = (obj.e_n'*obj.Du)';
102 obj.dv_n = kr(I_u,d1_r_v); 108 obj.dv_n = kr(I_u,d1_r_v);
103 109
104 % obj.x_u = x_u; 110 obj.DUU = sparse(obj.m_tot);
105 % obj.x_v = x_v; 111 obj.DVV = sparse(obj.m_tot);
106 % obj.y_u = y_u; 112 obj.ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot);
107 % obj.y_v = y_v; 113
108
109 obj.m = m; 114 obj.m = m;
110 obj.h = [h_u h_v]; 115 obj.h = [h_u h_v];
111 obj.order = order; 116 obj.order = order;
112 obj.grid = g; 117 obj.D = @(t)obj.d_fun(t);
113 118 obj.variable_update(0);
114 119 end
115 end 120
116 121 function [D] = d_fun(obj,t)
117 122 % obj.update_vairables(t); In driscretization?
118 function [D ]= d_fun(obj,t) 123 D = obj.Ji*(-1/2*(obj.b1*obj.Du-obj.b1_u+obj.Du*obj.b1) - 1/2*(obj.b2*obj.Dv - obj.b2_v +obj.Dv*obj.b2) + 1i*obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV));
119 % Metric derivatives 124
120 ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3(t),obj.p.p4(t)); 125 end
121 ti_tau = parametrization.Ti.points(obj.p_tau.p1(t),obj.p_tau.p2(t),obj.p_tau.p3(t),obj.p_tau.p4(t)); 126
122 127
123 lcoords=points(obj.grid); 128 function [D ]= variable_update(obj,t)
124 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); 129 % Metric derivatives
125 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); 130 if(obj.t_up == t)
126 x = reshape(obj.xm,obj.m_tot,1); 131 return
127 y = reshape(obj.ym,obj.m_tot,1); 132 else
128 obj.x=x; 133 ti = parametrization.Ti.points(obj.p{1}(t),obj.p{2}(t),obj.p{3}(t),obj.p{4}(t));
129 obj.y=y; 134 ti_tau = parametrization.Ti.points(obj.p{5}(t),obj.p{6}(t),obj.p{7}(t),obj.p{8}(t));
130 135
131 x_tau = reshape(x_tau,obj.m_tot,1); 136 lcoords=points(obj.grid);
132 y_tau = reshape(y_tau,obj.m_tot,1); 137 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
133 138 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
134 x_u = obj.Du*x; 139 x = reshape(obj.xm,obj.m_tot,1);
135 x_v = obj.Dv*x; 140 y = reshape(obj.ym,obj.m_tot,1);
136 y_u = obj.Du*y; 141 obj.x = x;
137 y_v = obj.Dv*y; 142 obj.y = y;
138 143
139 J = x_u.*y_v - x_v.*y_u; 144 x_tau = reshape(x_tau,obj.m_tot,1);
140 a11 = 1./J.* (x_v.^2 + y_v.^2); 145 y_tau = reshape(y_tau,obj.m_tot,1);
141 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); 146
142 a22 = 1./J .* (x_u.^2 + y_u.^2); 147 x_u = obj.Du*x;
143 148 x_v = obj.Dv*x;
144 obj.a11 = a11; 149 y_u = obj.Du*y;
145 obj.a12 = a12; 150 y_v = obj.Dv*y;
146 obj.a22 = a22; 151
147 152 J = x_u.*y_v - x_v.*y_u;
148 % Assemble full operators 153 a11 = 1./J.* (x_v.^2 + y_v.^2);
149 L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot); 154 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
150 Duv = obj.Du*L_12*obj.Dv; 155 a22 = 1./J .* (x_u.^2 + y_u.^2);
151 Dvu = obj.Dv*L_12*obj.Du; 156
152 157 obj.a11 = a11;
153 Duu = sparse(obj.m_tot); 158 obj.a12 = a12;
154 Dvv = sparse(obj.m_tot); 159 obj.a22 = a22;
155 ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot); 160
156 161 % Assemble full operators
157 162 L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot);
158 for i = 1:obj.m_v 163 obj.DUV = obj.Du*L_12*obj.Dv;
159 D = obj.D2_u(a11(ind(:,i))); 164 obj.DVU = obj.Dv*L_12*obj.Du;
160 p = ind(:,i); 165
161 Duu(p,p) = D; 166
167 for i = 1:obj.m_v
168 D = obj.D2_u(a11(obj.ind(:,i)));
169 p = obj.ind(:,i);
170 obj.DUU(p,p) = D;
171 end
172
173 for i = 1:obj.m_u
174 D = obj.D2_v(a22(obj.ind(i,:)));
175 p = obj.ind(i,:);
176 obj.DVV(p,p) = D;
177 end
178
179 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot);
180 obj.Ji = Ji;
181 obj.g_1 = x_v.*y_tau-x_tau.*y_v;
182 obj.g_2 = x_tau.*y_u - y_tau.*x_u;
183
184 obj.b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot);
185 obj.b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot);
186
187 obj.b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot);
188 obj.b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot);
189 obj.t_up=t;
162 end 190 end
163
164 for i = 1:obj.m_u
165 D = obj.D2_v(a22(ind(i,:)));
166 p = ind(i,:);
167 Dvv(p,p) = D;
168 end
169
170 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot);
171 obj.Ji=Ji;
172 obj.g_1 = x_v.*y_tau-x_tau.*y_v;
173 obj.g_2 = x_tau.*y_u - y_tau.*x_u;
174
175 b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot);
176 b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot);
177
178 b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot);
179 b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot);
180
181 %Add the flux splitting
182 % D = Ji*(-b1*obj.Du -b2*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv));
183 D = Ji*(-1/2*(b1*obj.Du-b1_u+obj.Du*b1) - 1/2*(b2*obj.Dv - b2_v +obj.Dv*b2) + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv));
184
185 % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
186 % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
187
188 end 191 end
189 192
190 % Closure functions return the opertors applied to the own doamin to close the boundary 193 % Closure functions return the opertors applied to the own doamin to close the boundary
191 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 194 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
192 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 195 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
193 % type is a string specifying the type of boundary condition if there are several. 196 % type is a string specifying the type of boundary condition if there are several.
194 % data is a function returning the data that should be applied at the boundary. 197 % data is a function returning the data that should be applied at the boundary.
195 % neighbour_scheme is an instance of Scheme that should be interfaced to. 198 % neighbour_scheme is an instance of Scheme that should be interfaced to.
196 % neighbour_boundary is a string specifying which boundary to interface to. 199 % neighbour_boundary is a string specifying which boundary to interface to.
197 function [closure, penalty] = boundary_condition(obj, boundary) 200 function [closure, penalty] = boundary_condition(obj, boundary,~)
198 [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); 201 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
199 202
200 a_t = spdiag(coeff_t); 203 a_t = spdiag(coeff_t);
201 F = (s * d_n' + s * a_t*d_t')'; 204 a_n = spdiag(coeff_n);
205 F = (s * a_n *d_n' + s * a_t*d_t')';
202 tau1 = 1; 206 tau1 = 1;
203 a = spdiag(g); 207 a = spdiag(g);
204 tau2 = (-1*s*a - abs(a))/4; 208 tau2 = (-1*s*a - abs(a))/4;
205 209
206 penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e; 210 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e;
207 penalty_parameter_2 = halfnorm_inv_n*e*tau2; 211 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2;
208 212
209 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' + obj.Ji* penalty_parameter_2*e'; 213 closure = @(t) obj.Ji*obj.c^2 * penalty_parameter_1(t)*e' + obj.Ji* penalty_parameter_2(t)*e';
210 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1*e'+ obj.Ji*penalty_parameter_2*e'; 214 penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'+ obj.Ji*penalty_parameter_2(t)*e';
211 215
212 end 216 end
213 217
214 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 218 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
215 end 219 % u denotes the solution in the own domain
216 220 % v denotes the solution in the neighbour domain
217 function [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary) 221 [e_u, d_n_u, d_t_u, coeff_t_u, coeff_n_u,s_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t,gamm_u, I_u] = obj.get_boundary_ops(boundary);
222 [e_v, d_n_v, d_t_v, coeff_t_v, coeff_n_v s_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t,gamm_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
223
224 a_n_u = spdiag(coeff_n_u);
225 a_t_u = spdiag(coeff_t_u);
226 a_n_v = spdiag(coeff_n_v);
227 a_t_v = spdiag(coeff_t_v);
228
229 F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')';
230 F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')';
231
232 a = spdiag(gamm_u);
233
234 u = obj;
235 v = neighbour_scheme;
236
237 tau = -1/2*1i;
238 sig = 1/2*1i;
239 gamm = (-1*s_u*a)/2;
240
241 penalty_parameter_1 =@(t) halfnorm_inv_u_n*(e_u*tau + sig*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u);
242 penalty_parameter_2 =@(t) halfnorm_inv_u_n * e_u * (-sig + gamm );
243
244 closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u');
245 penalty =@(t) obj.Ji*obj.c^2 * (-penalty_parameter_1(t)*e_v' + penalty_parameter_2(t)*F_v');
246 end
247
248
249 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary)
218 250
219 % gridMatrix = zeros(obj.m(2),obj.m(1)); 251 % gridMatrix = zeros(obj.m(2),obj.m(1));
220 % gridMatrix(:) = 1:numel(gridMatrix); 252 % gridMatrix(:) = 1:numel(gridMatrix);
221 253
222 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 254 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
228 d_t = obj.dv_w; 260 d_t = obj.dv_w;
229 s = -1; 261 s = -1;
230 262
231 I = ind(1,:); 263 I = ind(1,:);
232 coeff_t = obj.a12(I); 264 coeff_t = obj.a12(I);
265 coeff_n = obj.a11(I);
233 g=obj.g_1(I); 266 g=obj.g_1(I);
234 case 'e' 267 case 'e'
235 e = obj.e_e; 268 e = obj.e_e;
236 d_n = obj.du_e; 269 d_n = obj.du_e;
237 d_t = obj.dv_e; 270 d_t = obj.dv_e;
238 s = 1; 271 s = 1;
239 272
240 I = ind(end,:); 273 I = ind(end,:);
241 coeff_t = obj.a12(I); 274 coeff_t = obj.a12(I);
275 coeff_n = obj.a11(I);
242 g=obj.g_1(I); 276 g=obj.g_1(I);
243 case 's' 277 case 's'
244 e = obj.e_s; 278 e = obj.e_s;
245 d_n = obj.dv_s; 279 d_n = obj.dv_s;
246 d_t = obj.du_s; 280 d_t = obj.du_s;
247 s = -1; 281 s = -1;
248 282
249 I = ind(:,1)'; 283 I = ind(:,1)';
250 coeff_t = obj.a12(I); 284 coeff_t = obj.a12(I);
285 coeff_n = obj.a11(I);
251 g=obj.g_2(I); 286 g=obj.g_2(I);
252 case 'n' 287 case 'n'
253 e = obj.e_n; 288 e = obj.e_n;
254 d_n = obj.dv_n; 289 d_n = obj.dv_n;
255 d_t = obj.du_n; 290 d_t = obj.du_n;
256 s = 1; 291 s = 1;
257 292
258 I = ind(:,end)'; 293 I = ind(:,end)';
259 coeff_t = obj.a12(I); 294 coeff_t = obj.a12(I);
295 coeff_n = obj.a11(I);
260 g=obj.g_2(I); 296 g=obj.g_2(I);
261 otherwise 297 otherwise
262 error('No such boundary: boundary = %s',boundary); 298 error('No such boundary: boundary = %s',boundary);
263 end 299 end
264 300
276 end 312 end
277 313
278 function N = size(obj) 314 function N = size(obj)
279 N = prod(obj.m); 315 N = prod(obj.m);
280 end 316 end
281
282
283 end 317 end
284 end 318 end