Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 707:0de70ec8bf60 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
merge with feature/optim
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 10 Nov 2017 14:22:56 +0100
parents e89715fe6a6e
children acb58769610e
comparison
equal deleted inserted replaced
696:7c16b5af8d98 707:0de70ec8bf60
36 36
37 a11, a12, a22 37 a11, a12, a22
38 m_tot, m_u, m_v 38 m_tot, m_u, m_v
39 p 39 p
40 Ji, J 40 Ji, J
41 Hamiltonian
41 end 42 end
42 43
43 methods 44 methods
44 function obj = Schrodinger2dCurve(g ,order, opSet,p) 45 function obj = Schrodinger2dCurve(g ,order, opSet,p)
45 default_arg('opSet',@sbp.D2Variable); 46 default_arg('opSet',@sbp.D2Variable);
65 I_u = speye(obj.m_u); 66 I_u = speye(obj.m_u);
66 I_v = speye(obj.m_v); 67 I_v = speye(obj.m_v);
67 68
68 obj.D1_u = ops_u.D1; 69 obj.D1_u = ops_u.D1;
69 obj.D2_u = ops_u.D2; 70 obj.D2_u = ops_u.D2;
70 71
71 H_u = ops_u.H; 72 H_u = ops_u.H;
72 Hi_u = ops_u.HI; 73 Hi_u = ops_u.HI;
73 e_l_u = ops_u.e_l; 74 e_l_u = ops_u.e_l;
74 e_r_u = ops_u.e_r; 75 e_r_u = ops_u.e_r;
75 d1_l_u = ops_u.d1_l; 76 d1_l_u = ops_u.d1_l;
113 114
114 obj.m = m; 115 obj.m = m;
115 obj.h = [h_u h_v]; 116 obj.h = [h_u h_v];
116 obj.order = order; 117 obj.order = order;
117 obj.D = @(t)obj.d_fun(t); 118 obj.D = @(t)obj.d_fun(t);
119 obj.Hamiltonian = @(t)obj.h_fun(t);
118 obj.variable_update(0); 120 obj.variable_update(0);
119 end 121 end
120 122
121 function [D] = d_fun(obj,t) 123 function [D,n] = d_fun(obj,t)
122 % obj.update_vairables(t); In driscretization? 124 % obj.update_vairables(t); In driscretization?
123 % D = obj.Ji*(-1/2*(obj.b1*obj.Du-obj.b1_u+obj.Du*obj.b1) - 125 % D = obj.Ji*(-1/2*(obj.b1*obj.Du-obj.b1_u+obj.Du*obj.b1) -
124 % 1/2*(obj.b2*obj.Dv - obj.b2_v +obj.Dv*obj.b2) + 126 % 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 127 % 1i*obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV)); (ols
126 % not skew sym disc 128 % not skew sym disc
127 129
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); 130 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);
129 end 131 end
130 132
133 function [Hamiltonian] = h_fun(obj,t)
134 Hamiltonian = -sqrt(obj.Ji)*(obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV))*sqrt(obj.Ji);
135 end
131 136
132 function [D ]= variable_update(obj,t) 137 function [D ]= variable_update(obj,t)
133 % Metric derivatives 138 % Metric derivatives
134 if(obj.t_up == t) 139 if(obj.t_up == t)
135 return 140 return
202 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 207 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
203 % type is a string specifying the type of boundary condition if there are several. 208 % type is a string specifying the type of boundary condition if there are several.
204 % data is a function returning the data that should be applied at the boundary. 209 % data is a function returning the data that should be applied at the boundary.
205 % neighbour_scheme is an instance of Scheme that should be interfaced to. 210 % neighbour_scheme is an instance of Scheme that should be interfaced to.
206 % neighbour_boundary is a string specifying which boundary to interface to. 211 % neighbour_boundary is a string specifying which boundary to interface to.
207 function [closure, penalty] = boundary_condition(obj, boundary,~) 212 function [closure, penalty,closureHamiltonian,penaltyHamiltonian] = boundary_condition(obj, boundary,~)
208 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); 213 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
209 214
210 a_t = @(t) spdiag(coeff_t(t)); 215 a_t = @(t) spdiag(coeff_t(t));
211 a_n = @(t) spdiag(coeff_n(t)); 216 a_n = @(t) spdiag(coeff_n(t));
212 217
219 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e; 224 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e;
220 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t); 225 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t);
221 226
222 closure = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); 227 closure = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji);
223 penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); 228 penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji);
229
230 closureHamiltonian = @(t) 1i*sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji);
231 penaltyHamiltonian = @(t) -1i*sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji);
224 232
225 end 233 end
226 234
227 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 235 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
228 % u denotes the solution in the own domain 236 % u denotes the solution in the own domain