Mercurial > repos > public > sbplib
changeset 432:eca4ca84cf0a feature/quantumTriangles
Merge in feature/grids
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Wed, 08 Feb 2017 11:10:13 +0100 |
| parents | 5f4540e13f9b (diff) e56dbd9e4196 (current diff) |
| children | b13d44271ead |
| files | +sbp/D2Standard.m |
| diffstat | 3 files changed, 173 insertions(+), 21 deletions(-) [+] |
line wrap: on
line diff
--- a/+sbp/D2Standard.m Tue Feb 07 16:09:02 2017 +0100 +++ b/+sbp/D2Standard.m Wed Feb 08 11:10:13 2017 +0100 @@ -14,7 +14,6 @@ h % Step size x % grid borrowing % Struct with borrowing limits for different norm matrices - end methods @@ -68,6 +67,5 @@ function str = string(obj) str = [class(obj) '_' num2str(obj.order)]; end - end end
--- a/+scheme/Schrodinger.m Tue Feb 07 16:09:02 2017 +0100 +++ b/+scheme/Schrodinger.m Wed Feb 08 11:10:13 2017 +0100 @@ -4,7 +4,8 @@ h % Grid spacing x % Grid order % Order accuracy for the approximation - + grid + D % non-stabalized scheme operator H % Discrete norm M % Derivative norm @@ -23,25 +24,24 @@ % Solving SE in the form u_t = i*u_xx -i*V; function obj = Schrodinger(m,xlim,order,V) default_arg('V',0); - - [x, h] = util.get_grid(xlim{:},m); - - ops = sbp.Ordinary(m,h,order); - - obj.D2 = sparse(ops.derivatives.D2); - obj.H = sparse(ops.norms.H); - obj.Hi = sparse(ops.norms.HI); - obj.M = sparse(ops.norms.M); - obj.e_l = sparse(ops.boundary.e_1); - obj.e_r = sparse(ops.boundary.e_m); - obj.d1_l = sparse(ops.boundary.S_1); - obj.d1_r = sparse(ops.boundary.S_m); + ops = sbp.D2Standard(m,xlim,order); + + obj.x=ops.x; + obj.h=ops.h; + obj.D2 = ops.D2; + obj.H = ops.H; + obj.Hi = ops.HI; + obj.M = ops.M; + obj.e_l = ops.e_l; + obj.e_r = ops.e_r; + obj.d1_l = ops.d1_l; + obj.d1_r = ops.d1_r; if isa(V,'function_handle') - V_vec = V(x); + V_vec = V(obj.x); else - V_vec = x*0 + V; + V_vec = obj.x*0 + V; end V_mat = spdiags(V_vec,0,m,m); @@ -49,10 +49,7 @@ obj.D = 1i * obj.D2 - 1i * V_mat; obj.m = m; - obj.h = h; obj.order = order; - - obj.x = x; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger1dCurve.m Wed Feb 08 11:10:13 2017 +0100 @@ -0,0 +1,157 @@ +classdef Schrodinger1dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + xi % Grid + order % Order accuracy for the approximation + grid + + D % non-stabalized scheme operator + H % Discrete norm + M % Derivative norm + alpha + + V_mat + D1 + D2 + Hi + e_l + e_r + d1_l + d1_r + gamm + end + + methods + % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain; + function obj = Schrodinger1dCurve(m,order,V,constJi) + default_arg('V',0); + default_arg('constJi',false) + xilim={0 1}; + if constJi + ops = sbp.D2Standard(m,xilim,order); + else + ops = sbp.D4Variable(m,xilim,order); + end + + obj.xi=ops.x; + obj.h=ops.h; + obj.D2 = ops.D2; + obj.D1 = ops.D1; + obj.H = ops.H; + obj.Hi = ops.HI; + obj.M = ops.M; + obj.e_l = ops.e_l; + obj.e_r = ops.e_r; + obj.d1_l = ops.d1_l; + obj.d1_r = ops.d1_r; + + + if isa(V,'function_handle') + V_vec = V(obj.x); + else + V_vec = obj.xi*0 + V; + end + + obj.V_mat = spdiags(V_vec,0,m,m); + obj.D = @(a,a_xi,Ji) obj.d_fun(a, a_xi, Ji, constJi); + obj.m = m; + obj.order = order; + end + + + % Closure functions return the opertors appliedo to the own doamin to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % 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 [D] = d_fun(obj,a, a_xi , Ji , constJi) + if constJi + D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*Ji*obj.D2 + 1i*obj.V_mat; + else + D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*obj.D2(diag(Ji)) + 1i*obj.V_mat; + end + end + + function [closure, penalty] = boundary_condition(obj,boundary,type,data) + default_arg('type','dirichlet'); + default_arg('data',0); + + [e,d,s,p] = obj.get_boundary_ops(boundary); + + switch type + % Dirichlet boundary condition + case {'D','d','dirichlet'} + tau1 = s * 1i*d; + tau2 = @(a) (-1*s*a(p,p) - abs(a(p,p)))/4*e; + closure = @(a) obj.Hi*tau1*e' + obj.Hi*tau2(a)*e'; + + switch class(data) + case 'double' + penalty = @(a) -(obj.Hi*tau1*data+obj.Hi*tau2(a)*data); + % case 'function_handle' + % penalty = @(t)-obj.Hi*tau*data(t); + otherwise + error('Wierd data argument!') + end + + % Unknown, boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + 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_u,s_u] = obj.get_boundary_ops(boundary); + % [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + % a = -s_u* 1/2 * 1i ; + % b = a'; + + % tau = b*d_u; + % sig = -a*e_u; + + % closure = obj.Hi * (tau*e_u' + sig*d_u'); + % penalty = obj.Hi * (-tau*e_v' - sig*d_v'); + end + + % Ruturns the boundary ops and sign for the boundary specified by the string boundary. + % The right boundary is considered the positive boundary + function [e,d,s,p] = get_boundary_ops(obj,boundary) + switch boundary + case 'l' + e = obj.e_l; + d = obj.d1_l; + s = -1; + p=1; + case 'r' + e = obj.e_r; + d = obj.d1_r; + s = 1; + p=obj.m; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file