Mercurial > repos > public > sbplib
changeset 496:437fad4a47e1 feature/quantumTriangles
Add 1d interface for shrodinger in 1d
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Fri, 24 Feb 2017 13:58:18 +0100 |
| parents | b91d23271481 |
| children | 4905446f165e |
| files | +multiblock/stitchSchemes.m +scheme/Schrodinger1dCurve.m |
| diffstat | 2 files changed, 121 insertions(+), 56 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/stitchSchemes.m Fri Feb 24 09:04:02 2017 +0100 +++ b/+multiblock/stitchSchemes.m Fri Feb 24 13:58:18 2017 +0100 @@ -11,9 +11,9 @@ % Output parameters are cell arrays and cell matrices. % % Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) -function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound) +function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound,timeDep) default_arg('schmParam',[]); - + default_arg('timeDep','N'); n_blocks = numel(grids); % Creating Schemes @@ -47,46 +47,81 @@ for i = 1:n_blocks H{i,i} = schms{i}.H; end - + %% Total system matrix - + % Differentiation terms D = cell(n_blocks,n_blocks); for i = 1:n_blocks D{i,i} = schms{i}.D; end - + % Boundary penalty terms - for i = 1:n_blocks - if ~isstruct(bound{i}) - continue - end - - fn = fieldnames(bound{i}); - for j = 1:length(fn); - bc = bound{i}.(fn{j}); - if isempty(bc) - continue + switch timeDep + case {'n','no','N','No'} + for i = 1:n_blocks + if ~isstruct(bound{i}) + continue + end + + fn = fieldnames(bound{i}); + for j = 1:length(fn) + bc = bound{i}.(fn{j}); + if isempty(bc) + continue + end + + [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); + D{i,i} = D{i,i}+closure; + end end - - [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); - D{i,i} = D{i,i}+closure; - end - end - - % Interface penalty terms - for i = 1:n_blocks - for j = 1:n_blocks - intf = conn{i,j}; - if isempty(intf) - continue + + % Interface penalty terms + for i = 1:n_blocks + for j = 1:n_blocks + intf = conn{i,j}; + if isempty(intf) + continue + end + + [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2}); + D{i,i} = D{i,i} + uu; + D{i,j} = uv; + D{j,j} = D{j,j} + vv; + D{j,i} = vu; + end end - - [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2}); - D{i,i} = D{i,i} + uu; - D{i,j} = uv; - D{j,j} = D{j,j} + vv; - D{j,i} = vu; - end - end -end \ No newline at end of file + case {'y','yes','Y','Yes'} + for i = 1:n_blocks + if ~isstruct(bound{i}) + continue + end + + fn = fieldnames(bound{i}); + for j = 1:length(fn) + bc = bound{i}.(fn{j}); + if isempty(bc) + continue + end + + [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); + D{i,i} =@(t) D{i,i}(t) + closure(t); + end + end + + % Interface penalty terms + for i = 1:n_blocks + for j = 1:n_blocks + intf = conn{i,j}; + if isempty(intf) + continue + end + + [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2}); + D{i,i} = @(t) D{i,i}(t) + uu(t); + D{i,j} = uv; + D{j,j} = @(t)D{j,j}(t) + vv(t); + D{j,i} = vu; + end + end + end \ No newline at end of file
--- a/+scheme/Schrodinger1dCurve.m Fri Feb 24 09:04:02 2017 +0100 +++ b/+scheme/Schrodinger1dCurve.m Fri Feb 24 13:58:18 2017 +0100 @@ -10,6 +10,15 @@ H % Discrete norm M % Derivative norm alpha + x_r + x_l + ddt_x_r + ddt_x_l + a + a_xi + Ji + t_up + x V_mat D1 @@ -24,16 +33,22 @@ methods % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain; - function obj = Schrodinger1dCurve(g,m,order,V,constJi) + function obj = Schrodinger1dCurve(g,order,boundaries,V,constJi) default_arg('V',0); default_arg('constJi',false) xilim={0 1}; + m = N(g); if constJi ops = sbp.D2Standard(m,xilim,order); else ops = sbp.D4Variable(m,xilim,order); end + obj.x_l = boundaries{1}; + obj.x_r = boundaries{2}; + obj.ddt_x_l = boundaries{3}; + obj.ddt_x_r = boundaries{4}; + obj.xi=ops.x; obj.h=ops.h; obj.D2 = ops.D2; @@ -47,15 +62,14 @@ obj.d1_r = ops.d1_r; obj.grid = g; - 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.V_mat = spdiags(V_vec,0,m,m); + obj.D = @(t) obj.d_fun(t); obj.m = m; obj.order = order; end @@ -69,11 +83,27 @@ % 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; + function [D] = d_fun(obj,t) + % obj.update_vairables(t); In driscretization? + D= obj.Ji*(-0.5*(obj.D1*obj.a - obj.a_xi + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat); + + end + + + function [] = variable_update(obj,t) + if (t == obj.t_up) + return else - D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*obj.D2(diag(Ji)) + 1i*obj.V_mat; + x_r = obj.x_r(t); + x_l = obj.x_l(t); + ddt_x_r = obj. ddt_x_r(t); + ddt_x_l = obj.ddt_x_l(t); + obj.x = obj.xi*(x_r -x_l) + x_l; + obj.a = sparse(diag((-ddt_x_l*( x_r - x_l) - (obj.x-x_l)*(ddt_x_r-ddt_x_l))/(x_r-x_l))); + + obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x))); + obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x))); + obj.t_up = t; end end @@ -87,12 +117,12 @@ % 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'; + tau2 = @(t) obj.Ji*(-1*s*obj.a(p,p) - abs(obj.a(p,p)))/4*e; + closure = @(t) obj.Ji*obj.Hi*tau1*e' + obj.Hi*tau2(obj.a)*e'; switch class(data) case 'double' - penalty = @(a) -(obj.Hi*tau1*data+obj.Hi*tau2(a)*data); + penalty = @(t) -obj.Ji*(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data); % case 'function_handle' % penalty = @(t)-obj.Hi*tau*data(t); otherwise @@ -108,18 +138,18 @@ 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,p_u] = obj.get_boundary_ops(boundary); - [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary); + [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - a1 = -s_u* 1/2 * 1i ; - b1 = a1'; - gamma = -1*s*a(p_u,p_u)/2; + a1 = -s_u* 1/2 * 1i ; + b1 = a1'; + gamma = @(a) -s_u*a(p_u,p_u)/2*e_u; - tau = b1*d_u; - sig = -a1*e_u; + tau = b1*d_u; + sig = -a1*e_u; - closure = @(a) obj.Hi * (tau*e_u' + sig*d_u' + gamma(a)); - penalty = @(a) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(a)); + closure = @(t) obj.Hi * (tau*e_u' + sig*d_u' + gamma(obj.a)*e_u'); + penalty = @(t) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(obj.a)*e_v'); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary.