Mercurial > repos > public > sbplib
changeset 713:348d5bcf7daf feature/quantumTriangles
Merge with feature/frids
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Tue, 20 Feb 2018 15:00:30 +0100 |
| parents | c9147e05d228 (diff) facb8ffb5c34 (current diff) |
| children | |
| files | |
| diffstat | 19 files changed, 1114 insertions(+), 59 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DiffOpTimeDep.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,212 @@ +classdef DiffOpTimeDep < scheme.Scheme + properties + grid + order + diffOps + D + H + + blockmatrixDiv + end + + methods + function obj = DiffOpTimeDep(doHand, grid, order, doParam) + % 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, ...) + % Additional parameters for each doHand may be provided in + % the doParam input. + % grid -- a multiblock grid + % order -- integer specifying the order of accuracy + % doParam -- may either be a cell array or a cell array of cell arrays + % for each block. If it is a cell array with length equal + % to the number of blocks then each element is sent to the + % corresponding function handle as extra parameters: + % doHand(..., doParam{i}{:}) Otherwise doParam is sent as + % extra parameters to all doHand: doHand(..., doParam{:}) + default_arg('doParam', []) + + [getHand, getParam] = parseInput(doHand, grid, doParam); + obj.grid = grid; + nBlocks = grid.nBlocks(); + + obj.order = order; + + % Create the diffOps for each block + obj.diffOps = cell(1, nBlocks); + for i = 1:nBlocks + h = getHand(i); + p = getParam(i); + if ~iscell(p) + p = {p}; + 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 + for j = 1:nBlocks + D{i,j} = @(t) 0; + D{j,i} = @(t) 0; + end + end + + + 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} = @(t) D{i,j}(t) + ij(t); + + [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} = @(t) D{j,i}(t) + ji(t); + end + end + obj.D = 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') + getHand = @(i)doHand; + 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 + + % Get a boundary operator specified by opName for the given boundary/BoundaryGroup + function op = getBoundaryOperator(obj, opName, boundary) + switch class(boundary) + case 'cell' + localOpName = [opName '_' 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; + op = blockmatrix.toMatrix(blockOp); + return + case 'multiblock.BoundaryGroup' + op = sparse(size(obj.D,1),0); + for i = 1:length(boundary) + op = [op, obj.getBoundaryOperator(opName, boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') + 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 + % boundary of that block example: {1,'s'} or {3,'w'}. It + % can also be a boundary group + function [closure, penalty] = boundary_condition(obj, boundary, type) + switch class(boundary) + case 'cell' + [closure, penalty] = obj.singleBoundaryCondition(boundary, type); + case 'multiblock.BoundaryGroup' + nBlocks = obj.grid.nBlocks(); + [n,m] = size(obj.D{1,1}(0)); + %closure = sparse(n,m); + %penalty = sparse(n,0); + % closure =@(t)0; + % penalty = @(t)0; + for i = 1:nBlocks + for j = 1:nBlocks + closure{j,i} = @(t)0*sparse(n,m); + penalty{j,i} = @(t)0*sparse(n,1); + end + end + + + for i = 1:length(boundary) + boundaryPart = boundary{i}; + block = boundaryPart{1}; + [closurePart, penaltyPart] = obj.boundary_condition(boundaryPart, type); + closure{block,block} = @(t)closure{block,block}(t) + closurePart(t); + penalty{block,block} = @(t)penalty{block,block}(t) + penaltyPart(t); + end + otherwise + error('Unknown boundary indentifier') + end + + end + + function [blockClosure, blockPenalty] = singleBoundaryCondition(obj, boundary, type) + I = boundary{1}; + name = boundary{2}; + % Get the closure and penaly matrices + [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type); + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('not implemented') + end + + % Size returns the number of degrees of freedom + function N = size(obj) + N = 0; + for i = 1:length(obj.diffOps) + N = N + obj.diffOps{i}.size(); + end + end + end +end
--- a/+multiblock/stitchSchemes.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+multiblock/stitchSchemes.m Tue Feb 20 15:00:30 2018 +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,penalty,J] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound,timeDep) default_arg('schmParam',[]); - + default_arg('timeDep','N'); n_blocks = numel(grids); % Creating Schemes @@ -49,44 +49,84 @@ end %% Total system matrix - + % 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 - + % 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 + + % 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 - - [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 + 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, 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 - - [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 + + % 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/+noname/animate.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+noname/animate.m Tue Feb 20 15:00:30 2018 +0100 @@ -84,7 +84,7 @@ if ~do_step - pause + % pause anim.animate(@G, Tstart, Tend, time_modifier); else pause
--- a/+sbp/D1Upwind.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+sbp/D1Upwind.m Tue Feb 20 15:00:30 2018 +0100 @@ -53,7 +53,7 @@ obj.borrowing = []; end - + function str = string(obj) str = [class(obj) '_' num2str(obj.order)]; end
--- a/+sbp/D2Standard.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+sbp/D2Standard.m Tue Feb 20 15:00:30 2018 +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/Hypsyst3d.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+scheme/Hypsyst3d.m Tue Feb 20 15:00:30 2018 +0100 @@ -7,6 +7,7 @@ X, Y, Z% Values of x and y for each grid point Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces order % Order accuracy for the approximation + grid D % non-stabalized scheme operator A, B, C, E % Symbolic coefficient matrices
--- a/+scheme/Hypsyst3dCurve.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+scheme/Hypsyst3dCurve.m Tue Feb 20 15:00:30 2018 +0100 @@ -28,6 +28,7 @@ I_xi,I_eta,I_zeta, I_N,onesN e_w, e_e, e_s, e_n, e_b, e_t index_w, index_e,index_s,index_n, index_b, index_t + grid params %parameters for the coeficient matrice end @@ -215,6 +216,8 @@ obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta); obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); + Hi = obj.Hxii*obj.Hetai*obj.Hzetai; + obj.H = inv(Hi); obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1); obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1);
--- a/+scheme/Schrodinger.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+scheme/Schrodinger.m Tue Feb 20 15:00:30 2018 +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 @@ -21,27 +22,29 @@ methods % Solving SE in the form u_t = i*u_xx -i*V; - function obj = Schrodinger(m,xlim,order,V) + function obj = Schrodinger(g,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); + obj.grid = g; + m = N(obj.grid); + obj.x = points(obj.grid); + ops = sbp.D2Standard(m,{obj.x(1) obj.x(end)},order); + + + 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 +52,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 Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,183 @@ +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 + x_r + x_l + ddt_x_r + ddt_x_l + a + a_xi + Ji + J + t_up + x + + 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(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; + 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; + 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 = @(t) obj.d_fun(t); + 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,t) + obj.variable_update(t); % In driscretization? + D = sqrt(obj.Ji)*(-0.5*(obj.D1*obj.a + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat)*sqrt(obj.Ji); +% D = (-0.5*(obj.D1*obj.a -obj.a_xi+ obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat); + % D= obj.Ji*(-sqrt(obj.a)*obj.D1*sqrt(obj.a) + 0.5*obj.a_xi + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat); + end + + + function [] = variable_update(obj,t) + if (t == obj.t_up) + return + else + 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.J = sparse(x_r -x_l); + obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x))); + obj.t_up = t; + 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 = @(t) s * 1i*obj.Ji(p,p)*d; + tau2 = @(t) (1*s*obj.a(p,p))/2*e; + closure = @(t)obj.Hi*sqrt(obj.Ji)*(tau1(t)*e' + tau2(obj.a)*e')*sqrt(obj.Ji); + penalty = @(t) -obj.Hi*sqrt(obj.Ji)*(tau1(t)*e' + tau2(obj.a)*e')*sqrt(obj.Ji); + % 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,p_u] = obj.get_boundary_ops(boundary); + [e_v,d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + a1 = s_u* 1/2 * 1i ; + b1 = -s_u* 1/2 * 1i; + gamma = @(a) s_u*a(p_u,p_u)/2*e_u; + + tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u; + sig = b1*e_u; + + closure = @(t) obj.Hi * (tau(t)*e_u' + sig*obj.Ji(p_u,p_u)^2*d_u' + obj.Ji(p_u,p_u)*gamma(obj.a)*e_u'); + penalty = @(t) obj.Hi * (-tau(t)*e_v' - sig*obj.Ji(p_u,p_u)^2*d_v' - obj.Ji(p_u,p_u)*gamma(obj.a)*e_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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger2dCurve.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,345 @@ +classdef Schrodinger2dCurve < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + xm, ym + + order % Order accuracy for the approximation + + D % non-stabalized scheme operator + M % Derivative norm + H % Discrete norm + Hi + H_u, H_v % Norms in the x and y directions + Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. + Hi_u, Hi_v + Hiu, Hiv + D1_v, D1_u + 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 + du_s, dv_s + du_n, dv_n + g_1, g_2 + c + ind + t_up + + a11, a12, a22 + m_tot, m_u, m_v + p + Ji, J + Hamiltonian + end + + methods + function obj = Schrodinger2dCurve(g ,order, opSet,p) + default_arg('opSet',@sbp.D2Variable); + default_arg('c', 1); + + obj.p=p; + 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); + + % Operators + ops_u = opSet(obj.m_u, {0, 1}, order); + ops_v = opSet(obj.m_v, {0, 1}, order); + + I_u = speye(obj.m_u); + I_v = speye(obj.m_v); + + obj.D1_u = ops_u.D1; + obj.D2_u = ops_u.D2; + + H_u = ops_u.H; + Hi_u = ops_u.HI; + e_l_u = ops_u.e_l; + e_r_u = ops_u.e_r; + d1_l_u = ops_u.d1_l; + d1_r_u = ops_u.d1_r; + + obj.D1_v = ops_v.D1; + obj.D2_v = ops_v.D2; + H_v = ops_v.H; + Hi_v = ops_v.HI; + e_l_v = ops_v.e_l; + e_r_v = ops_v.e_r; + d1_l_v = ops_v.d1_l; + d1_r_v = ops_v.d1_r; + + 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); + obj.e_n = kr(I_u,e_r_v); + obj.du_w = kr(d1_l_u,I_v); + obj.dv_w = (obj.e_w'*obj.Dv)'; + obj.du_e = kr(d1_r_u,I_v); + obj.dv_e = (obj.e_e'*obj.Dv)'; + obj.du_s = (obj.e_s'*obj.Du)'; + 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.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.D = @(t)obj.d_fun(t); + obj.Hamiltonian = @(t)obj.h_fun(t); + obj.variable_update(0); + end + + function [D,n] = 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)); (ols + % not skew sym disc + + 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); + end + + function [Hamiltonian] = h_fun(obj,t) + Hamiltonian = -sqrt(obj.Ji)*(obj.c^2*(obj.DUU + obj.DUV + obj.DVU + obj.DVV))*sqrt(obj.Ji); + 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_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; + 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); + obj.J = spdiags(J, 0, obj.m_tot, obj.m_tot); + + Ji = spdiags(1./J, 0, obj.m_tot, obj.m_tot); + obj.Ji = Ji; + + a11 = Ji* (x_v.^2 + y_v.^2); + a12 = -Ji* (x_u.*x_v + y_u.*y_v); + a22 = Ji* (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 + + 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 + end + + % Closure functions return the opertors applied 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 [closure, penalty,closureHamiltonian,penaltyHamiltonian] = 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 = @(t) spdiag(coeff_t(t)); + a_n = @(t) spdiag(coeff_n(t)); + + F = @(t)(s * a_n(t)*d_n' + s * a_t(t) *d_t')'; + tau1 = 1; + a = @(t)spdiag(g(t)); + tau2 = @(t) (1*s*a(t))/2; + + 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) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); + penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); + + closureHamiltonian = @(t) 1i*sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji); + penaltyHamiltonian = @(t) -1i*sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji); + + 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,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,Ji_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + + 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)); + + + F_u = @(t)(a_n_u(t)*d_n_u' + a_t_u(t)*d_t_u')'; + F_v = @(t)(a_n_v(t)*d_n_v' + a_t_v(t)*d_t_v')'; + + 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) ); + + + closure =@(t) sqrt(Ji_u)*obj.c^2 * ( penalty_parameter_1(t)*e_u'... + + penalty_parameter_2(t)*F_u(t)' + penalty_parameter_3(t)*e_u')*sqrt(Ji_u); + penalty =@(t) sqrt(Ji_u)*obj.c^2 * ( -penalty_parameter_1(t)*e_v'... + - penalty_parameter_2(t)*F_v(t)' - penalty_parameter_3(t)*e_v')*sqrt(Ji_v); + end + + + function [e, d_n, d_t, coeff_t,coeff_n, s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g,Ji] = get_boundary_ops(obj, boundary) + + % gridMatrix = zeros(obj.m(2),obj.m(1)); + % gridMatrix(:) = 1:numel(gridMatrix); + + ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); + + switch boundary + case 'w' + e = obj.e_w; + d_n = obj.du_w; + d_t = obj.dv_w; + s = -1; + + I = ind(1,:); + 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; + d_t = obj.dv_e; + s = 1; + + I = ind(end,:); + 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; + d_t = obj.du_s; + s = -1; + + I = ind(:,1)'; + coeff_t = @(t)obj.a12(I); + coeff_n = @(t)obj.a22(I); + g = @(t)obj.g_2(I); + case 'n' + e = obj.e_n; + d_n = obj.dv_n; + d_t = obj.du_n; + s = 1; + + I = ind(:,end)'; + coeff_t = @(t)obj.a12(I); + coeff_n = @(t)obj.a22(I); + g = @(t)obj.g_2(I); + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','e'} + halfnorm_inv_n = obj.Hiu; + halfnorm_inv_t = obj.Hiv; + halfnorm_t = obj.Hv; + + case {'s','n'} + halfnorm_inv_n = obj.Hiv; + halfnorm_inv_t = obj.Hiu; + halfnorm_t = obj.Hu; + end + Ji = obj.Ji; + 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
--- a/+scheme/Utux.m Tue Dec 19 17:07:07 2017 +0100 +++ b/+scheme/Utux.m Tue Feb 20 15:00:30 2018 +0100 @@ -12,6 +12,7 @@ Hi e_l e_r + grid v0 end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/Arnoldi.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,37 @@ +function y = Arnoldi(A,v,t,tol,m) + +n = size (v,1); +V = zeros(n ,m+1); +H = zeros(m+1,m); + +beta = norm(v); +V(:,1) = v/beta; +resnorm = inf; + +j=0; + +while resnorm > tol + j = j+1; + z = A*V(:,j); + for i=1:j + H(i,j) = z'*V(:,i); + z = z - H(i,j)*V(:,i); + end + H(j+1,j) = norm(z); + e1 = zeros(j,1); e1(1) = 1; + ej = zeros(j,1); ej(j) = 1; + s = [0.01, 1/3, 2/3, 1]*t; + for q=1:length(s) + u = expm(-s(q)*H(1:j,1:j))*e1; + beta_j(q) = -H(j+1,j)* (ej'*u); + end + resnorm = norm(beta_j,'inf'); + fprintf('j = %d, resnorm = %.2e\n',j,resnorm); + if resnorm<=toler + break + elseif j==m + disp('warning: no convergence within m steps'); + end + V(:,j+1) = z/H(j+1,j); +end +y = V(:,1:j)*(beta*u);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/Magnus_4.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,29 @@ +% Takes one time step of size k using a fourth order magnus integrator +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = Magnus_4(v, D, t , k , matrixexp ,tol) + + + +if isa(D,'function_handle') + c1 = 1/2 - sqrt(3)/6; + c2 = 1/2 + sqrt(3)/6; + + A1 = D(t +c1*k); + A2 = D(t + c2*k); + Omega = 1/2*(A1 + A2) + sqrt(3)*k/12*(A1*A2-A2*A1); +else + Omega = D; +end + + +switch matrixexp + case 'expm' + v = expm(k*Omega)*v; + case 'Arnoldi' + v = time.expint.expm_Arnoldi(-Omega,v,k,tol,100); + otherwise + error('No such matrix exponential evaluation') + +end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/Magnus_mp.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,21 @@ +% Takes one time step of size k using the magnus midpoinr +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = Magnus_mp(v,D, t , k,matrixexp,tol) + +if isa(D,'function_handle') + Omega = D(t +k/2); +else + Omega = D; +end + +switch matrixexp + case 'expm' + v = expm(k*Omega)*v; + case 'Arnoldi' + v = time.expint.expm_Arnoldi(-Omega,v,k,tol,1000); + otherwise + error('No such matrix exponential evaluation') + +end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/Magnus_mp.m~ Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,21 @@ +% Takes one time step of size k using the magnus midpoinr +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = Magnus_mp(v,D, t , k,matrixexp,tol) + +if isa(D,'function_handle') + Omega = D(t +k/2); +else + Omega = D; +end + +switch matrixexp + case 'expm' + v = expm(k*Omega)*v; + case 'Arnoldi' + [v, i] = time.expint.expm_Arnoldi(-Omega,v,k,tol,1000); + otherwise + error('No such matrix exponential evaluation') + +end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/expm_Arnoldi.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,59 @@ +function y = expm_Arnoldi(A,v,t,toler,m) + +global iter + +% y = expm_Arnoldi(A,v,t,toler,m) +% +% computes $y = \exp(-t A) v$ +% input: A (n x n)-matrix, v n-vector, t>0 time interval, +% toler>0 tolerance, m maximal Krylov dimension +% +% Copyright (c) 2012 by M.A. Botchev +% Permission to copy all or part of this work is granted, +% provided that the copies are not made or distributed +% for resale, and that the copyright notice and this +% notice are retained. +% +% THIS WORK IS PROVIDED ON AN "AS IS" BASIS. THE AUTHOR +% PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED, +% REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS +% MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE. +% +n = size (v,1); +V = zeros(n ,m+1); +H = zeros(m+1,m); + +beta = norm(v); +V(:,1) = v/beta; +resnorm = inf; + +j=0; +iter = 0; + + +while resnorm > toler + iter = iter +1; + j = j+1; + w = A*V(:,j); + for i=1:j + H(i,j) = w'*V(:,i); + w = w - H(i,j)*V(:,i); + end + H(j+1,j) = norm(w); + e1 = zeros(j,1); e1(1) = 1; + ej = zeros(j,1); ej(j) = 1; + s = [0.01, 1/3, 2/3, 1]*t; + for q=1:length(s) + u = expm(-s(q)*H(1:j,1:j))*e1; + beta_j(q) = -H(j+1,j)* (ej'*u); + end + resnorm = norm(beta_j,'inf'); + % fprintf('j = %d, resnorm = %.2e\n',j,resnorm); + if resnorm<=toler + break + elseif j==m + disp('warning: no convergence within m steps'); + end + V(:,j+1) = w/H(j+1,j); +end +y = V(:,1:j)*(beta*u);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/Magnus4.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,49 @@ +classdef Magnus4 < time.Timestepper + properties + D + S + F + k + t + v + m + n + matrixexp + tol + end + + + methods + function obj = Magnus4(D, k, t0, v0,matrixexp,tol) + default_arg('matrixexp','expm') + default_arg('tol',1e-6) + obj.D = D; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.m = length(v0); + obj.n = 0; + obj.matrixexp = matrixexp; + obj.tol = tol; + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + obj.v = time.expint.Magnus_4(obj.v,obj.D, obj.t, obj.k, obj.matrixexp, obj.tol); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + k = rk4.get_rk4_time_step(lambda); + end + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/MagnusMP.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,49 @@ +classdef MagnusMP < time.Timestepper + properties + D + S + F + k + t + v + m + n + matrixexp + tol + end + + + methods + function obj = MagnusMP(D, k ,t0,v0, matrixexp,tol) + default_arg('matrixexp','expm') + default_arg('tol',1e-6) + obj.D = D; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.m = length(v0); + obj.n = 0; + obj.matrixexp = matrixexp; + obj.tol = tol; + end + + function [v,t] = getV(obj) + v = obj.v; + t = obj.t; + end + + function obj = step(obj) + obj.v = time.expint.Magnus_mp(obj.v,obj.D, obj.t, obj.k,obj.matrixexp,obj.tol); + obj.t = obj.t + obj.k; + obj.n = obj.n + 1; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + k = rk4.get_rk4_time_step(lambda); + end + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diffSymfun.m Tue Feb 20 15:00:30 2018 +0100 @@ -0,0 +1,7 @@ +% Differentiates a symbolic function like diff does, but keeps the function as a symfun +function g = diffSymfun(f, varargin) + assertType(f, 'symfun'); + + args = argnames(f); + g = symfun(diff(f,varargin{:}), args); +end