Mercurial > repos > public > sbplib
changeset 431:5f4540e13f9b feature/quantumTriangles
Meged with default
| author | Ylva Rydin <ylva.rydin@telia.com> |
|---|---|
| date | Wed, 08 Feb 2017 09:18:08 +0100 |
| parents | 25053554524b (current diff) b6a5dc423990 (diff) |
| children | eca4ca84cf0a |
| files | +sbp/D2Standard.m |
| diffstat | 9 files changed, 437 insertions(+), 100 deletions(-) [+] |
line wrap: on
line diff
diff -r 25053554524b -r 5f4540e13f9b +sbp/+implementations/d1_gauss_4.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_gauss_4.m Wed Feb 08 09:18:08 2017 +0100 @@ -0,0 +1,63 @@ +function [D1,H,x,h,e_l,e_r] = d1_gauss_4(L) + +% L: Domain length +% N: Number of grid points +if(nargin < 2) + L = 1; +end + +N = 4; + +% Quadrature nodes on interval [-1, 1] +x = [ -0.8611363115940526; -0.3399810435848563; 0.3399810435848563; 0.8611363115940526]; + +% Shift nodes to [0,L] +x = (x+1)/2*L; + +% Boundary extrapolation operators +e_l = [1.5267881254572668; -0.8136324494869273; 0.4007615203116504; -0.1139171962819899]; +e_r = flipud(e_l); +e_l = sparse(e_l); +e_r = sparse(e_r); + +%%%% Compute approximate h %%%%%%%%%% +h = L/(N-1); +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Norm matrix on [-1,1] %%%%%%%% +P = sparse(N,N); +P(1,1) = 0.3478548451374539; +P(2,2) = 0.6521451548625461; +P(3,3) = 0.6521451548625461; +P(4,4) = 0.3478548451374539; +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% Norm matrix on [0,L] %%%%%%%% +H = P*L/2; +%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%% D1 on [-1,1] %%%%%%%% +D1 = sparse(N,N); +D1(1,1) = -3.3320002363522817; +D1(1,2) = 4.8601544156851962; +D1(1,3) = -2.1087823484951789; +D1(1,4) = 0.5806281691622644; + +D1(2,1) = -0.7575576147992339; +D1(2,2) = -0.3844143922232086; +D1(2,3) = 1.4706702312807167; +D1(2,4) = -0.3286982242582743; + +D1(3,1) = 0.3286982242582743; +D1(3,2) = -1.4706702312807167; +D1(3,3) = 0.3844143922232086; +D1(3,4) = 0.7575576147992339; + +D1(4,1) = -0.5806281691622644; +D1(4,2) = 2.1087823484951789; +D1(4,3) = -4.8601544156851962; +D1(4,4) = 3.3320002363522817; +%%%%%%%%%%%%%%%%%%%%%%%%% + +% D1 on [0,L] +D1 = D1*2/L; \ No newline at end of file
diff -r 25053554524b -r 5f4540e13f9b +sbp/D1Gauss.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D1Gauss.m Wed Feb 08 09:18:08 2017 +0100 @@ -0,0 +1,41 @@ +classdef D1Gauss < sbp.OpSet + % Diagonal-norm SBP operators based on the Gauss quadrature formula + % with m nodes, which is of degree 2m-1. Hence, The operator D1 is + % accurate of order m. + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D1Gauss(m,lim) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + + switch m + case 4 + [obj.D1,obj.H,obj.x,obj.h,obj.e_l,obj.e_r] = ... + sbp.implementations.d1_gauss_4(L); + otherwise + error('Invalid number of points: %d.', m); + end + + + obj.x = obj.x + x_l; + obj.HI = inv(obj.H); + obj.Q = obj.H*obj.D1 - obj.e_r*obj.e_r' + obj.e_l*obj.e_l'; + + obj.borrowing = []; + end + end +end
diff -r 25053554524b -r 5f4540e13f9b +sbp/D1Nonequidistant.m --- a/+sbp/D1Nonequidistant.m Wed Feb 08 09:07:07 2017 +0100 +++ b/+sbp/D1Nonequidistant.m Wed Feb 08 09:18:08 2017 +0100 @@ -11,23 +11,23 @@ x % grid borrowing % Struct with borrowing limits for different norm matrices end - + methods function obj = D1Nonequidistant(m,lim,order,option) - + default_arg('option','Accurate'); % 'Accurate' operators are optimized for accuracy % 'Minimal' operators have the smallest possible boundary % closure - + x_l = lim{1}; x_r = lim{2}; L = x_r-x_l; - + switch option - + case {'Accurate','accurate','A'} - + if order == 4 [obj.D1,obj.H,obj.x,obj.h] = ... sbp.implementations.d1_noneq_4(m,L); @@ -46,9 +46,9 @@ else error('Invalid operator order %d.',order); end - + case {'Minimal','minimal','M'} - + if order == 4 [obj.D1,obj.H,obj.x,obj.h] = ... sbp.implementations.d1_noneq_minimal_4(m,L); @@ -67,28 +67,20 @@ else error('Invalid operator order %d.',order); end - + end - + obj.x = obj.x + x_l; - + obj.e_l = sparse(m,1); obj.e_r = sparse(m,1); obj.e_l(1) = 1; obj.e_r(m) = 1; - + obj.HI = inv(obj.H); obj.Q = obj.H*obj.D1 - obj.e_r*obj.e_r' + obj.e_l*obj.e_l'; - + obj.borrowing = []; - end end - - end - - - - -
diff -r 25053554524b -r 5f4540e13f9b +sbp/D2Standard.m --- a/+sbp/D2Standard.m Wed Feb 08 09:07:07 2017 +0100 +++ b/+sbp/D2Standard.m Wed Feb 08 09:18:08 2017 +0100 @@ -14,7 +14,7 @@ h % Step size x % grid borrowing % Struct with borrowing limits for different norm matrices - + end methods @@ -63,15 +63,12 @@ end obj.m = m; - end function str = string(obj) str = [class(obj) '_' num2str(obj.order)]; end end - - end
diff -r 25053554524b -r 5f4540e13f9b +time/SBPInTime.m --- a/+time/SBPInTime.m Wed Feb 08 09:07:07 2017 +0100 +++ b/+time/SBPInTime.m Wed Feb 08 09:18:08 2017 +0100 @@ -1,89 +1,92 @@ classdef SBPInTime < time.Timestepper % The SBP in time method. % Implemented for v_t = A*v + f(t) - % k_local -- time-step - % Nblock -- number of points in each block - % nodes -- points such that t_n + nodes are the points in block n. - % Each "step" takes one block step and thus advances - % k = k_local*(Nblock-1) in time. - % M -- matrix used in every solve. - % [L,U,P,Q] = lu(M); + % + % Each "step" takes one block step and thus advances + % k = k_local*(blockSize-1) in time. properties - M - L - U - P - Q + M % System matrix + L,U,P,Q % LU factorization of M A Et_r penalty f - k_local - k + k_local % step size within a block + k % Time size of a block k/(blockSize-1) = k_local t v m n - Nblock + blockSize % number of points in each block order nodes end methods - function obj = SBPInTime(A, f, k, order, Nblock, t0, v0, TYPE) - default_arg('TYPE','equidistant'); - default_arg('Nblock',time.SBPInTime.smallestBlockSize(order,TYPE)); - + function obj = SBPInTime(A, f, k, t0, v0, TYPE, order, blockSize) + + default_arg('TYPE','gauss'); + + if(strcmp(TYPE,'gauss')) + default_arg('order',4) + default_arg('blockSize',4) + else + default_arg('order', 8); + default_arg('blockSize',time.SBPInTime.smallestBlockSize(order,TYPE)); + end + obj.A = A; obj.f = f; - obj.k_local = k; - obj.k = k*(Nblock-1); - obj.Nblock = Nblock; + obj.k_local = k/(blockSize-1); + obj.k = k; + obj.blockSize = blockSize; obj.t = t0; obj.m = length(v0); obj.n = 0; - + %==== Build the time discretization matrix =====% switch TYPE case 'equidistant' - ops = sbp.D2Standard(Nblock,{0,obj.k},order); + ops = sbp.D2Standard(blockSize,{0,obj.k},order); case 'optimal' - ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order); + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order); case 'minimal' - ops = sbp.D1Nonequidistant(Nblock,{0,obj.k},order,'minimal'); + ops = sbp.D1Nonequidistant(blockSize,{0,obj.k},order,'minimal'); + case 'gauss' + ops = sbp.D1Gauss(blockSize,{0,obj.k}); end - + D1 = ops.D1; HI = ops.HI; e_l = ops.e_l; e_r = ops.e_r; obj.nodes = ops.x; - + Ix = speye(size(A)); - It = speye(Nblock,Nblock); - - obj.Et_r = kron(e_r,Ix); - + It = speye(blockSize,blockSize); + + obj.Et_r = kron(e_r,Ix); + % Time derivative + penalty tau = 1; - Mt = D1 + tau*HI*(e_l*e_l'); - + Mt = D1 + tau*HI*(e_l*e_l'); + % penalty to impose "data" penalty = tau*HI*e_l; obj.penalty = kron(penalty,Ix); - + Mx = kron(It,A); - Mt = kron(Mt,Ix); + Mt = kron(Mt,Ix); obj.M = Mt - Mx; %==============================================% - + % LU factorization [obj.L,obj.U,obj.P,obj.Q] = lu(obj.M); - + % Pretend that the initial condition is the last level % of a previous step. - obj.v = obj.Et_r * v0; - + obj.v = 1/(e_r'*e_r) * obj.Et_r * v0; + end function [v,t] = getV(obj) @@ -93,39 +96,20 @@ function obj = step(obj) obj.v = time.sbp.sbpintime(obj.v, obj.t, obj.nodes,... - obj.penalty, obj.f, obj.Nblock,... + obj.penalty, obj.f, obj.blockSize,... obj.Et_r,... obj.L, obj.U, obj.P, obj.Q); obj.t = obj.t + obj.k; - obj.n = obj.n + obj.Nblock-1; + obj.n = obj.n + 1; end end - - + methods(Static) - - % - function [k,numberOfBlocks] = alignedTimeStep(k,Tend,Nblock) - - % input k is the desired time-step - % Nblock is the number of points per block. - - % Make sure that we reach the final time by advancing - % an integer number of blocks - kblock = (Nblock-1)*k; - numberOfBlocks = ceil(Tend/kblock); - kblock = Tend/(numberOfBlocks); + function N = smallestBlockSize(order,TYPE) + default_arg('TYPE','gauss') - % Corrected time step - k = kblock/(Nblock-1); - - end - - function N = smallestBlockSize(order,TYPE) - default_arg('TYPE','equidistant') - switch TYPE - + case 'equidistant' switch order case 2 @@ -143,9 +127,9 @@ otherwise error('Operator does not exist'); end - + case 'optimal' - + switch order case 4 N = 8; @@ -160,9 +144,9 @@ otherwise error('Operator does not exist'); end - + case 'minimal' - + switch order case 4 N = 6; @@ -177,13 +161,9 @@ otherwise error('Operator does not exist'); end - + case 'gauss' + N = 4; end - end - end - - - end
diff -r 25053554524b -r 5f4540e13f9b +time/SBPInTimeSecondOrderForm.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/SBPInTimeSecondOrderForm.m Wed Feb 08 09:18:08 2017 +0100 @@ -0,0 +1,67 @@ +classdef SBPInTimeSecondOrderForm < time.Timestepper + properties + A,B,C + M, f + + n + t + k + + firstOrderTimeStepper + end + + methods + % Solves u_tt = Au + Bu_t + C + % A, B can either both be constants or both be function handles, + % They can also be omitted by setting them equal to the empty matrix. + function obj = SBPInTimeSecondOrderForm(A, B, C, k, t0, v0, v0t, TYPE, order, blockSize) + default_arg('TYPE', []); + default_arg('order', []); + default_arg('blockSize',[]); + + m = length(v0); + + default_arg('A', sparse(m, m)); + default_arg('B', sparse(m, m)); + default_arg('C', sparse(m, 1)); + + I = speye(m); + O = sparse(m,m); + + obj.M = [ + O, I; + A, B; + ]; + obj.f = @(t)[ + sparse(m,1); + C; + ]; + + w0 = [v0; v0t]; + + obj.k = k; + obj.t = t0; + obj.n = 0; + + obj.firstOrderTimeStepper = time.SBPInTime(obj.M, obj.f, obj.k, obj.t, w0, TYPE, order, blockSize); + end + + function [v,t] = getV(obj) + w = obj.firstOrderTimeStepper.getV(); + v = w(1:end/2); + t = obj.t; + end + + function [vt,t] = getVt(obj) + w = obj.firstOrderTimeStepper.getV(); + vt = w(end/2+1:end); + t = obj.t; + end + + function obj = step(obj) + obj.firstOrderTimeStepper.step(); + obj.t = obj.firstOrderTimeStepper.t; + obj.n = obj.firstOrderTimeStepper.n; + end + end +end
diff -r 25053554524b -r 5f4540e13f9b Map.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/Map.m Wed Feb 08 09:18:08 2017 +0100 @@ -0,0 +1,91 @@ +classdef Map < handle + properties + map + end + + % can we support multi map using varargin? + % probably a bad idea. For example it complicates keys(); + + methods + function obj = Map() + obj.map = containers.Map(); + end + + function set(obj, k, v) + keyByteStream = getByteStreamFromArray(k); + + obj.map(char(keyByteStream)) = v; + end + + function v = get(obj, k) + keyByteStream = getByteStreamFromArray(k); + + v = obj.map(char(keyByteStream)); + end + + function b = isKey(obj, k) + keyByteStream = getByteStreamFromArray(k); + b = obj.map.isKey(char(keyByteStream)); + end + + function c = keys(obj) + keyByteStreams = obj.map.keys; + + n = length(keyByteStreams); + + c = cell(1, n); + for i = 1:n + c{i} = getArrayFromByteStream(uint8(keyByteStreams{i})); + end + end + + function l = length(obj) + l = obj.map.length; + end + + function remove(obj, k) + keyByteStream = getByteStreamFromArray(k); + obj.map.remove(char(keyByteStream)); + end + + function s = size(obj) + s = obj.map.size; + end + + function c = values(obj) + c = obj.map.values; + end + + function v = subsref(obj, S) + switch S(1).type + case '()' + k = S.subs{1}; + try + v = get(obj, k); + catch ME + if strcmp(ME.identifier,'MATLAB:Containers:Map:NoKey') + error('Reference to non-existent entry %s',toString(S.subs)); + else + throw(ME); + end + end + otherwise + try + v = builtin('subsref', obj, S); + catch e + error('You can''t use dot notation for this because Matlab(TM). What is this piece of shit software anyway?') + end + end + end + + function obj = subsasgn(obj, S, v); + switch S(1).type + case '()' + k = S.subs{1}; + set(obj, k, v); + otherwise + error('You can''t use dot notation because Matlab(TM). What is this piece of shit software anyway?') + end + end + end +end
diff -r 25053554524b -r 5f4540e13f9b MapTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MapTest.m Wed Feb 08 09:18:08 2017 +0100 @@ -0,0 +1,97 @@ +function tests = MatTest() + tests = functiontests(localfunctions); +end + +function kvp = getKeyValuePairs() + kvp = { + {1,3},1; + struct(), [1; 3; 4]; + [1,2; 4 3], struct(); + 'Hej', struct('lol', 6); + 0, 'Nej'; + }; +end + +function testSetAndGet(testCase) + keyValuePairs = getKeyValuePairs(); + + map = Map(); + + % Insert + for i = 1:length(keyValuePairs) + map(keyValuePairs{i,1}) = keyValuePairs{i,2}; + end + + % Validate output + for i = 1:length(keyValuePairs) + v = map(keyValuePairs{i,1}); + testCase.verifyEqual(v, keyValuePairs{i,2}); + end +end + +function map = exampleMap() + keyValuePairs = getKeyValuePairs(); + + map = Map(); + + % Insert + for i = 1:length(keyValuePairs) + map(keyValuePairs{i,1}) = keyValuePairs{i,2}; + end +end + +function testLength(testCase) + map = Map(); + testCase.verifyEqual(map.length, 0); + + map = exampleMap(); + testCase.verifyEqual(map.length, 5) +end + + +function testIsKey(testCase) + map = exampleMap(); + + keyValuePairs = getKeyValuePairs(); + keys = keyValuePairs(:,1); + + for i = 1:length(keys) + testCase.verifyTrue(map.isKey(keys{i})); + end + + notKeys = { + 'hej', + [], + 1, + {2,5}, + }; + + for i = 1:length(notKeys) + testCase.verifyFalse(map.isKey(notKeys{i})); + end +end + + +function testRemove(testCase) + map = exampleMap(); + + remove(map, struct()); + + testCase.verifyFalse(map.isKey(struct())); +end + +% function testValues(testCase) +% keyValuePairs = getKeyValuePairs(); + +% map = exampleMap(); + +% testCase.verifyEqual(values(map), keyValuePairs(:,2)'); +% end + +% function testKeys(testCase) +% keyValuePairs = getKeyValuePairs(); + +% map = exampleMap(); + +% testCase.verifyEqual(keys(map), keyValuePairs(:,1)'); +% end
diff -r 25053554524b -r 5f4540e13f9b diags.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diags.m Wed Feb 08 09:18:08 2017 +0100 @@ -0,0 +1,9 @@ +function A = diags(B,d,m,n) + assert(size(B,1) == m); + + A = repmat(B(:,1)*0, [1, n]); + + for i = 1:size(B,2) + A(:,d(i)+ (1:m)) = A(:,d(i)+ (1:m)) + diag(B(:,i)); + end +end