Mercurial > repos > public > sbplib
changeset 497:4905446f165e feature/quantumTriangles
Added 2D interface to shrodinger
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Sat, 25 Feb 2017 12:44:01 +0100 |
| parents | 437fad4a47e1 |
| children | 324c927d8b1d |
| files | +anim/make_movie.sh +scheme/Schrodinger2dCurve.m |
| diffstat | 1 files changed, 129 insertions(+), 95 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Schrodinger2dCurve.m Fri Feb 24 13:58:18 2017 +0100 +++ b/+scheme/Schrodinger2dCurve.m Sat Feb 25 12:44:01 2017 +0100 @@ -20,7 +20,10 @@ D2_v, D2_u Du, Dv x,y - + b1, b2 + b1_u,b2_v + DU, DV, DUU, DUV, DVU, DVV + e_w, e_e, e_s, e_n du_w, dv_w du_e, dv_e @@ -28,26 +31,29 @@ du_n, dv_n g_1, g_2 c + ind + t_up + a11, a12, a22 m_tot, m_u, m_v - p,p_tau + p Ji end methods - function obj = Schrodinger2dCurve(g ,order, opSet,p,p_tau) + function obj = Schrodinger2dCurve(g ,order, opSet,p) default_arg('opSet',@sbp.D2Variable); default_arg('c', 1); obj.p=p; - obj.p_tau=p_tau; obj.c=1; m = g.size(); obj.m_u = m(1); obj.m_v = m(2); obj.m_tot = g.N(); - + obj.grid=g; + h = g.scaling(); h_u = h(1); h_v = h(2); @@ -80,14 +86,14 @@ obj.Du = kr(obj.D1_u,I_v); obj.Dv = kr(I_u,obj.D1_v); - + obj.H = kr(H_u,H_v); obj.Hi = kr(Hi_u,Hi_v); obj.Hu = kr(H_u,I_v); obj.Hv = kr(I_u,H_v); obj.Hiu = kr(Hi_u,I_v); obj.Hiv = kr(I_u,Hi_v); - + obj.e_w = kr(e_l_u,I_v); obj.e_e = kr(e_r_u,I_v); obj.e_s = kr(I_u,e_l_v); @@ -100,91 +106,88 @@ obj.dv_s = kr(I_u,d1_l_v); obj.du_n = (obj.e_n'*obj.Du)'; obj.dv_n = kr(I_u,d1_r_v); - -% obj.x_u = x_u; -% obj.x_v = x_v; -% obj.y_u = y_u; -% obj.y_v = y_v; - + + obj.DUU = sparse(obj.m_tot); + obj.DVV = sparse(obj.m_tot); + obj.ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot); + obj.m = m; obj.h = [h_u h_v]; obj.order = order; - obj.grid = g; - - + obj.D = @(t)obj.d_fun(t); + obj.variable_update(0); end - - function [D ]= d_fun(obj,t) - % Metric derivatives - ti = parametrization.Ti.points(obj.p.p1(t),obj.p.p2(t),obj.p.p3(t),obj.p.p4(t)); - 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)); - - lcoords=points(obj.grid); - [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); - [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); - x = reshape(obj.xm,obj.m_tot,1); - y = reshape(obj.ym,obj.m_tot,1); - obj.x=x; - obj.y=y; - - x_tau = reshape(x_tau,obj.m_tot,1); - y_tau = reshape(y_tau,obj.m_tot,1); - - x_u = obj.Du*x; - x_v = obj.Dv*x; - y_u = obj.Du*y; - y_v = obj.Dv*y; - - J = x_u.*y_v - x_v.*y_u; - a11 = 1./J.* (x_v.^2 + y_v.^2); - a12 = -1./J .* (x_u.*x_v + y_u.*y_v); - a22 = 1./J .* (x_u.^2 + y_u.^2); - - obj.a11 = a11; - obj.a12 = a12; - obj.a22 = a22; + function [D] = d_fun(obj,t) + % obj.update_vairables(t); In driscretization? + 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)); - % Assemble full operators - L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot); - Duv = obj.Du*L_12*obj.Dv; - Dvu = obj.Dv*L_12*obj.Du; - - Duu = sparse(obj.m_tot); - Dvv = sparse(obj.m_tot); - ind = grid.funcToMatrix(obj.grid, 1:obj.m_tot); - - - for i = 1:obj.m_v - D = obj.D2_u(a11(ind(:,i))); - p = ind(:,i); - Duu(p,p) = D; + end + + + function [D ]= variable_update(obj,t) + % Metric derivatives + if(obj.t_up == t) + return + else + ti = parametrization.Ti.points(obj.p{1}(t),obj.p{2}(t),obj.p{3}(t),obj.p{4}(t)); + ti_tau = parametrization.Ti.points(obj.p{5}(t),obj.p{6}(t),obj.p{7}(t),obj.p{8}(t)); + + lcoords=points(obj.grid); + [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); + [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_u,2)); + x = reshape(obj.xm,obj.m_tot,1); + y = reshape(obj.ym,obj.m_tot,1); + obj.x = x; + obj.y = y; + + x_tau = reshape(x_tau,obj.m_tot,1); + y_tau = reshape(y_tau,obj.m_tot,1); + + x_u = obj.Du*x; + x_v = obj.Dv*x; + y_u = obj.Du*y; + y_v = obj.Dv*y; + + J = x_u.*y_v - x_v.*y_u; + a11 = 1./J.* (x_v.^2 + y_v.^2); + a12 = -1./J .* (x_u.*x_v + y_u.*y_v); + a22 = 1./J .* (x_u.^2 + y_u.^2); + + obj.a11 = a11; + obj.a12 = a12; + obj.a22 = a22; + + % Assemble full operators + L_12 = spdiags(a12, 0, obj.m_tot, obj.m_tot); + obj.DUV = obj.Du*L_12*obj.Dv; + obj.DVU = obj.Dv*L_12*obj.Du; + + + for i = 1:obj.m_v + D = obj.D2_u(a11(obj.ind(:,i))); + p = obj.ind(:,i); + obj.DUU(p,p) = D; + end + + for i = 1:obj.m_u + D = obj.D2_v(a22(obj.ind(i,:))); + p = obj.ind(i,:); + obj.DVV(p,p) = D; + end + + Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); + obj.Ji = Ji; + obj.g_1 = x_v.*y_tau-x_tau.*y_v; + obj.g_2 = x_tau.*y_u - y_tau.*x_u; + + obj.b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot); + obj.b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot); + + obj.b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot); + obj.b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot); + obj.t_up=t; end - - for i = 1:obj.m_u - D = obj.D2_v(a22(ind(i,:))); - p = ind(i,:); - Dvv(p,p) = D; - end - - Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); - obj.Ji=Ji; - obj.g_1 = x_v.*y_tau-x_tau.*y_v; - obj.g_2 = x_tau.*y_u - y_tau.*x_u; - - b1 = spdiags(obj.g_1, 0, obj.m_tot, obj.m_tot); - b2 = spdiags(obj.g_2, 0, obj.m_tot, obj.m_tot); - - b1_u = spdiags(obj.Du*obj.g_1, 0, obj.m_tot, obj.m_tot); - b2_v = spdiags(obj.Dv*obj.g_2, 0, obj.m_tot, obj.m_tot); - - %Add the flux splitting - % D = Ji*(-b1*obj.Du -b2*obj.Dv + 1i*obj.c^2*(Duu + Duv + Dvu + Dvv)); - 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)); - -% obj.gamm_u = h_u*ops_u.borrowing.M.d1; -% obj.gamm_v = h_v*ops_v.borrowing.M.d1; - end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -194,27 +197,56 @@ % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj, boundary) - [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); + function [closure, penalty] = boundary_condition(obj, boundary,~) + [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); a_t = spdiag(coeff_t); - F = (s * d_n' + s * a_t*d_t')'; + a_n = spdiag(coeff_n); + F = (s * a_n *d_n' + s * a_t*d_t')'; tau1 = 1; a = spdiag(g); tau2 = (-1*s*a - abs(a))/4; - penalty_parameter_1 = 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e; - penalty_parameter_2 = halfnorm_inv_n*e*tau2; + penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F*e'*halfnorm_t*e; + penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2; - closure = obj.Ji*obj.c^2 * penalty_parameter_1*e' + obj.Ji* penalty_parameter_2*e'; - penalty = -obj.Ji*obj.c^2 * penalty_parameter_1*e'+ obj.Ji*penalty_parameter_2*e'; + closure = @(t) obj.Ji*obj.c^2 * penalty_parameter_1(t)*e' + obj.Ji* penalty_parameter_2(t)*e'; + penalty = @(t) -obj.Ji*obj.c^2 * penalty_parameter_1(t)*e'+ obj.Ji*penalty_parameter_2(t)*e'; end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + [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); + [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); + + a_n_u = spdiag(coeff_n_u); + a_t_u = spdiag(coeff_t_u); + a_n_v = spdiag(coeff_n_v); + a_t_v = spdiag(coeff_t_v); + + F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; + F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; + + a = spdiag(gamm_u); + + u = obj; + v = neighbour_scheme; + + tau = -1/2*1i; + sig = 1/2*1i; + gamm = (-1*s_u*a)/2; + + 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); + penalty_parameter_2 =@(t) halfnorm_inv_u_n * e_u * (-sig + gamm ); + + closure =@(t) obj.Ji*obj.c^2 * ( penalty_parameter_1(t)*e_u' + penalty_parameter_2(t)*F_u'); + penalty =@(t) obj.Ji*obj.c^2 * (-penalty_parameter_1(t)*e_v' + penalty_parameter_2(t)*F_v'); end - function [e, d_n, d_t, coeff_t, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I] = get_boundary_ops(obj, boundary) + + 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) % gridMatrix = zeros(obj.m(2),obj.m(1)); % gridMatrix(:) = 1:numel(gridMatrix); @@ -230,6 +262,7 @@ I = ind(1,:); coeff_t = obj.a12(I); + coeff_n = obj.a11(I); g=obj.g_1(I); case 'e' e = obj.e_e; @@ -239,6 +272,7 @@ I = ind(end,:); coeff_t = obj.a12(I); + coeff_n = obj.a11(I); g=obj.g_1(I); case 's' e = obj.e_s; @@ -248,6 +282,7 @@ I = ind(:,1)'; coeff_t = obj.a12(I); + coeff_n = obj.a11(I); g=obj.g_2(I); case 'n' e = obj.e_n; @@ -257,6 +292,7 @@ I = ind(:,end)'; coeff_t = obj.a12(I); + coeff_n = obj.a11(I); g=obj.g_2(I); otherwise error('No such boundary: boundary = %s',boundary); @@ -278,7 +314,5 @@ function N = size(obj) N = prod(obj.m); end - - end end \ No newline at end of file