Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 492:6b95a894cbd7 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fixed a bug in the metric coefficients but something is wrong at the boundaries
author Ylva Rydin <ylva.rydin@telia.com>
date Wed, 15 Feb 2017 11:21:04 +0100
parents 26125cfefe11
children 6b8297f66c91
comparison
equal deleted inserted replaced
491:26125cfefe11 492:6b95a894cbd7
17 Hi_u, Hi_v 17 Hi_u, Hi_v
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 22 x,y
23 23
24 e_w, e_e, e_s, e_n 24 e_w, e_e, e_s, e_n
25 du_w, dv_w 25 du_w, dv_w
26 du_e, dv_e 26 du_e, dv_e
27 du_s, dv_s 27 du_s, dv_s
123 lcoords=points(obj.grid); 123 lcoords=points(obj.grid);
124 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); 124 [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
125 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); 125 [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2));
126 x = reshape(obj.xm,obj.m_tot,1); 126 x = reshape(obj.xm,obj.m_tot,1);
127 y = reshape(obj.ym,obj.m_tot,1); 127 y = reshape(obj.ym,obj.m_tot,1);
128 obj.x=x;
129 obj.y=y;
128 130
129 x_tau = reshape(x_tau,obj.m_tot,1); 131 x_tau = reshape(x_tau,obj.m_tot,1);
130 y_tau = reshape(y_tau,obj.m_tot,1); 132 y_tau = reshape(y_tau,obj.m_tot,1);
131 133
132 x_u = obj.Du*x; 134 x_u = obj.Du*x;
165 Dvv(p,p) = D; 167 Dvv(p,p) = D;
166 end 168 end
167 169
168 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); 170 Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot);
169 obj.Ji=Ji; 171 obj.Ji=Ji;
170 obj.g_1 = x_tau.*y_v-y_tau.*x_v; 172 obj.g_1 = x_v.*y_tau-x_tau.*y_v;
171 obj.g_2 = -x_tau.*y_u + y_tau.*x_u; 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);
172 180
173 %Add the flux splitting 181 %Add the flux splitting
174 D = Ji*(obj.g_1.*obj.Du + obj.g_2.*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); %% g_1' och g_2'? 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));
175 184
176 % obj.gamm_u = h_u*ops_u.borrowing.M.d1; 185 % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
177 % obj.gamm_v = h_v*ops_v.borrowing.M.d1; 186 % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
178 187
179 end 188 end
192 F = (s * d_n' + s * a_t*d_t')'; 201 F = (s * d_n' + s * a_t*d_t')';
193 tau1 = 1; 202 tau1 = 1;
194 a = spdiag(g); 203 a = spdiag(g);
195 tau2 = (-1*s*a - abs(a))/4; 204 tau2 = (-1*s*a - abs(a))/4;
196 205
197 penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*e*F'*halfnorm_t*e; 206 penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e;
198 penalty_parameter_2 = halfnorm_inv_n*e*tau2; 207 penalty_parameter_2 = halfnorm_inv_n*e*tau2;
199 208
200 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' +obj.Ji* penalty_parameter_2*e'; 209 closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' + obj.Ji* penalty_parameter_2*e';
201 % penalty = -obj.Ji*obj.c^2 * penalty_parameter_2; 210 penalty = -obj.Ji*obj.c^2 * penalty_parameter_1*e'+ obj.Ji*penalty_parameter_2*e';
202 penalty=0;
203 211
204 end 212 end
205 213
206 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 214 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
207 end 215 end