Mercurial > repos > public > sbplib
changeset 1300:196123459178
Merge in feature/boundary_optimized_grids
| author | Vidar Stiernström <vidar.stiernstrom@it.uu.se> |
|---|---|
| date | Wed, 08 Jul 2020 18:22:54 +0200 |
| parents | 8ec777fb473e (current diff) 73e52c74baac (diff) |
| children | 8978521b0f06 |
| files | |
| diffstat | 16 files changed, 372 insertions(+), 524 deletions(-) [+] |
line wrap: on
line diff
diff -r 8ec777fb473e -r 196123459178 +grid/boundaryOptimized.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryOptimized.m Wed Jul 08 18:22:54 2020 +0200 @@ -0,0 +1,54 @@ +% Creates a Cartesian grid of dimension length(m) +% over the domain xlim, ylim, ... +% The grid is non-equidistant in the boundary regions, +% with node placement based on boundary-optimized SBP operators. +% Examples: +% g = grid.boundaryOptimized([mx, my], xlim, ylim, order, opt) +% g = grid.boundaryOptimized([10, 15], {0,1}, {0,2}, 4) - defaults to 'accurate' stencils +% g = grid.boundaryOptimized([10, 15], {0,1}, {0,2}, 4, 'minimal') +function g = boundaryOptimized(m, varargin) + n = length(m); + + % Check that parameters matches dimensions + matchingParams = false; + if length(varargin) == n+1 % Minimal number of arguments + matchingParams = iscell([varargin{1:n}]) && ... + isfloat([varargin{n+1}]); + elseif length(varargin) == n+2 % Stencil options supplied + matchingParams = iscell([varargin{1:n}]) && ... + isfloat([varargin{n+1}]) && ... + ischar([varargin{n+2}]); + end + assert(matchingParams,'grid:boundaryOptimized:NonMatchingParameters','The number of parameters per dimensions do not match.'); + + % Check that stencil options are passed correctly (if supplied) + if length(varargin) == n+2 % Stencil options supplied + availabe_opts = ["Accurate","accurate","A","Minimal","minimal","M"]; + assert(any(varargin{n+2} == availabe_opts), ... + 'grid:boundaryOptimized:InvalidOption',"The operator option must be 'accurate' or 'minimal.'"); + else %If not passed, populate varargin with default option 'accurate' + varargin(n+2) = {'accurate'}; + end + + % Specify generating function + switch varargin{n+2} + case {'Accurate','accurate','A'} + gridgenerator = @sbp.grid.accurateBoundaryOptimizedGrid; + case {'Minimal','minimal','M'} + gridgenerator = @sbp.grid.minimalBoundaryOptimizedGrid; + end + + X = {}; + h = []; + for i = 1:n + try + [X{i},h(i)] = gridgenerator(varargin{i},m(i),varargin{n+1}); + catch exception % Propagate any errors in the grid generation functions. + msgText = getReport(exception); + error('grid:boundaryOptimized:InvalidParameter',msgText) + end + end + + g = grid.Cartesian(X{:}); + g.h = h; +end \ No newline at end of file
diff -r 8ec777fb473e -r 196123459178 +grid/boundaryOptimizedTest.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/boundaryOptimizedTest.m Wed Jul 08 18:22:54 2020 +0200 @@ -0,0 +1,111 @@ +function tests = boundaryOptimizedTest() + tests = functiontests(localfunctions); +end + +function testErrorInvalidParam(testCase) + in = { + %Invalid order + {[10 10],{0,1},{0,2},3}, + %Invalid grid size + {5, {0,1}, 4}, + {[10 5],{0,1},{0,2},4}, + {[10 5],{0,1},{0,2},6,'M'}, + %Invalid limits + {10,{1},4}, + {[10,10],{0,1},{1},4}, + {[10,10],{1},{1,0},4}, + {10,{1,0},4}, + {[10, 5],{1,0},{0,-1},4}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:InvalidParameter',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorInvalidOption(testCase) + in = { + {[8 8],{0,1},{0,2},4,'acrurate'}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:InvalidOption',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testErrorNonMatchingParam(testCase) + in = { + {[],{1},4}, + {[],{0,1},{0,1},4}, + {[5,5],{0,1},{0,1},{0,1},4}, + {[5,5,4],{0,1},{0,1},4,'accurate'} + {[5,5,4],{0,1},{0,1},{0,1},4,4}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.boundaryOptimized(in{i}{:}),'grid:boundaryOptimized:NonMatchingParameters',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +% Tests that the expected grid points are obtained for a boundary optimized grid with a 4th order +% accurate stencil and 8th order minimal stencil. +% The boundary grid point distance weights are taken from the D1Nonequidistant operators and +% grid spacing is calculated according to Mattsson et al 2018. +function testCompiles(testCase) + + %% 1D 4th order accurate stencil + % Boundary weights, number of non-equidistantly spaced points for 4th order accurate stencil + bw = [0.0000000000000e+00 6.8764546205559e-01 1.8022115125776e+00]; + n = length(bw)-1; + xi_n = bw(end); + + % Grid points in x-direction. + Lx = 1; + mx = 8; + hx_4 = Lx/(2*xi_n+(mx-2*n-1)); + + bp_l = hx_4*bw; + bp_r = Lx-flip(hx_4*bw); + interior = [hx_4*(xi_n+1) hx_4*(xi_n+2)]; + x_4 = [bp_l interior bp_r]; + + % Boundary weights, number of non-equidistantly spaced points for 8th order minimal stencil + bw = [0.0000000000000e+00, 4.9439570885261e-01, 1.4051531374839e+00]; + n = length(bw)-1; + xi_n = bw(end); + + %% 2D 8th order minimal stencil + % Grid points in x-direction. + hx_8 = Lx/(2*xi_n+(mx-2*n-1)); + + bp_l = hx_8*bw; + bp_r = Lx-flip(hx_8*bw); + interior = [hx_8*(xi_n+1) hx_8*(xi_n+2)]; + x_8 = [bp_l interior bp_r]; + + % Grid points in y-direction. + Ly = 2; + my = 9; + hy = Ly/(2*xi_n+(my-2*n-1)); + + bp_l = hy*bw; + bp_r = Ly-flip(hy*bw); + interior = [hy*(xi_n+1) hy*(xi_n+2) hy*(xi_n+3)]; + y = [bp_l interior bp_r]; + + in = { + {mx, {0,Lx},4}, + {[mx, my],{0,Lx},{0,Ly},8,'M'}, + }; + + out = { + {[x_4'],hx_4} + {[kr(x_8',ones(size(y'))),kr(ones(size(x_8')),y')],[hx_8, hy]} + }; + + for i = 1:length(in) + g = grid.boundaryOptimized(in{i}{:}); + testCase.verifyEqual(g.points(),out{i}{1},'AbsTol', 1e-14, 'RelTol', 1e-14); + testCase.verifyEqual(g.scaling(),out{i}{2},'AbsTol', 1e-14, 'RelTol', 1e-14); + end +end \ No newline at end of file
diff -r 8ec777fb473e -r 196123459178 +multiblock/+domain/Rectangle.m --- a/+multiblock/+domain/Rectangle.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+multiblock/+domain/Rectangle.m Wed Jul 08 18:22:54 2020 +0200 @@ -93,7 +93,7 @@ N_id = flat_index(m,j,1); S{j} = {S_id,'s'}; N{j} = {N_id,'n'}; - end + end boundaryGroups.E = multiblock.BoundaryGroup(E); boundaryGroups.W = multiblock.BoundaryGroup(W); boundaryGroups.S = multiblock.BoundaryGroup(S); @@ -117,6 +117,20 @@ % Returns a multiblock.Grid given some parameters % ms: cell array of [mx, my] vectors % For same [mx, my] in every block, just input one vector. + % Currently defaults to equidistant grid if varargin is empty. + % If varargin is non-empty, the first argument should supply the grid type, followed by + % additional arguments required to construct the grid. + % Grid types: + % 'equidist' - equidistant grid + % Additional argumets: none + % 'boundaryopt' - boundary optimized grid based on boundary + % optimized SBP operators + % Additional arguments: order, stencil option + % Example: g = getGrid() - the local blocks are 21x21 equidistant grids. + % g = getGrid(ms,) - block i is an equidistant grid with size given by ms{i}. + % g = getGrid(ms,'equidist') - block i is an equidistant grid with size given by ms{i}. + % g = getGrid(ms,'boundaryopt',4,'minimal') - block i is a Cartesian grid with size given by ms{i} + % and nodes placed according to the boundary optimized minimal 4th order SBP operator. function g = getGrid(obj, ms, varargin) default_arg('ms',[21,21]) @@ -129,10 +143,19 @@ ms{i} = m; end end - + if isempty(varargin) || strcmp(varargin{1},'equidist') + gridgenerator = @(m,xlim,ylim)grid.equidistant(m,xlim,ylim); + elseif strcmp(varargin{1},'boundaryopt') + order = varargin{2}; + stenciloption = varargin{3}; + gridgenerator = @(m,xlim,ylim)grid.boundaryOptimized(m,xlim,ylim,... + order,stenciloption); + else + error('No grid type supplied!'); + end grids = cell(1, obj.nBlocks); for i = 1:obj.nBlocks - grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); + grids{i} = gridgenerator(ms{i}, obj.xlims{i}, obj.ylims{i}); end g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups);
diff -r 8ec777fb473e -r 196123459178 +sbp/+grid/accurateBoundaryOptimizedGrid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+grid/accurateBoundaryOptimizedGrid.m Wed Jul 08 18:22:54 2020 +0200 @@ -0,0 +1,67 @@ +% Computes the grid points x and grid spacing h used by the boundary optimized SBP operators +% with improved boundary accuracy, presented in +% 'Boundary optimized diagonal-norm SBP operators - Mattsson, Almquist, van der Weide 2018'. +% +% lim - cell array with domain limits +% N - Number of grid points +% order - order of accuracy of sbp operator. +function [x,h] = accurateBoundaryOptimizedGrid(lim,N,order) + assert(iscell(lim) && numel(lim) == 2,'The limit should be cell array with 2 elements.'); + L = lim{2} - lim{1}; + assert(L>0,'Limits must be given in increasing order.'); + %%%% Non-equidistant grid points %%%%% + xb = boundaryPoints(order); + m = length(xb)-1; % Number of non-equidistant points + assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Compute h %%%%%%%%%% + h = L/(2*xb(end) + N-1-2*m); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Define grid %%%%%%%% + x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; + x = x + lim{1}; + %%%%%%%%%%%%%%%%%%%%%%%%% +end +function xb = boundaryPoints(order) + switch order + case 4 + x0 = 0.0000000000000e+00; + x1 = 6.8764546205559e-01; + x2 = 1.8022115125776e+00; + xb = [x0 x1 x2]'; + case 6 + x0 = 0.0000000000000e+00; + x1 = 4.4090263368623e-01; + x2 = 1.2855984345073e+00; + x3 = 2.2638953951239e+00; + xb = [x0 x1 x2 x3]'; + case 8 + x0 = 0.0000000000000e+00; + x1 = 3.8118550247622e-01; + x2 = 1.1899550868338e+00; + x3 = 2.2476300175641e+00; + x4 = 3.3192851303204e+00; + xb = [x0 x1 x2 x3 x4]'; + case 10 + x0 = 0.0000000000000e+00; + x1 = 3.5902433622052e-01; + x2 = 1.1436659188355e+00; + x3 = 2.2144895894456e+00; + x4 = 3.3682742337736e+00; + x5 = 4.4309689056870e+00; + xb = [x0 x1 x2 x3 x4 x5]'; + case 12 + x0 = 0.0000000000000e+00; + x1 = 3.6098032343909e-01; + x2 = 1.1634317168086e+00; + x3 = 2.2975905356987e+00; + x4 = 3.6057529790929e+00; + x5 = 4.8918275675510e+00; + x6 = 6.0000000000000e+00; + xb = [x0 x1 x2 x3 x4 x5 x6]'; + otherwise + error('Invalid operator order %d.',order); + end +end \ No newline at end of file
diff -r 8ec777fb473e -r 196123459178 +sbp/+grid/minimalBoundaryOptimizedGrid.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+grid/minimalBoundaryOptimizedGrid.m Wed Jul 08 18:22:54 2020 +0200 @@ -0,0 +1,60 @@ +% Computes the grid points x and grid spacing h used by the boundary optimized SBP operators +% with minimal number of non-equidistant boundary points, presented in +% 'Boundary optimized diagonal-norm SBP operators - Mattsson, Almquist, van der Weide 2018'. +% +% lim - cell array with domain limits +% N - Number of grid points +% order - order of accuracy of sbp operator. +function [x,h] = minimalBoundaryOptimizedGrid(lim,N,order) + assert(iscell(lim) && numel(lim) == 2,'The limit should be a cell array with 2 elements.'); + L = lim{2} - lim{1}; + assert(L>0,'Limits must be given in increasing order.'); + %%%% Non-equidistant grid points %%%%% + xb = boundaryPoints(order); + m = length(xb)-1; % Number of non-equidistant points + assert(N-2*(m+1)>=0,'Not enough grid points to contain the boundary region. Requires at least %d points.',2*(m+1)); + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Compute h %%%%%%%%%% + h = L/(2*xb(end) + N-1-2*m); + %%%%%%%%%%%%%%%%%%%%%%%%% + + %%%% Define grid %%%%%%%% + x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; + x = x + lim{1}; + %%%%%%%%%%%%%%%%%%%%%%%%% +end + +function xb = boundaryPoints(order) + switch order + case 4 + x0 = 0.0000000000000e+00; + x1 = 7.7122987842562e-01; + xb = [x0 x1]'; + case 6 + x0 = 0.0000000000000e+00; + x1 = 4.0842950991998e-01; + x2 = 1.1968523189207e+00; + xb = [x0 x1 x2]'; + case 8 + x0 = 0.0000000000000e+00; + x1 = 4.9439570885261e-01; + x2 = 1.4051531374839e+00; + xb = [x0 x1 x2]'; + case 10 + x0 = 0.0000000000000e+00; + x1 = 5.8556160757529e-01; + x2 = 1.7473267488572e+00; + x3 = 3.0000000000000e+00; + xb = [x0 x1 x2 x3]'; + case 12 + x0 = 0.0000000000000e+00; + x1 = 4.6552112904489e-01; + x2 = 1.4647984306493e+00; + x3 = 2.7620429464763e+00; + x4 = 4.0000000000000e+00; + xb = [x0 x1 x2 x3 x4]'; + otherwise + error('Invalid operator order %d.',order); + end +end \ No newline at end of file
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_10.m --- a/+sbp/+implementations/d1_noneq_10.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_10.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,48 +1,12 @@ -function [D1,H,x,h] = d1_noneq_10(N,L) +function [D1,H] = d1_noneq_10(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<20) error('Operator requires at least 20 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 10; -m = 5; -order = 10; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 3.5902433622052e-01; -x2 = 1.1436659188355e+00; -x3 = 2.2144895894456e+00; -x4 = 3.3682742337736e+00; -x5 = 4.4309689056870e+00; -x6 = 5.4309689056870e+00; -x7 = 6.4309689056870e+00; -x8 = 7.4309689056870e+00; -x9 = 8.4309689056870e+00; -x10 = 9.4309689056870e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -69,22 +33,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 10; +d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_12.m --- a/+sbp/+implementations/d1_noneq_12.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_12.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,50 +1,12 @@ -function [D1,H,x,h] = d1_noneq_12(N,L) +function [D1,H] = d1_noneq_12(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<24) error('Operator requires at least 24 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 12; -m = 6; -order = 12; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 3.6098032343909e-01; -x2 = 1.1634317168086e+00; -x3 = 2.2975905356987e+00; -x4 = 3.6057529790929e+00; -x5 = 4.8918275675510e+00; -x6 = 6.0000000000000e+00; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; -x9 = 9.0000000000000e+00; -x10 = 1.0000000000000e+01; -x11 = 1.1000000000000e+01; -x12 = 1.2000000000000e+01; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -73,22 +35,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 12; +d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_4.m --- a/+sbp/+implementations/d1_noneq_4.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_4.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,42 +1,12 @@ -function [D1,H,x,h] = d1_noneq_4(N,L) +function [D1,H] = d1_noneq_4(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<8) error('Operator requires at least 8 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 4; -m = 2; -order = 4; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 6.8764546205559e-01; -x2 = 1.8022115125776e+00; -x3 = 2.8022115125776e+00; -x4 = 3.8022115125776e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -57,22 +27,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 4; +d = [1/12,-2/3,0,2/3,-1/12]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_6.m --- a/+sbp/+implementations/d1_noneq_6.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_6.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,44 +1,12 @@ -function [D1,H,x,h] = d1_noneq_6(N,L) +function [D1,H] = d1_noneq_6(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<12) error('Operator requires at least 12 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 6; -m = 3; -order = 6; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.4090263368623e-01; -x2 = 1.2855984345073e+00; -x3 = 2.2638953951239e+00; -x4 = 3.2638953951239e+00; -x5 = 4.2638953951239e+00; -x6 = 5.2638953951239e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -61,22 +29,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 6; +d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_8.m --- a/+sbp/+implementations/d1_noneq_8.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_8.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,46 +1,12 @@ -function [D1,H,x,h] = d1_noneq_8(N,L) +function [D1,H] = d1_noneq_8(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<16) error('Operator requires at least 16 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 8; -m = 4; -order = 8; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 3.8118550247622e-01; -x2 = 1.1899550868338e+00; -x3 = 2.2476300175641e+00; -x4 = 3.3192851303204e+00; -x5 = 4.3192851303204e+00; -x6 = 5.3192851303204e+00; -x7 = 6.3192851303204e+00; -x8 = 7.3192851303204e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -65,22 +31,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 8; +d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_minimal_10.m --- a/+sbp/+implementations/d1_noneq_minimal_10.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_minimal_10.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,46 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_10(N,L) +function [D1,H] = d1_noneq_minimal_10(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<16) error('Operator requires at least 16 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 8; -m = 3; -order = 10; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 5.8556160757529e-01; -x2 = 1.7473267488572e+00; -x3 = 3.0000000000000e+00; -x4 = 4.0000000000000e+00; -x5 = 5.0000000000000e+00; -x6 = 6.0000000000000e+00; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -65,22 +31,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 10; +d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_minimal_12.m --- a/+sbp/+implementations/d1_noneq_minimal_12.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_minimal_12.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,48 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_12(N,L) +function [D1,H] = d1_noneq_minimal_12(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<20) error('Operator requires at least 20 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 10; -m = 4; -order = 12; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.6552112904489e-01; -x2 = 1.4647984306493e+00; -x3 = 2.7620429464763e+00; -x4 = 4.0000000000000e+00; -x5 = 5.0000000000000e+00; -x6 = 6.0000000000000e+00; -x7 = 7.0000000000000e+00; -x8 = 8.0000000000000e+00; -x9 = 9.0000000000000e+00; -x10 = 1.0000000000000e+01; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -69,22 +33,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 12; +d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_minimal_4.m --- a/+sbp/+implementations/d1_noneq_minimal_4.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_minimal_4.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,41 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_4(N,L) +function [D1,H] = d1_noneq_minimal_4(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<6) error('Operator requires at least 6 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 3; -m = 1; -order = 4; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 7.7122987842562e-01; -x2 = 1.7712298784256e+00; -x3 = 2.7712298784256e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -55,22 +26,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 4; +d = [1/12,-2/3,0,2/3,-1/12]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_minimal_6.m --- a/+sbp/+implementations/d1_noneq_minimal_6.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_minimal_6.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,43 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_6(N,L) +function [D1,H] = d1_noneq_minimal_6(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<10) error('Operator requires at least 10 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 5; -m = 2; -order = 6; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.0842950991998e-01; -x2 = 1.1968523189207e+00; -x3 = 2.1968523189207e+00; -x4 = 3.1968523189207e+00; -x5 = 4.1968523189207e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -59,22 +28,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 6; +d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/+implementations/d1_noneq_minimal_8.m --- a/+sbp/+implementations/d1_noneq_minimal_8.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/+implementations/d1_noneq_minimal_8.m Wed Jul 08 18:22:54 2020 +0200 @@ -1,44 +1,12 @@ -function [D1,H,x,h] = d1_noneq_minimal_8(N,L) +function [D1,H] = d1_noneq_minimal_8(N,h) -% L: Domain length % N: Number of grid points -if(nargin < 2) - L = 1; -end - if(N<12) error('Operator requires at least 12 grid points'); end % BP: Number of boundary points -% m: Number of nonequidistant spacings -% order: Accuracy of interior stencil BP = 6; -m = 2; -order = 8; - -%%%% Non-equidistant grid points %%%%% -x0 = 0.0000000000000e+00; -x1 = 4.9439570885261e-01; -x2 = 1.4051531374839e+00; -x3 = 2.4051531374839e+00; -x4 = 3.4051531374839e+00; -x5 = 4.4051531374839e+00; -x6 = 5.4051531374839e+00; - -xb = sparse(m+1,1); -for i = 0:m - xb(i+1) = eval(['x' num2str(i)]); -end -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Compute h %%%%%%%%%% -h = L/(2*xb(end) + N-1-2*m); -%%%%%%%%%%%%%%%%%%%%%%%%% - -%%%% Define grid %%%%%%%% -x = h*[xb; linspace(xb(end)+1,L/h-xb(end)-1,N-2*(m+1))'; L/h-flip(xb) ]; -%%%%%%%%%%%%%%%%%%%%%%%%% %%%% Norm matrix %%%%%%%% P = sparse(BP,1); @@ -61,22 +29,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%% %%%% Q matrix %%%%%%%%%%% - % interior stencil -switch order - case 2 - d = [-1/2,0,1/2]; - case 4 - d = [1/12,-2/3,0,2/3,-1/12]; - case 6 - d = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60]; - case 8 - d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; - case 10 - d = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260]; - case 12 - d = [1/5544,-1/385,1/56,-5/63,15/56,-6/7,0,6/7,-15/56,5/63,-1/56,1/385,-1/5544]; -end +order = 8; +d = [1/280,-4/105,1/5,-4/5,0,4/5,-1/5,4/105,-1/280]; d = repmat(d,N,1); Q = spdiags(d,-order/2:order/2,N,N);
diff -r 8ec777fb473e -r 196123459178 +sbp/D1Nonequidistant.m --- a/+sbp/D1Nonequidistant.m Wed Nov 20 22:24:06 2019 +0000 +++ b/+sbp/D1Nonequidistant.m Wed Jul 08 18:22:54 2020 +0200 @@ -19,58 +19,53 @@ % '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'} - + [x,h] = sbp.util.accurateBoundaryOptimizedGrid(lim,m,order); if order == 4 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_4(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_4(m,h); elseif order == 6 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_6(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_6(m,h); elseif order == 8 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_8(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_8(m,h); elseif order == 10 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_10(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_10(m,h); elseif order == 12 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_12(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_12(m,h); else error('Invalid operator order %d.',order); end case {'Minimal','minimal','M'} - + [x,h] = sbp.util.minimalBoundaryOptimizedGrid(lim,m,order); if order == 4 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_4(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_4(m,h); elseif order == 6 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_6(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_6(m,h); elseif order == 8 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_8(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_8(m,h); elseif order == 10 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_10(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_10(m,h); elseif order == 12 - [obj.D1,obj.H,obj.x,obj.h] = ... - sbp.implementations.d1_noneq_minimal_12(m,L); + [obj.D1,obj.H] = ... + sbp.implementations.d1_noneq_minimal_12(m,h); else error('Invalid operator order %d.',order); end end - - obj.x = obj.x + x_l; + obj.h = h; + obj.x = x; obj.e_l = sparse(m,1); obj.e_r = sparse(m,1);