Mercurial > repos > public > sbplib

changeset 431:5f4540e13f9b feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
--- /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
--- /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
--- 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
-
-
-
-
-
--- 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
 
 
--- 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
--- /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
--- /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
--- /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
--- /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