Mercurial > repos > public > sbplib
changeset 498:324c927d8b1d feature/quantumTriangles
chnaged sbp interfacein 1d among many things
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Fri, 03 Mar 2017 16:19:41 +0100 |
| parents | 4905446f165e |
| children | f1465e6aeb26 |
| files | +multiblock/DiffOp.m +multiblock/stitchSchemes.m +scheme/Schrodinger1dCurve.m +scheme/Schrodinger2dCurve.m +scheme/Utux.m |
| diffstat | 5 files changed, 152 insertions(+), 106 deletions(-) [+] |
line wrap: on
line diff
diff -r 4905446f165e -r 324c927d8b1d +multiblock/DiffOp.m --- a/+multiblock/DiffOp.m Sat Feb 25 12:44:01 2017 +0100 +++ b/+multiblock/DiffOp.m Fri Mar 03 16:19:41 2017 +0100 @@ -5,12 +5,12 @@ diffOps D H - + blockmatrixDiv end - + methods - function obj = DiffOp(doHand, grid, order, doParam) + function obj = DiffOp(doHand, grid, order, doParam,timeDep) % doHand -- may either be a function handle or a cell array of % function handles for each grid. The function handle(s) % should be on the form do = doHand(grid, order, ...) @@ -25,13 +25,13 @@ % doHand(..., doParam{i}{:}) Otherwise doParam is sent as % extra parameters to all doHand: doHand(..., doParam{:}) default_arg('doParam', []) - + [getHand, getParam] = parseInput(doHand, grid, doParam); - + nBlocks = grid.nBlocks(); - + obj.order = order; - + % Create the diffOps for each block obj.diffOps = cell(1, nBlocks); for i = 1:nBlocks @@ -42,48 +42,73 @@ end obj.diffOps{i} = h(grid.grids{i}, order, p{:}); end - - + + % Build the norm matrix H = cell(nBlocks, nBlocks); for i = 1:nBlocks H{i,i} = obj.diffOps{i}.H; end obj.H = blockmatrix.toMatrix(H); - - + + % Build the differentiation matrix - obj.blockmatrixDiv = {grid.Ns, grid.Ns}; - D = blockmatrix.zero(obj.blockmatrixDiv); - for i = 1:nBlocks - D{i,i} = obj.diffOps{i}.D; - end - - for i = 1:nBlocks - for j = 1:nBlocks - intf = grid.connections{i,j}; - if isempty(intf) - continue + switch timeDep + case {'n','no','N','No'} + obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + D = blockmatrix.zero(obj.blockmatrixDiv); + for i = 1:nBlocks + D{i,i} = obj.diffOps{i}.D; end - - - [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); - D{i,i} = D{i,i} + ii; - D{i,j} = D{i,j} + ij; - - [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); - D{j,j} = D{j,j} + jj; - D{j,i} = D{j,i} + ji; - end + + for i = 1:nBlocks + for j = 1:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + D{i,i} = D{i,i} + ii; + D{i,j} = D{i,j} + ij; + + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = D{j,j} + jj; + D{j,i} = D{j,i} + ji; + end + end + obj.D = blockmatrix.toMatrix(D); + case {'y','yes','Y','Yes'} + for i = 1:nBlocks + D{i,i} = @(t)obj.diffOps{i}.D(t); + end + + for i = 1:nBlocks + for j = 1:nBlocks + intf = grid.connections{i,j}; + if isempty(intf) + continue + end + + [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2}); + D{i,i} = @(t)D{i,i}(t) + ii(t); + D{i,j} = ij; + + [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1}); + D{j,j} = @(t)D{j,j}(t) + jj(t); + D{j,i} = ji; + end + end + obj.D = D; end - obj.D = blockmatrix.toMatrix(D); - - + + function [getHand, getParam] = parseInput(doHand, grid, doParam) if ~isa(grid, 'multiblock.Grid') error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); end - + if iscell(doHand) && length(doHand) == grid.nBlocks() getHand = @(i)doHand{i}; elseif isa(doHand, 'function_handle') @@ -91,41 +116,41 @@ else error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks'); end - + if isempty(doParam) getParam = @(i){}; return end - + if ~iscell(doParam) getParam = @(i)doParam; return end - + % doParam is a non-empty cell-array - + if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) % doParam is a cell-array of cell-arrays getParam = @(i)doParam{i}; return end - + getParam = @(i)doParam; end end - + function ops = splitOp(obj, op) % Splits a matrix operator into a cell-matrix of matrix operators for % each grid. ops = sparse2cell(op, obj.NNN); end - + function op = getBoundaryOperator(obj, op, boundary) if iscell(boundary) localOpName = [op '_' boundary{2}]; blockId = boundary{1}; localOp = obj.diffOps{blockId}.(localOpName); - + div = {obj.blockmatrixDiv{1}, size(localOp,2)}; blockOp = blockmatrix.zero(div); blockOp{blockId,1} = localOp; @@ -135,7 +160,7 @@ % Boundary är en sträng med en boundary group i. end end - + % Creates the closure and penalty matrix for a given boundary condition, % boundary -- the name of the boundary on the form {id,name} where % id is the number of a block and name is the name of a @@ -143,10 +168,10 @@ function [closure, penalty] = boundary_condition(obj, boundary, type) I = boundary{1}; name = boundary{2}; - + % Get the closure and penaly matrices [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); - + % Expand to matrix for full domain. div = obj.blockmatrixDiv; if ~iscell(blockClosure) @@ -160,7 +185,7 @@ closure{i} = blockmatrix.toMatrix(temp); end end - + div{2} = size(blockPenalty, 2); % Penalty is a column vector if ~iscell(blockPenalty) p = blockmatrix.zero(div); @@ -174,11 +199,11 @@ end end end - + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - + end - + % Size returns the number of degrees of freedom function N = size(obj) N = 0;
diff -r 4905446f165e -r 324c927d8b1d +multiblock/stitchSchemes.m --- a/+multiblock/stitchSchemes.m Sat Feb 25 12:44:01 2017 +0100 +++ b/+multiblock/stitchSchemes.m Fri Mar 03 16:19:41 2017 +0100 @@ -11,7 +11,7 @@ % 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,timeDep) +function [schms, D, H,penalty] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound,timeDep) default_arg('schmParam',[]); default_arg('timeDep','N'); n_blocks = numel(grids); @@ -52,7 +52,11 @@ % Differentiation terms D = cell(n_blocks,n_blocks); + penalty = cell(n_blocks,n_blocks); + + for i = 1:n_blocks + penalty{i,i}= @(t)0; D{i,i} = schms{i}.D; end @@ -104,8 +108,9 @@ continue end - [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); + [closure, penaltyi] = schms{i}.boundary_condition(fn{j},bc{:}); D{i,i} =@(t) D{i,i}(t) + closure(t); + penalty{i,i} = @(t) penalty{i,i}(t) + penaltyi(t); end end
diff -r 4905446f165e -r 324c927d8b1d +scheme/Schrodinger1dCurve.m --- a/+scheme/Schrodinger1dCurve.m Sat Feb 25 12:44:01 2017 +0100 +++ b/+scheme/Schrodinger1dCurve.m Fri Mar 03 16:19:41 2017 +0100 @@ -116,9 +116,9 @@ switch type % Dirichlet boundary condition case {'D','d','dirichlet'} - tau1 = s * 1i*d; + tau1 = @(t) s * 1i*obj.Ji(p,p)^2*d; 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'; + closure = @(t) obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e'; switch class(data) case 'double' @@ -141,15 +141,15 @@ [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 = @(a) -s_u*a(p_u,p_u)/2*e_u; + a1 = s_u* 1/2 * 1i ; + b1 = -s_u* 1/2 * 1i; + gamma = @(a) -obj.Ji*s_u*a(p_u,p_u)/2*e_u; - tau = b1*d_u; - sig = -a1*e_u; + tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u; + sig = b1*e_u; - 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'); + closure = @(t) obj.Hi * (tau(t)*e_u' + sig*obj.Ji(p_u,p_u)^2*d_u' + gamma(obj.a)*e_u'); + penalty = @(t) obj.Hi * (-tau(t)*e_v' - sig*obj.Ji(p_u,p_u)^2*d_v' - gamma(obj.a)*e_v'); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
diff -r 4905446f165e -r 324c927d8b1d +scheme/Schrodinger2dCurve.m --- a/+scheme/Schrodinger2dCurve.m Sat Feb 25 12:44:01 2017 +0100 +++ b/+scheme/Schrodinger2dCurve.m Fri Mar 03 16:19:41 2017 +0100 @@ -134,8 +134,8 @@ 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)); + [obj.xm,obj.ym]= ti.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2)); + [x_tau,y_tau]= ti_tau.map(lcoords(1:obj.m_v:end,1),lcoords(1:obj.m_v,2)); x = reshape(obj.xm,obj.m_tot,1); y = reshape(obj.ym,obj.m_tot,1); obj.x = x; @@ -198,55 +198,59 @@ % 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, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); + [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g,I,Ji] = obj.get_boundary_ops(boundary); - a_t = spdiag(coeff_t); - a_n = spdiag(coeff_n); - F = (s * a_n *d_n' + s * a_t*d_t')'; + a_t = @(t) spdiag(coeff_t(t)); + a_n = @(t) spdiag(coeff_n(t)); + Ji = spdiag(Ji); + + F = @(t)(s * (d_n * Ji)' + s * a_t(t) *d_t')'; tau1 = 1; - a = spdiag(g); - tau2 = (-1*s*a - abs(a))/4; + a = @(t)spdiag(g(t)); + tau2 = @(t) (-1*s*a(t) - abs(a(t)))/4; - 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; + penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e; + penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t); 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'; + 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); + [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,JI_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,JI_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); + a_n_u = @(t) spdiag(coeff_n_u(t)); + a_t_u = @(t) spdiag(coeff_t_u(t)); + a_n_v = @(t) spdiag(coeff_n_v(t)); + a_t_v = @(t) spdiag(coeff_t_v(t)); + + Ji_u = spdiag(JI_u); + Ji_v = spdiag(JI_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')'; + F_u = @(t)((d_n_u*Ji_u)' + a_t_u(t)*d_t_u')'; + F_v = @(t)((d_n_v*Ji_v)' + a_t_v(t)*d_t_v')'; - a = spdiag(gamm_u); - - u = obj; - v = neighbour_scheme; + a = @(t)spdiag(gamm_u(t)); + + tau = s_u*1*1i/2; + sig = -s_u*1*1i/2; + gamm = @(t) (-s_u*a(t))/2; + + penalty_parameter_1 = @(t) halfnorm_inv_u_n*(tau*halfnorm_inv_u_t*F_u(t)*e_u'*halfnorm_u_t*e_u); + penalty_parameter_2 = @(t) halfnorm_inv_u_n * e_u * (sig ); + penalty_parameter_3 = @(t) halfnorm_inv_u_n * e_u * (gamm(t) ); - 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'); + 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'); + 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'); end - 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) + function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g, I,Ji] = get_boundary_ops(obj, boundary) % gridMatrix = zeros(obj.m(2),obj.m(1)); % gridMatrix(:) = 1:numel(gridMatrix); @@ -261,9 +265,9 @@ s = -1; I = ind(1,:); - coeff_t = obj.a12(I); - coeff_n = obj.a11(I); - g=obj.g_1(I); + coeff_t = @(t)obj.a12(I); + coeff_n =@(t) obj.a11(I); + g = @(t)obj.g_1(I); case 'e' e = obj.e_e; d_n = obj.du_e; @@ -271,9 +275,9 @@ s = 1; I = ind(end,:); - coeff_t = obj.a12(I); - coeff_n = obj.a11(I); - g=obj.g_1(I); + coeff_t = @(t)obj.a12(I); + coeff_n = @(t)obj.a11(I); + g = @(t)obj.g_1(I); case 's' e = obj.e_s; d_n = obj.dv_s; @@ -281,9 +285,9 @@ s = -1; I = ind(:,1)'; - coeff_t = obj.a12(I); - coeff_n = obj.a11(I); - g=obj.g_2(I); + coeff_t = @(t)obj.a12(I); + coeff_n = @(t)obj.a11(I); + g = @(t)obj.g_2(I); case 'n' e = obj.e_n; d_n = obj.dv_n; @@ -291,9 +295,9 @@ s = 1; I = ind(:,end)'; - coeff_t = obj.a12(I); - coeff_n = obj.a11(I); - g=obj.g_2(I); + coeff_t = @(t)obj.a12(I); + coeff_n = @(t)obj.a11(I); + g = @(t)obj.g_2(I); otherwise error('No such boundary: boundary = %s',boundary); end @@ -309,10 +313,21 @@ halfnorm_inv_t = obj.Hiu; halfnorm_t = obj.Hu; end + fis = diag(obj.Ji); + Ji = fis(I); end function N = size(obj) N = prod(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