Mercurial > repos > public > sbplib

diff -r c9147e05d228 -r 348d5bcf7daf +noname/calculateErrors.m
--- a/+noname/calculateErrors.m	Tue Feb 20 14:58:26 2018 +0100
+++ b/+noname/calculateErrors.m	Tue Feb 20 15:00:30 2018 +0100
@@ -4,27 +4,40 @@
 % m are grid size parameters.
 % N are number of timesteps to use for each gird size
 % timeOpt are options for the timeStepper
+% errorFun is a function_handle taking 2 or 3 arguments, errorFun(trueSolution, approxSolution), errorFun(trueSolution, approxSolution, discr)
 function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt)
+    %TODO: Ability to choose paralell or not
     assertType(schemeFactory, 'function_handle');
     assertNumberOfArguments(schemeFactory, 1);
     assertScalar(T);
     assert(length(m) == length(N), 'Vectors m and N must have the same length');
     assertType(errorFun, 'function_handle');
-    assertNumberOfArguments(errorFun, 2);
-    default_arg('timeOpt');
+
+    if ~ismember(nargin(errorFun), [2,3])
+        error('sbplib:noname:calculateErrors:wrongNumberOfArguments', '"%s" must have 2 or 3, found %d', toString(errorFun), nargin(errorFun));
+    end

-    e = [];
-    for i = 1:length(m)
+    default_arg('timeOpt', struct());
+
+
+    e = zeros(1,length(m));
+    parfor i = 1:length(m)
         done = timeTask('m = %3d ', m(i));

         [discr, trueSolution] = schemeFactory(m(i));

-        timeOpt.k = T/N(i);
-        ts = discr.getTimestepper(timeOpt);
+        timeOptTemp = timeOpt;
+        timeOptTemp.k = T/N(i);
+        ts = discr.getTimestepper(timeOptTemp);
         ts.stepTo(N(i), true);
         approxSolution = discr.getTimeSnapshot(ts);

-        e(i) = errorFun(trueSolution, approxSolution);
+        switch nargin(errorFun)
+            case 2
+                e(i) = errorFun(trueSolution, approxSolution);
+            case 3
+                e(i) = errorFun(trueSolution, approxSolution, discr);
+        end

         fprintf('e = %.4e', e(i))
         done()
diff -r c9147e05d228 -r 348d5bcf7daf +scheme/bcSetup.m
--- a/+scheme/bcSetup.m	Tue Feb 20 14:58:26 2018 +0100
+++ b/+scheme/bcSetup.m	Tue Feb 20 15:00:30 2018 +0100
@@ -22,7 +22,7 @@
         [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
         closure = closure + localClosure;

-        if isempty(bc{i}.data)
+        if ~isfield(bc,'data') || isempty(bc{i}.data)
             continue
         end
         assertType(bc{i}.data, 'function_handle');
diff -r c9147e05d228 -r 348d5bcf7daf Color.m
--- a/Color.m	Tue Feb 20 14:58:26 2018 +0100
+++ b/Color.m	Tue Feb 20 15:00:30 2018 +0100
@@ -10,6 +10,10 @@
         black     = [0.000 0.000 0.000];
         white     = [1.000 1.000 1.000];
         colors = { Color.blue, Color.red, Color.yellow, Color.green, Color.purple, Color.lightblue, Color.darkred, Color.black, Color.white};
+        markers = {'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'};
+        lineStyles = {'-', '--', ':', '-.'};
+
+        solidMarkers = {'o', 'square', 'diamond', 'v', 'pentagram', '^', '>', '<', 'hexagram'};

         notabilityYellow     = [100.0   99.0    22.0    ]/100;
         notabilityOrange     = [97.0    61.0    15.0    ]/100;
@@ -34,13 +38,11 @@

     methods(Static)
         function sample()
-            markers ={'+', 'o', '*', '.', 'x', 'square', 'diamond', 'v', '^', '>', '<', 'pentagram', 'hexagram'};
             % Filled and non-filled markers?
-            lineStyles = {'-', '--', ':', '-.'};


             function showMarkers(x0, y0, lx, ly, color, filled)
-                n = length(markers);
+                n = length(Color.markers);
                 s = ceil(sqrt(n));

                 x = linspace(x0, x0 + lx, s);
@@ -50,7 +52,7 @@

                 for i = 1:n
                     lh = line(X(i),Y(i));
-                    lh.Marker = markers{i};
+                    lh.Marker = Color.markers{i};
                     lh.MarkerSize = 12;
                     lh.Color = color;

@@ -79,13 +81,13 @@
             end

             function showLines(y0, ly, A, w)
-                n = length(lineStyles);
+                n = length(Color.lineStyles);
                 x = linspace(0,1,100);
                 y = linspace(y0, y0+ ly, n);
                 for i = 1:n
                     lh = line(x, y(i) + A*sin(pi*x*w));
                     lh.LineWidth = 2;
-                    lh.LineStyle = lineStyles{i};
+                    lh.LineStyle = Color.lineStyles{i};
                 end
             end
diff -r c9147e05d228 -r 348d5bcf7daf TextTable.m
--- a/TextTable.m	Tue Feb 20 14:58:26 2018 +0100
+++ b/TextTable.m	Tue Feb 20 15:00:30 2018 +0100
@@ -4,28 +4,36 @@
         fmtArray
         vertDiv
         horzDiv
-
-        nCols
-        nRows
     end

     methods
-        function obj = TextTable(data, vertDiv, horzDiv);
+        function obj = TextTable(data, vertDiv, horzDiv)
             default_arg('vertDiv', []);
             default_arg('horzDiv', []);

-
             obj.data = data;
             obj.vertDiv = vertDiv;
             obj.horzDiv = horzDiv;

-            [obj.nRows, obj.nCols] = size(data);
             obj.fmtArray = cell(size(data));
             obj.formatAll('%s');

         end

+        function n = nRows(obj)
+            n = size(obj.data, 1);
+        end
+
+        function n = nCols(obj)
+            n = size(obj.data, 2);
+        end
+
+        function print(obj)
+            disp(obj.toString());
+        end
+
         function formatAll(obj, fmt)
+            obj.fmtArray = cell(size(obj.data));
             obj.fmtArray(:,:) = {fmt};
         end

@@ -58,28 +66,31 @@

             str = '';

+            N = size(strArray, 2);
+
             % First horzDiv
-            if ismember(0, obj.horzDiv)
+            if isDiv(0, obj.horzDiv, N);
                 str = [str, obj.getHorzDiv(widths)];
             end

             for i = 1:obj.nRows
-                str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv, obj.horzDiv)];
+                str = [str, TextTable.rowToString(strArray(i,:), widths, obj.vertDiv)];

                 % Interior horzDiv
-                if ismember(i, obj.horzDiv)
+                if isDiv(i, obj.horzDiv, N)
                     str = [str, obj.getHorzDiv(widths)];
                 end
             end
         end

         function str = getHorzDiv(obj, widths)
-            str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv, obj.horzDiv);
+            str = TextTable.rowToString(cell(1,obj.nCols), widths, obj.vertDiv);
             str(find(' ' == str)) = '-';
             str(find('|' == str)) = '+';
         end

         function strArray = getStringArray(obj)
+            assert(all(size(obj.data) == size(obj.fmtArray)), 'Sizes of format matrix and data matrix do not match.')
             strArray = cell(size(obj.data));

             for i = 1:obj.nRows
@@ -91,32 +102,42 @@
     end

     methods (Static)
-        function str = rowToString(strs, widths, vertDiv, horzDiv)
+        function str = rowToString(strs, widths, vertDiv)
+            N = length(strs);
+
             % First vertDiv
-            if ismember(0, vertDiv)
-                str = '| ';
+            if isDiv(0, vertDiv, N)
+                prefix = '| ';
             else
-                str = ' ';
+                prefix = ' ';
             end

-            % Interior cols
-            for j = 1:length(strs) - 1
-                str = [str, sprintf('%*s ', widths(j), strs{j})];
+            % Pad strings
+            for i = 1:N
+                strs{i} = sprintf('%*s', widths(i), strs{i});
+            end

-                % Interior vertDiv
-                if ismember(j, vertDiv)
-                    str = [str, '| '];
+            % Column delimiters
+            delims = cell(1,N-1);
+            for i = 1:length(delims)
+                if isDiv(i, vertDiv, N);
+                    delims{i} = '| ';
+                else
+                    delims{i} = ' ';
                 end
             end

-            % Last col
-            str = [str, sprintf('%*s ', widths(end), strs{end})];
-
-            if ismember(length(strs), vertDiv)
-                str = [str, '|'];
+            if isDiv(N, vertDiv, N);
+                suffix = '|';
+            else
+                suffix = '';
             end

-            str = [str, sprintf('\n')];
+            str = [prefix, strjoin(strs, delims), suffix, sprintf('\n')];
         end
     end
+end
+
+function b = isDiv(i, div, N)
+    b = ismember(i, div) || ismember(i, N+div+1);
 end
\ No newline at end of file
diff -r c9147e05d228 -r 348d5bcf7daf convergencePlot.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/convergencePlot.m	Tue Feb 20 15:00:30 2018 +0100
@@ -0,0 +1,55 @@
+function hand = convergencePlot(orders, h, e)
+    N = length(orders);
+
+    fh = figure();
+    ah = axes();
+    ah.XScale = 'log';
+    ah.YScale = 'log';
+    hold on
+    ph = {};
+    phc = {};
+    legends = {};
+    for i = 1:N
+        ph{i} = loglog(h{i}, e{i});
+        phc{i} = plotConvergenceFit(orders{i}, h{i}, e{i});
+
+        ph{i}.LineStyle = 'none';
+        ph{i}.Marker = Color.solidMarkers{i};
+        ph{i}.MarkerSize = 12;
+        ph{i}.Color = Color.colors{i};
+        ph{i}.MarkerFaceColor = Color.colors{i};
+
+        legends{i} = sprintf('$o = %d$', orders{i});
+    end
+    hold off
+
+    lh = legend([ph{:}], legends);
+    lh.Interpreter = 'latex';
+    lh.Location = 'SouthEast';
+
+    for i = 1:N
+        uistack(phc{i}, 'bottom');
+    end
+
+    xlabel('$h$', 'interpreter', 'latex')
+    ylabel('Error', 'interpreter', 'latex')
+
+    % xlim([0.7e-2, 1e-1])
+    % ylim([3e-5, 4])
+
+    grid on
+
+    ah = gca();
+    ah.TickLabelInterpreter = 'latex';
+    setFontSize(fh);
+
+    % if savePngs
+    %     savepng(fh, 'fig/conv/conv',600)
+    % end
+
+    hand = struct();
+    hand.fig = fh;
+    hand.data = ph;
+    hand.fits = phc;
+    hand.legend = lh;
+end
diff -r c9147e05d228 -r 348d5bcf7daf gaussian.m
--- a/gaussian.m	Tue Feb 20 14:58:26 2018 +0100
+++ b/gaussian.m	Tue Feb 20 15:00:30 2018 +0100
@@ -1,3 +1,3 @@
 function z = gaussian(x,x0,d)
-    z = exp(-norm(x-x0).^2/d^2);
+    z = exp(-sum((x-x0).^2,2)/d^2);
 end
\ No newline at end of file
diff -r c9147e05d228 -r 348d5bcf7daf mononomial.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/mononomial.m	Tue Feb 20 15:00:30 2018 +0100
@@ -0,0 +1,8 @@
+function y = mononomial(x, k)
+    if k < 0
+        y = x*0;
+        return
+    end
+    y = x.^k/factorial(k);
+end
+
diff -r c9147e05d228 -r 348d5bcf7daf rickerWavelet.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rickerWavelet.m	Tue Feb 20 15:00:30 2018 +0100
@@ -0,0 +1,3 @@
+function y = rickerWavelet(x, x0, A)
+    y = (1-2*pi^2*A^2*(x-x0).^2).*exp(-pi^2*A^2*(x-x0).^2);
+end
diff -r c9147e05d228 -r 348d5bcf7daf stencilEquation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/stencilEquation.m	Tue Feb 20 15:00:30 2018 +0100
@@ -0,0 +1,13 @@
+% Find the equation for the stencil for d^k/dx^k
+function [A,b] = stencilEquation(k, offsets, order)
+    q = sym('q', [1, length(offsets)]);
+
+    p = 0:(order-1+k);
+
+    v     = vandermonde(offsets, p);
+    vdiff = vandermonde(      0, p-k);
+
+    eq = q*v == vdiff;
+
+    [A,b] = equationsToMatrix(eq, q);
+end
diff -r c9147e05d228 -r 348d5bcf7daf vandermonde.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/vandermonde.m	Tue Feb 20 15:00:30 2018 +0100
@@ -0,0 +1,10 @@
+% Create vandermonde matrix for points x and polynomials of order p
+% x and p are vectors
+% v is a length(x) by length(p) matrix
+function V = vandermonde(x, p)
+    V = sym(zeros(length(x), length(p))); % Is there a way to make this work for both double and sym
+
+    for i = 1:length(p)
+        V(:, i) = mononomial(x,p(i));
+    end
+end