Mercurial > repos > public > sbplib
changeset 755:14f0058356f2 feature/d1_staggered
Merge with feature/grids
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Fri, 15 Jun 2018 18:10:26 -0700 |
| parents | 5264ce57b573 (current diff) 87436a107d8a (diff) |
| children | f891758ad7a4 |
| files | |
| diffstat | 11 files changed, 179 insertions(+), 53 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Fri Jun 15 14:01:13 2018 -0700 +++ b/+multiblock/DiffOp.m Fri Jun 15 18:10:26 2018 -0700 @@ -10,13 +10,13 @@ end methods - function obj = DiffOp(doHand, grid, order, doParam) + function obj = DiffOp(doHand, g, 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 + % g -- 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 @@ -26,9 +26,9 @@ % extra parameters to all doHand: doHand(..., doParam{:}) default_arg('doParam', []) - [getHand, getParam] = parseInput(doHand, grid, doParam); + [getHand, getParam] = parseInput(doHand, g, doParam); - nBlocks = grid.nBlocks(); + nBlocks = g.nBlocks(); obj.order = order; @@ -40,7 +40,7 @@ if ~iscell(p) p = {p}; end - obj.diffOps{i} = h(grid.grids{i}, order, p{:}); + obj.diffOps{i} = h(g.grids{i}, order, p{:}); end @@ -53,7 +53,7 @@ % Build the differentiation matrix - obj.blockmatrixDiv = {grid.Ns, grid.Ns}; + obj.blockmatrixDiv = {g.Ns, g.Ns}; D = blockmatrix.zero(obj.blockmatrixDiv); for i = 1:nBlocks D{i,i} = obj.diffOps{i}.D; @@ -61,7 +61,7 @@ for i = 1:nBlocks for j = 1:nBlocks - intf = grid.connections{i,j}; + intf = g.connections{i,j}; if isempty(intf) continue end @@ -77,14 +77,15 @@ end end obj.D = blockmatrix.toMatrix(D); + obj.grid = g; - function [getHand, getParam] = parseInput(doHand, grid, doParam) - if ~isa(grid, 'multiblock.Grid') + function [getHand, getParam] = parseInput(doHand, g, doParam) + if ~isa(g, 'multiblock.Grid') error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.'); end - if iscell(doHand) && length(doHand) == grid.nBlocks() + if iscell(doHand) && length(doHand) == g.nBlocks() getHand = @(i)doHand{i}; elseif isa(doHand, 'function_handle') getHand = @(i)doHand; @@ -104,7 +105,7 @@ % doParam is a non-empty cell-array - if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam)) + if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam)) % doParam is a cell-array of cell-arrays getParam = @(i)doParam{i}; return @@ -116,7 +117,7 @@ function ops = splitOp(obj, op) % Splits a matrix operator into a cell-matrix of matrix operators for - % each grid. + % each g. ops = sparse2cell(op, obj.NNN); end
--- a/+noname/calculateErrors.m Fri Jun 15 14:01:13 2018 -0700 +++ b/+noname/calculateErrors.m Fri Jun 15 18:10:26 2018 -0700 @@ -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()
--- a/+scheme/bcSetup.m Fri Jun 15 14:01:13 2018 -0700 +++ b/+scheme/bcSetup.m Fri Jun 15 18:10:26 2018 -0700 @@ -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{i},'data') || isempty(bc{i}.data) continue end assertType(bc{i}.data, 'function_handle');
--- a/Color.m Fri Jun 15 14:01:13 2018 -0700 +++ b/Color.m Fri Jun 15 18:10:26 2018 -0700 @@ -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
--- a/TextTable.m Fri Jun 15 14:01:13 2018 -0700 +++ b/TextTable.m Fri Jun 15 18:10:26 2018 -0700 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/convergencePlot.m Fri Jun 15 18:10:26 2018 -0700 @@ -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
--- a/gaussian.m Fri Jun 15 14:01:13 2018 -0700 +++ b/gaussian.m Fri Jun 15 18:10:26 2018 -0700 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/mononomial.m Fri Jun 15 18:10:26 2018 -0700 @@ -0,0 +1,8 @@ +function y = mononomial(x, k) + if k < 0 + y = x*0; + return + end + y = x.^k/factorial(k); +end +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rickerWavelet.m Fri Jun 15 18:10:26 2018 -0700 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stencilEquation.m Fri Jun 15 18:10:26 2018 -0700 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/vandermonde.m Fri Jun 15 18:10:26 2018 -0700 @@ -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