comparison +scheme/Schrodinger2dCurve.m @ 515:ba92b27da5a0 feature/quantumTriangles

Tried to make skewsym. Values left on the diagonal when time-dep!
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 07 Jul 2017 14:22:02 +0200
parents 83734c26b8e3
children afff85574ddb
comparison
equal deleted inserted replaced
514:32a24485f3e8 515:ba92b27da5a0
118 obj.variable_update(0); 118 obj.variable_update(0);
119 end 119 end
120 120
121 function [D] = d_fun(obj,t) 121 function [D] = d_fun(obj,t)
122 % obj.update_vairables(t); In driscretization? 122 % obj.update_vairables(t); In driscretization?
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)); 123 % D = obj.Ji*(-1/2*(obj.b1*obj.Du-obj.b1_u+obj.Du*obj.b1) -
124 124 % 1/2*(obj.b2*obj.Dv - obj.b2_v +obj.Dv*obj.b2) +
125 % 1i*obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV)); (ols
126 % not skew sym disc
127
128 D = sqrt(obj.Ji)*(-1/2*(obj.b1*obj.Du + obj.Du*obj.b1) - 1/2*(obj.b2*obj.Dv + obj.Dv*obj.b2) + 1i*obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV))*sqrt(obj.Ji);
125 end 129 end
126 130
127 131
128 function [D ]= variable_update(obj,t) 132 function [D ]= variable_update(obj,t)
129 % Metric derivatives 133 % Metric derivatives
212 tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4; 216 tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4;
213 217
214 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e; 218 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e;
215 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t); 219 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t);
216 220
217 closure = @(t) (obj.Ji*obj.c^2 * penalty_parameter_1(t)*e' + obj.Ji* penalty_parameter_2(t)*e'); 221 closure = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji);
218 penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'- obj.Ji*penalty_parameter_2(t)*e'; 222 penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' - penalty_parameter_2(t)*e')*sqrt(obj.Ji);
219 223
220 end 224 end
221 225
222 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 226 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
223 % u denotes the solution in the own domain 227 % u denotes the solution in the own domain
243 penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u); 247 penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u);
244 penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig ); 248 penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig );
245 penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) ); 249 penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) );
246 250
247 251
248 closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u'); 252 closure =@(t) sqrt(obj.Ji)*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u')*sqrt(obj.Ji);
249 penalty =@(t) obj.Ji*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v'); 253 penalty =@(t) sqrt(obj.Ji)*obj.c^2 * ( -penalty_parameter_1(t)*e_v' - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v')*sqrt(obj.Ji);
250 end 254 end
251 255
252 256
253 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = get_boundary_ops(obj, boundary) 257 function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = get_boundary_ops(obj, boundary)
254 258