Mercurial > repos > public > sbplib
changeset 758:0ef96fcdc028 feature/d1_staggered
Merge with feature grids.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 17 Jun 2018 15:29:46 -0700 |
| parents | 179f234f6cbf (diff) 43e64f5b388e (current diff) |
| children | 2645188489f6 |
| files | +grid/Cartesian.m |
| diffstat | 60 files changed, 4933 insertions(+), 143 deletions(-) [+] |
line wrap: on
line diff
--- a/+grid/Cartesian.m Mon Jun 04 17:48:19 2018 +0200 +++ b/+grid/Cartesian.m Sun Jun 17 15:29:46 2018 -0700 @@ -5,6 +5,7 @@ m % Number of points in each direction x % Cell array of vectors with node placement for each dimension. h % Spacing/Scaling + lim % Cell array of left and right boundaries for each dimension. end % General d dimensional grid with n points @@ -27,6 +28,7 @@ end obj.h = []; + obj.lim = []; end % n returns the number of points in the grid function o = N(obj)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/Staggered1d.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,99 @@ +classdef Staggered1d < grid.Structured + properties + grids % Cell array of grids + Ngrids % Number of grids + h % Interior grid spacing + d % Number of dimensions + end + + methods + + % Accepts multiple grids and combines them into a staggered grid + function obj = Staggered1d(varargin) + + obj.d = 1; + + obj.Ngrids = length(varargin); + obj.grids = cell(obj.Ngrids, 1); + for i = 1:obj.Ngrids + obj.grids{i} = varargin{i}; + end + + obj.h = []; + end + + % N returns the total number of points + function o = N(obj) + o = 0; + for i=1:obj.Ngrids + o = o+obj.grids{i}.size(); + end + end + + % D returns the spatial dimension of the grid + function o = D(obj) + o = obj.d; + end + + % size returns the number of points on the primal grid + function m = size(obj) + m = obj.grids{1}.size(); + end + + % points returns an n x 1 vector containing the coordinates for all points. + function X = points(obj) + X = []; + for i = 1:obj.Ngrids + X = [X; obj.grids{i}.points()]; + end + end + + % matrices returns a cell array with coordinates in matrix form. + % For 2d case these will have to be transposed to work with plotting routines. + function X = matrices(obj) + + % There is no 1d matrix data type in matlab, handle special case + X{1} = reshape(obj.points(), [obj.m 1]); + + end + + function h = scaling(obj) + if isempty(obj.h) + error('grid:Staggered1d:NoScalingSet', 'No scaling set') + end + + h = obj.h; + end + + % Restricts the grid function gf on obj to the subgrid g. + % Only works for even multiples + function gf = restrictFunc(obj, gf, g) + error('grid:Staggered1d:NotImplemented','This method does not exist yet') + end + + % Projects the grid function gf on obj to the grid g. + function gf = projectFunc(obj, gf, g) + error('grid:Staggered1d:NotImplemented','This method does not exist yet') + end + + % Return the names of all boundaries in this grid. + function bs = getBoundaryNames(obj) + switch obj.D() + case 1 + bs = {'l', 'r'}; + case 2 + bs = {'w', 'e', 's', 'n'}; + case 3 + bs = {'w', 'e', 's', 'n', 'd', 'u'}; + otherwise + error('not implemented'); + end + end + + % Return coordinates for the given boundary + function X = getBoundary(obj, name) + % Use boundaries of first grid + X = obj.grids{1}.getBoundary(name); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/blockEvalOn.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,30 @@ +% Useful for evaulating forcing functions with different functional expressions for each block +% f: cell array of function handles fi +% f_i = f_i(x1,y,...,t) +% t: time point. If not specified, it is assumed that the functions take only spatial arguments. +function gf = blockEvalOn(g, f, t) +default_arg('t',[]); + +grids = g.grids; +nBlocks = length(grids); +gf = cell(nBlocks,1); + +if isempty(t) + for i = 1:nBlocks + grid.evalOn(grids{i}, f{i} ); + end +else + dim = nargin(f{1}) - 1; + for i = 1:nBlocks + switch dim + case 1 + gf{i} = grid.evalOn(grids{i}, @(x)f{i}(x,t) ); + case 2 + gf{i} = grid.evalOn(grids{i}, @(x,y)f{i}(x,y,t) ); + case 3 + gf{i} = grid.evalOn(grids{i}, @(x,y,z)f{i}(x,y,z,t) ); + end + end +end + +gf = blockmatrix.toMatrix(gf); \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/primalDual1D.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,30 @@ +% Creates a 1D staggered grid of dimension length(m). +% over the interval xlims +% Primal grid: equidistant with m points. +% Dual grid: m + 1 points, h/2 spacing first point. +% Examples +% g = grid.primal_dual_1D(m, xlim) +% g = grid.primal_dual_1D(11, {0,1}) +function [g_primal, g_dual] = primalDual1D(m, xlims) + + if ~iscell(xlims) || numel(xlims) ~= 2 + error('grid:primalDual1D:InvalidLimits','The limits should be cell arrays with 2 elements.'); + end + + if xlims{1} > xlims{2} + error('grid:primalDual1D:InvalidLimits','The elements of the limit must be increasing.'); + end + + xl = xlims{1}; + xr = xlims{2}; + h = (xr-xl)/(m-1); + + % Primal grid + g_primal = grid.equidistant(m, xlims); + g_primal.h = h; + + % Dual grid + x = [xl; linspace(xl+h/2, xr-h/2, m-1)'; xr]; + g_dual = grid.Cartesian(x); + g_dual.h = h; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+grid/primalDual1DTest.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,32 @@ +function tests = primalDual1DTest() + tests = functiontests(localfunctions); +end + + +function testErrorInvalidLimits(testCase) + in = { + {10,{1}}, + {10,[0,1]}, + {10,{1,0}}, + }; + + for i = 1:length(in) + testCase.verifyError(@()grid.primalDual1D(in{i}{:}),'grid:primalDual1D:InvalidLimits',sprintf('in(%d) = %s',i,toString(in{i}))); + end +end + +function testCompiles(testCase) + in = { + {5, {0,1}}, + }; + + out = { + {[0; 0.25; 0.5; 0.75; 1], [0; 0.125; 0.375; 0.625; 0.875; 1]}, + }; + + for i = 1:length(in) + [gp, gd] = grid.primalDual1D(in{i}{:}); + testCase.verifyEqual(gp.points(),out{i}{1}); + testCase.verifyEqual(gd.points(),out{i}{2}); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Circle.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,98 @@ +classdef Circle < multiblock.DefCurvilinear + properties + r, c + + hs + r_arc + omega + end + + methods + function obj = Circle(r, c, hs) + default_arg('r', 1); + default_arg('c', [0; 0]); + default_arg('hs', 0.435); + + + % alpha = 0.75; + % hs = alpha*r/sqrt(2); + + % Square should not be a square, it should be an arc. The arc radius + % is chosen so that the three angles of the meshes are all equal. + % This gives that the (half)arc opening angle of should be omega = pi/12 + omega = pi/12; + r_arc = hs*(2*sqrt(2))/(sqrt(3)-1); % = hs* 1/sin(omega) + c_arc = c - [(1/(2-sqrt(3))-1)*hs; 0]; + + cir = parametrization.Curve.circle(c,r,[-pi/4 pi/4]); + + c2 = cir(0); + c3 = cir(1); + + s1 = [-hs; -hs]; + s2 = [ hs; -hs]; + s3 = [ hs; hs]; + s4 = [-hs; hs]; + + sp2 = parametrization.Curve.line(s2,c2); + sp3 = parametrization.Curve.line(s3,c3); + + Se1 = parametrization.Curve.circle(c_arc,r_arc,[-omega, omega]); + Se2 = Se1.rotate(c,pi/2); + Se3 = Se2.rotate(c,pi/2); + Se4 = Se3.rotate(c,pi/2); + + + S = parametrization.Ti(Se1,Se2,Se3,Se4).rotate_edges(-1); + + A = parametrization.Ti(sp2, cir, sp3.reverse, Se1.reverse); + B = A.rotate(c,1*pi/2).rotate_edges(-1); + C = A.rotate(c,2*pi/2).rotate_edges(-1); + D = A.rotate(c,3*pi/2).rotate_edges(0); + + blocks = {S,A,B,C,D}; + blocksNames = {'S','A','B','C','D'}; + + conn = cell(5,5); + conn{1,2} = {'e','w'}; + conn{1,3} = {'n','s'}; + conn{1,4} = {'w','s'}; + conn{1,5} = {'s','w'}; + + conn{2,3} = {'n','e'}; + conn{3,4} = {'w','e'}; + conn{4,5} = {'w','s'}; + conn{5,2} = {'n','s'}; + + boundaryGroups = struct(); + boundaryGroups.E = multiblock.BoundaryGroup({2,'e'}); + boundaryGroups.N = multiblock.BoundaryGroup({3,'n'}); + boundaryGroups.W = multiblock.BoundaryGroup({4,'n'}); + boundaryGroups.S = multiblock.BoundaryGroup({5,'e'}); + boundaryGroups.all = multiblock.BoundaryGroup({{2,'e'},{3,'n'},{4,'n'},{5,'e'}}); + + obj = obj@multiblock.DefCurvilinear(blocks, conn, boundaryGroups, blocksNames); + + obj.r = r; + obj.c = c; + obj.hs = hs; + obj.r_arc = r_arc; + obj.omega = omega; + end + + function ms = getGridSizes(obj, m) + m_S = m; + + % m_Radial + s = 2*obj.hs; + innerArc = obj.r_arc*obj.omega; + outerArc = obj.r*pi/2; + shortSpoke = obj.r-s/sqrt(2); + x = (1/(2-sqrt(3))-1)*obj.hs; + longSpoke = (obj.r+x)-obj.r_arc; + m_R = parametrization.equal_step_size((innerArc+outerArc)/2, m_S, (shortSpoke+longSpoke)/2); + + ms = {[m_S m_S], [m_R m_S], [m_S m_R], [m_S m_R], [m_R m_S]}; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Rectangle.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,188 @@ +classdef Rectangle < multiblock.Definition + properties + + blockTi % Transfinite interpolation objects used for plotting + xlims + ylims + blockNames % Cell array of block labels + nBlocks + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + + end + + + methods + % Creates a divided rectangle + % x and y are vectors of boundary and interface positions. + % blockNames: cell array of labels. The id is default. + function obj = Rectangle(x,y,blockNames) + default_arg('blockNames',[]); + + n = length(y)-1; % number of blocks in the y direction. + m = length(x)-1; % number of blocks in the x direction. + N = n*m; % number of blocks + + if ~issorted(x) + error('The elements of x seem to be in the wrong order'); + end + if ~issorted(flip(y)) + error('The elements of y seem to be in the wrong order'); + end + + % Dimensions of blocks and number of points + blockTi = cell(N,1); + xlims = cell(N,1); + ylims = cell(N,1); + for i = 1:n + for j = 1:m + p1 = [x(j), y(i+1)]; + p2 = [x(j+1), y(i)]; + I = flat_index(m,j,i); + blockTi{I} = parametrization.Ti.rectangle(p1,p2); + xlims{I} = {x(j), x(j+1)}; + ylims{I} = {y(i+1), y(i)}; + end + end + + % Interface couplings + conn = cell(N,N); + for i = 1:n + for j = 1:m + I = flat_index(m,j,i); + if i < n + J = flat_index(m,j,i+1); + conn{I,J} = {'s','n'}; + end + + if j < m + J = flat_index(m,j+1,i); + conn{I,J} = {'e','w'}; + end + end + end + + % Block names (id number as default) + if isempty(blockNames) + obj.blockNames = cell(1, N); + for i = 1:N + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == N); + obj.blockNames = blockNames; + end + nBlocks = N; + + % Boundary groups + boundaryGroups = struct(); + nx = m; + ny = n; + E = cell(1,ny); + W = cell(1,ny); + S = cell(1,nx); + N = cell(1,nx); + for i = 1:ny + E_id = flat_index(m,nx,i); + W_id = flat_index(m,1,i); + E{i} = {E_id,'e'}; + W{i} = {W_id,'w'}; + end + for j = 1:nx + S_id = flat_index(m,j,ny); + N_id = flat_index(m,j,1); + S{j} = {S_id,'s'}; + N{j} = {N_id,'n'}; + end + boundaryGroups.E = multiblock.BoundaryGroup(E); + boundaryGroups.W = multiblock.BoundaryGroup(W); + boundaryGroups.S = multiblock.BoundaryGroup(S); + boundaryGroups.N = multiblock.BoundaryGroup(N); + boundaryGroups.all = multiblock.BoundaryGroup([E,W,S,N]); + boundaryGroups.WS = multiblock.BoundaryGroup([W,S]); + boundaryGroups.WN = multiblock.BoundaryGroup([W,N]); + boundaryGroups.ES = multiblock.BoundaryGroup([E,S]); + boundaryGroups.EN = multiblock.BoundaryGroup([E,N]); + + obj.connections = conn; + obj.nBlocks = nBlocks; + obj.boundaryGroups = boundaryGroups; + obj.blockTi = blockTi; + obj.xlims = xlims; + obj.ylims = ylims; + + end + + + % 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. + function g = getGrid(obj, ms, varargin) + + default_arg('ms',[21,21]) + + % Extend ms if input is a single vector + if (numel(ms) == 2) && ~iscell(ms) + m = ms; + ms = cell(1,obj.nBlocks); + for i = 1:obj.nBlocks + ms{i} = m; + end + end + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistant(ms{i}, obj.xlims{i}, obj.ylims{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockTi{i}.show(2,2); + end + axis equal + return + end + + if gridLines + m = 10; + for i = 1:obj.nBlocks + obj.blockTi{i}.show(m,m); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockTi{i}, labels{i}); + end + + axis equal + end + + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + function ms = getGridSizes(obj, varargin) + end + end +end
--- a/+multiblock/Def.m Mon Jun 04 17:48:19 2018 +0200 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,101 +0,0 @@ -classdef Def - properties - nBlocks - blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation - blockNames - connections % Cell array specifying connections between blocks - boundaryGroups % Structure of boundaryGroups - end - - methods - % Defines a multiblock setup for transfinite interpolation blocks - % TODO: How to bring in plotting of points? - function obj = Def(blockMaps, connections, boundaryGroups, blockNames) - default_arg('boundaryGroups', struct()); - default_arg('blockNames',{}); - - nBlocks = length(blockMaps); - - obj.nBlocks = nBlocks; - - obj.blockMaps = blockMaps; - - assert(all(size(connections) == [nBlocks, nBlocks])); - obj.connections = connections; - - - if isempty(blockNames) - obj.blockNames = cell(1, nBlocks); - for i = 1:length(blockMaps) - obj.blockNames{i} = sprintf('%d', i); - end - else - assert(length(blockNames) == nBlocks); - obj.blockNames = blockNames; - end - - obj.boundaryGroups = boundaryGroups; - end - - function g = getGrid(obj, varargin) - ms = obj.getGridSizes(varargin{:}); - - grids = cell(1, obj.nBlocks); - for i = 1:obj.nBlocks - grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); - end - - g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); - end - - function show(obj, label, gridLines, varargin) - default_arg('label', 'name') - default_arg('gridLines', false); - - if isempty('label') && ~gridLines - for i = 1:obj.nBlocks - obj.blockMaps{i}.show(2,2); - end - axis equal - return - end - - if gridLines - ms = obj.getGridSizes(varargin{:}); - for i = 1:obj.nBlocks - obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); - end - end - - - switch label - case 'name' - labels = obj.blockNames; - case 'id' - labels = {}; - for i = 1:obj.nBlocks - labels{i} = num2str(i); - end - otherwise - axis equal - return - end - - for i = 1:obj.nBlocks - parametrization.Ti.label(obj.blockMaps{i}, labels{i}); - end - - axis equal - end - end - - methods (Abstract) - % Returns the grid size of each block in a cell array - % The input parameters are determined by the subclass - ms = getGridSizes(obj, varargin) - % end - end - -end - -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/DefCurvilinear.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,101 @@ +classdef DefCurvilinear < multiblock.Definition + properties + nBlocks + blockMaps % Maps from logical blocks to physical blocks build from transfinite interpolation + blockNames + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + end + + methods + % Defines a multiblock setup for transfinite interpolation blocks + % TODO: How to bring in plotting of points? + function obj = DefCurvilinear(blockMaps, connections, boundaryGroups, blockNames) + default_arg('boundaryGroups', struct()); + default_arg('blockNames',{}); + + nBlocks = length(blockMaps); + + obj.nBlocks = nBlocks; + + obj.blockMaps = blockMaps; + + assert(all(size(connections) == [nBlocks, nBlocks])); + obj.connections = connections; + + + if isempty(blockNames) + obj.blockNames = cell(1, nBlocks); + for i = 1:length(blockMaps) + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == nBlocks); + obj.blockNames = blockNames; + end + + obj.boundaryGroups = boundaryGroups; + end + + function g = getGrid(obj, varargin) + ms = obj.getGridSizes(varargin{:}); + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistantCurvilinear(obj.blockMaps{i}.S, ms{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + function show(obj, label, gridLines, varargin) + default_arg('label', 'name') + default_arg('gridLines', false); + + if isempty('label') && ~gridLines + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(2,2); + end + axis equal + return + end + + if gridLines + ms = obj.getGridSizes(varargin{:}); + for i = 1:obj.nBlocks + obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); + end + end + + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + for i = 1:obj.nBlocks + parametrization.Ti.label(obj.blockMaps{i}, labels{i}); + end + + axis equal + end + end + + methods (Abstract) + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + ms = getGridSizes(obj, varargin) + % end + end + +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/Definition.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,11 @@ +classdef Definition + methods (Abstract) + + % Returns a multiblock.Grid given some parameters + g = getGrid(obj, varargin) + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + show(obj, label, gridLines, varargin) + end +end
--- a/+multiblock/Grid.m Mon Jun 04 17:48:19 2018 +0200 +++ b/+multiblock/Grid.m Sun Jun 17 15:29:46 2018 -0700 @@ -77,7 +77,7 @@ % Collect number of points in each block N = zeros(1,nBlocks); for i = 1:nBlocks - N(i) = obj.grids{i}.N(); + N(i) = obj.grids{i}.N()*nComponents; end gfs = blockmatrix.fromMatrix(gf, {N,1});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_2.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,71 @@ +function [xp, xm, Pp, Pm, Qp, Qm] = d1_staggered_2(n, L) +% [xp, xm, Pp, Pm, Qp, Qm] = sbp_staggered_2(n, L) +% Input arguments: +% n : number of grid points: n on xp grid and n+1 on xm grid +% L : interval length + +h = L/(n-1); + +n = n - 1; +assert(n >= 3, 'Not enough grid points'); + +qm20 = -1/4; +qm01 = 1/2; +qp11 = -1/4; +qm11 = 1/4; +qp02 = 1/4; +pp0 = 1/2; +pm1 = 1/4; +qm00 = -1/2; +qp10 = -1/2; +qp12 = 3/4; +qp00 = -1/2; +qm21 = -3/4; +pp1 = 1; +qp01 = 1/4; +pm0 = 1/2; +pm2 = 5/4; +qm10 = -1/4; + +% Q+ and Q-, top-left corner +QpL = [... +qp00, qp01, qp02; + qp10, qp11, qp12 +]; +QmL = [... +qm00, qm01; + qm10, qm11; + qm20, qm21 +]; + +% Q+ and Q- +b = 2; +w = b; +s = [-1, 1]'; +Qp = spdiags(repmat(-s(end:-1:1)',[n+2 1]), -(w/2-1):w/2, n+2, n+2); +Qm = spdiags(repmat(s(:)',[n+2 1]), -(w/2-1)-1:w/2-1, n+2, n+2); +Qp(end,:) = []; +Qm(:,end) = []; + +% Add SBP boundary closures +Qp(1:b,1:b+1) = QpL; +Qp(end-b+1:end,end-b:end) = -fliplr(flipud(QpL)); +Qm(1:b+1,1:b) = QmL; +Qm(end-b:end,end-b+1:end) = -fliplr(flipud(QmL)); + +% P+ and P- +Pp = ones(n+1,1); +Pm = ones(n+2,1); + +Pp(1:b) = [pp0, pp1]; +Pp(end-b+1:end) = Pp(b:-1:1); +Pm(1:b+1) = [pm0, pm1, pm2]; +Pm(end-b:end) = Pm(b+1:-1:1); +Pp = spdiags(Pp,0,n+1,n+1); +Pm = spdiags(Pm,0,n+2,n+2); + +Pp = h*Pp; +Pm = h*Pm; + +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_4.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,58 @@ +function [xp, xm, Pp, Pm, Qp, Qm] = d1_staggered_4(n, L) +% [xp, xm, Pp, Pm, Qp, Qm] = sbp_staggered_4(n, L) +% Input arguments: +% n : number of grid points: n on xp grid and n+1 on xm grid +% L : interval length + +h = L/(n-1); + +n = n - 1; +assert(n >= 7, 'Not enough grid points'); + +% Q+ and Q-, top-left corner +QpL = [... +-323/378, 2783/3072, -1085/27648, -17/9216, -667/64512; + -5/21, -847/1024, 1085/1024, -37/1024, 899/21504; + 5/42, -121/1024, -1085/1024, 3575/3072, -753/7168; + -5/189, 121/3072, 1085/27648, -10759/9216, 74635/64512 +]; +QmL = [... +-55/378, 5/21, -5/42, 5/189; + -2783/3072, 847/1024, 121/1024, -121/3072; + 1085/27648, -1085/1024, 1085/1024, -1085/27648; + 17/9216, 37/1024, -3575/3072, 10759/9216; + 667/64512, -899/21504, 753/7168, -74635/64512 +]; + +% Q+ and Q- +w = 4; +s = [ 1/24, -9/8, 9/8, -1/24]; +Qp = spdiags(repmat(-s(end:-1:1),[n+2 1]), -(w/2-1):w/2, n+2, n+2); +Qm = spdiags(repmat(s(:)',[n+2 1]), -(w/2-1)-1:w/2-1, n+2, n+2); +Qp(end,:) = []; +Qm(:,end) = []; + +% Add SBP boundary closures +bp = 4; +bm = 5; +Qp(1:bp,1:bm) = double(QpL); +Qp(end-bp+1:end,end-bm+1:end) = -fliplr(flipud(QpL)); +Qm(1:bm,1:bp) = double(QmL); +Qm(end-bm+1:end,end-bp+1:end) = -fliplr(flipud(QmL)); + +% P+ and P- +Pp = ones(n+1,1); +Pm = ones(n+2,1); + +Pp(1:bp) = [407/1152, 473/384, 343/384, 1177/1152]; +Pp(end-bp+1:end) = Pp(bp:-1:1); +Pm(1:bm) = [5/63, 121/128, 1085/1152, 401/384, 2659/2688]; +Pm(end-bm+1:end) = Pm(bm:-1:1); +Pp = spdiags(Pp,0,n+1,n+1); +Pm = spdiags(Pm,0,n+2,n+2); + +Pp = h*Pp; +Pm = h*Pm; + +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_6.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,103 @@ +function [xp, xm, Pp, Pm, Qp, Qm] = d1_staggered_6(n, L) +% [xp, xm, Pp, Pm, Qp, Qm] = sbp_staggered_6(n, L) +% Input arguments: +% n : number of grid points: n on xp grid and n+1 on xm grid +% L : interval length + +h = L/(n-1); + +n = n - 1; +assert(n >= 9, 'Not enough grid points'); + +% Coefficients determined such that the SBP property is satisfied +qp44 = -144662434693018337/131093333881896960; +pp3 = 793/720; +qp22 = -729260819980499/780317463582720; +pm5 = 1262401/1244160; +pp4 = 157/160; +qp03 = -1056618403248889/17835827739033600; +qp32 = 17518644727427/2340952390748160; +pm4 = 224201/241920; +pm2 = 90101/103680; +qp40 = -572605772036807/27311111225395200; +qp30 = 3572646778087103/61450000257139200; +pm3 = 26369/23040; +qp41 = 935926732952051/29964190601576448; +qp02 = -67953322736279/2340952390748160; +pp0 = 95/288; +qp01 = 143647626126052207/149820953007882240; +pp1 = 317/240; +qp12 = 444281115722813/468190478149632; +qp25 = 312883915256963/24970158834647040; +qp11 = -36400588465707331/37455238251970560; +qp20 = -137012145866261/13655555612697600; +qp10 = -1008199398207353/6827777806348800; +qp15 = -1088192686015471/37455238251970560; +qp04 = -239061311631881/131093333881896960; +qp35 = -2480203967217637/112365714755911680; +qp21 = 1658590350384343/24970158834647040; +qp05 = 4843362440877649/449462859023646720; +qp24 = -372310079307409/4369777796063232; +qp33 = -4809914536549211/4458956934758400; +qp23 = 2826535075774271/2972637956505600; +qp42 = 4202889834071/585238097687040; +qp14 = 492703999899311/32773333470474240; +pm1 = 133801/138240; +qp34 = 36544940032227659/32773333470474240; +qp00 = -216175739231516543/245800001028556800; +qp43 = 128340370711031/17835827739033600; +qp13 = 823082791653949/4458956934758400; +pm0 = 5251/68040; +pp2 = 23/30; +qp45 = 170687698457070187/149820953007882240; +qp31 = -3169112007572299/37455238251970560; + +% Q+ and Q-, top-left corner +QpL = [... +-216175739231516543/245800001028556800, 143647626126052207/149820953007882240, -67953322736279/2340952390748160, -1056618403248889/17835827739033600, -239061311631881/131093333881896960, 4843362440877649/449462859023646720; + -1008199398207353/6827777806348800, -36400588465707331/37455238251970560, 444281115722813/468190478149632, 823082791653949/4458956934758400, 492703999899311/32773333470474240, -1088192686015471/37455238251970560; + -137012145866261/13655555612697600, 1658590350384343/24970158834647040, -729260819980499/780317463582720, 2826535075774271/2972637956505600, -372310079307409/4369777796063232, 312883915256963/24970158834647040; + 3572646778087103/61450000257139200, -3169112007572299/37455238251970560, 17518644727427/2340952390748160, -4809914536549211/4458956934758400, 36544940032227659/32773333470474240, -2480203967217637/112365714755911680; + -572605772036807/27311111225395200, 935926732952051/29964190601576448, 4202889834071/585238097687040, 128340370711031/17835827739033600, -144662434693018337/131093333881896960, 170687698457070187/149820953007882240 +]; +QmL = [... +-29624261797040257/245800001028556800, 1008199398207353/6827777806348800, 137012145866261/13655555612697600, -3572646778087103/61450000257139200, 572605772036807/27311111225395200; + -143647626126052207/149820953007882240, 36400588465707331/37455238251970560, -1658590350384343/24970158834647040, 3169112007572299/37455238251970560, -935926732952051/29964190601576448; + 67953322736279/2340952390748160, -444281115722813/468190478149632, 729260819980499/780317463582720, -17518644727427/2340952390748160, -4202889834071/585238097687040; + 1056618403248889/17835827739033600, -823082791653949/4458956934758400, -2826535075774271/2972637956505600, 4809914536549211/4458956934758400, -128340370711031/17835827739033600; + 239061311631881/131093333881896960, -492703999899311/32773333470474240, 372310079307409/4369777796063232, -36544940032227659/32773333470474240, 144662434693018337/131093333881896960; + -4843362440877649/449462859023646720, 1088192686015471/37455238251970560, -312883915256963/24970158834647040, 2480203967217637/112365714755911680, -170687698457070187/149820953007882240 +]; + +% Q+ and Q- +w = 6; +s = [ -3/640, 25/384, -75/64, 75/64, -25/384, 3/640]; +Qp = spdiags(repmat(-s(end:-1:1),[n+2 1]), -(w/2-1):w/2, n+2, n+2); +Qm = spdiags(repmat(s(:)',[n+2 1]), -(w/2-1)-1:w/2-1, n+2, n+2); +Qp(end,:) = []; +Qm(:,end) = []; + +% Add SBP boundary closures +bp = 5; +bm = 6; +Qp(1:bp,1:bm) = double(QpL); +Qp(end-bp+1:end,end-bm+1:end) = -fliplr(flipud(QpL)); +Qm(1:bm,1:bp) = double(QmL); +Qm(end-bm+1:end,end-bp+1:end) = -fliplr(flipud(QmL)); + +% P+ and P- +Pp = ones(n+1,1); +Pm = ones(n+2,1); + +Pp(1:bp) = [95/288, 317/240, 23/30, 793/720, 157/160]; +Pp(end-bp+1:end) = Pp(bp:-1:1); +Pm(1:bm) = [5251/68040, 133801/138240, 90101/103680, 26369/23040, 224201/241920, 1262401/1244160]; +Pm(end-bm+1:end) = Pm(bm:-1:1); +Pp = spdiags(Pp,0,n+1,n+1); +Pm = spdiags(Pm,0,n+2,n+2); + +Pp = h*Pp; +Pm = h*Pm; + +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_upwind_2.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,73 @@ +function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_2(n,L) + +assert(n-1 >= 6,'Not enough grid points'); + +np=n; +nm=n+1; +h = L/(n-1); + + +%H- norm +Hm_U=[0.88829e5 / 0.599299e6 0 0; 0 0.8056187e7 / 0.9588784e7 0; 0 0 0.9700117e7 / 0.9588784e7;]; +%H+ norm +Hp_U=[0.7489557e7 / 0.19177568e8 0 0; 0 0.1386089e7 / 0.1198598e7 0; 0 0 0.18276939e8 / 0.19177568e8;]; +% Upper part of Qpp. Notice that the B matrix is not included here. +% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0 +Qpp_U=[-0.1e1 / 0.32e2 0.12886609e8 / 0.19177568e8 -0.1349263e7 / 0.9588784e7; -0.10489413e8 / 0.19177568e8 -0.3e1 / 0.16e2 0.16482403e8 / 0.19177568e8; 0.187491e6 / 0.2397196e7 -0.9290815e7 / 0.19177568e8 -0.11e2 / 0.32e2;]; +% Upper part of Qmp. Notice that the B matrix is not included here. +% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0 +Qmp_U=[-0.2e1 / 0.21e2 0.12495263e8 / 0.16780372e8 -0.7520839e7 / 0.50341116e8; -0.7700871e7 / 0.16780372e8 -0.31e2 / 0.112e3 0.57771939e8 / 0.67121488e8; 0.2726447e7 / 0.50341116e8 -0.31402783e8 / 0.67121488e8 -0.113e3 / 0.336e3;]; +% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0 +Qp_U=[-0.801195e6 / 0.2397196e7 0.16507959e8 / 0.19177568e8 -0.509615e6 / 0.19177568e8; -0.219745e6 / 0.1198598e7 -0.2112943e7 / 0.2397196e7 0.2552433e7 / 0.2397196e7; 0.42087e5 / 0.2397196e7 0.395585e6 / 0.19177568e8 -0.19909849e8 / 0.19177568e8;]; +Hp=spdiags(ones(np,1),0,np,np); +Hp(1:3,1:3)=Hp_U; +Hp(np-2:np,np-2:np)=fliplr(flipud(Hp_U)); +Hp=Hp*h; +HIp=inv(Hp); + +Hm=spdiags(ones(nm,1),0,nm,nm); +Hm(1:3,1:3)=Hm_U; +Hm(nm-2:nm,nm-2:nm)=fliplr(flipud(Hm_U)); +Hm=Hm*h; +HIm=inv(Hm); + +Qpp=spdiags(repmat([-0.3e1 / 0.8e1 -0.3e1 / 0.8e1 0.7e1 / 0.8e1 -0.1e1 / 0.8e1;],[np,1]),-1:2,np,np); +Qpp(1:3,1:3)=Qpp_U; +Qpp(np-2:np,np-2:np)=flipud( fliplr(Qpp_U(1:3,1:3) ) )'; + +Qpm=-Qpp'; + +Qmp=spdiags(repmat([-0.3e1 / 0.8e1 -0.3e1 / 0.8e1 0.7e1 / 0.8e1 -0.1e1 / 0.8e1;],[nm,1]),-1:2,nm,nm); +Qmp(1:3,1:3)=Qmp_U; +Qmp(nm-2:nm,nm-2:nm)=flipud( fliplr(Qmp_U(1:3,1:3) ) )'; + +Qmm=-Qmp'; + + +Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1; +Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1; + +Dpp=HIp*(Qpp+1/2*Bpp) ; +Dpm=HIp*(Qpm+1/2*Bpp) ; + + +Dmp=HIm*(Qmp+1/2*Bmp) ; +Dmm=HIm*(Qmm+1/2*Bmp) ; + + +%%% Start with the staggered +Qp=spdiags(repmat([-1 1],[np,1]),0:1,np,nm); +Qp(1:3,1:3)=Qp_U; +Qp(np-2:np,nm-2:nm)=flipud( fliplr(-Qp_U(1:3,1:3) ) ); +Qm=-Qp'; + +Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1; +Bm=Bp'; + +Dp=HIp*(Qp+1/2*Bp) ; + +Dm=HIm*(Qm+1/2*Bm) ; + +% grids +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_upwind_4.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,76 @@ +function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_4(n,L) + +assert(n-1 >= 8,'Not enough grid points'); + +np=n; +nm=n+1; +h = L/(n-1); + + +%H- norm +Hm_U=[0.91e2 / 0.720e3 0 0 0; 0 0.325e3 / 0.384e3 0 0; 0 0 0.595e3 / 0.576e3 0; 0 0 0 0.1909e4 / 0.1920e4;]; +%H+ norm +Hp_U=[0.805e3 / 0.2304e4 0 0 0; 0 0.955e3 / 0.768e3 0 0; 0 0 0.677e3 / 0.768e3 0; 0 0 0 0.2363e4 / 0.2304e4;]; +% Upper part of Qpp. Notice that the B matrix is not included here. +% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0 +Qpp_U=[-0.11e2 / 0.1536e4 0.1493e4 / 0.2304e4 -0.571e3 / 0.4608e4 -0.13e2 / 0.768e3; -0.697e3 / 0.1152e4 -0.121e3 / 0.1536e4 0.97e2 / 0.128e3 -0.473e3 / 0.4608e4; 0.373e3 / 0.4608e4 -0.139e3 / 0.256e3 -0.319e3 / 0.1536e4 0.1999e4 / 0.2304e4; 0.1e1 / 0.32e2 -0.121e3 / 0.4608e4 -0.277e3 / 0.576e3 -0.143e3 / 0.512e3;]; + +% Upper part of Qmp. Notice that the B matrix is not included here. +% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0 +Qmp_U=[-0.209e3 / 0.5970e4 0.13439e5 / 0.17910e5 -0.13831e5 / 0.47760e5 0.10637e5 / 0.143280e6; -0.44351e5 / 0.71640e5 -0.20999e5 / 0.152832e6 0.230347e6 / 0.229248e6 -0.70547e5 / 0.254720e6; 0.3217e4 / 0.15920e5 -0.86513e5 / 0.114624e6 -0.15125e5 / 0.76416e5 0.1087241e7 / 0.1146240e7; -0.1375e4 / 0.28656e5 0.36117e5 / 0.254720e6 -0.655601e6 / 0.1146240e7 -0.211717e6 / 0.764160e6;]; + +% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0 +Qp_U=[-0.338527e6 / 0.1004160e7 0.4197343e7 / 0.4819968e7 0.1423e4 / 0.803328e6 -0.854837e6 / 0.24099840e8; -0.520117e6 / 0.3012480e7 -0.492581e6 / 0.535552e6 0.2476673e7 / 0.2409984e7 0.520117e6 / 0.8033280e7; -0.50999e5 / 0.3012480e7 0.117943e6 / 0.1606656e7 -0.2476673e7 / 0.2409984e7 0.2712193e7 / 0.2677760e7; 0.26819e5 / 0.1004160e7 -0.117943e6 / 0.4819968e7 -0.1423e4 / 0.803328e6 -0.26119411e8 / 0.24099840e8;]; + +Hp=spdiags(ones(np,1),0,np,np); +Hp(1:4,1:4)=Hp_U; +Hp(np-3:np,np-3:np)=fliplr(flipud(Hp_U)); +Hp=Hp*h; +HIp=inv(Hp); + +Hm=spdiags(ones(nm,1),0,nm,nm); +Hm(1:4,1:4)=Hm_U; +Hm(nm-3:nm,nm-3:nm)=fliplr(flipud(Hm_U)); +Hm=Hm*h; +HIm=inv(Hm); + +Qpp=spdiags(repmat([0.7e1 / 0.128e3 -0.67e2 / 0.128e3 -0.55e2 / 0.192e3 0.61e2 / 0.64e2 -0.29e2 / 0.128e3 0.11e2 / 0.384e3;],[np,1]),-2:3,np,np); +Qpp(1:4,1:4)=Qpp_U; +Qpp(np-3:np,np-3:np)=flipud( fliplr(Qpp_U(1:4,1:4) ) )'; + +Qpm=-Qpp'; + +Qmp=spdiags(repmat([0.7e1 / 0.128e3 -0.67e2 / 0.128e3 -0.55e2 / 0.192e3 0.61e2 / 0.64e2 -0.29e2 / 0.128e3 0.11e2 / 0.384e3;],[nm,1]),-2:3,nm,nm); +Qmp(1:4,1:4)=Qmp_U; +Qmp(nm-3:nm,nm-3:nm)=flipud( fliplr(Qmp_U(1:4,1:4) ) )'; + +Qmm=-Qmp'; + + +Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1; +Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1; + +Dpp=HIp*(Qpp+1/2*Bpp) ; +Dpm=HIp*(Qpm+1/2*Bpp) ; + + +Dmp=HIm*(Qmp+1/2*Bmp) ; +Dmm=HIm*(Qmm+1/2*Bmp) ; + + +%%% Start with the staggered +Qp=spdiags(repmat([1/24 -9/8 9/8 -1/24],[np,1]),-1:2,np,nm); +Qp(1:4,1:4)=Qp_U; +Qp(np-3:np,nm-3:nm)=flipud( fliplr(-Qp_U(1:4,1:4) ) ); +Qm=-Qp'; + +Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1; +Bm=Bp'; + +Dp=HIp*(Qp+1/2*Bp) ; + +Dm=HIm*(Qm+1/2*Bm) ; + +% grids +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_upwind_6.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,78 @@ +function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_6(n,L) + +assert(n-1 >= 12,'Not enough grid points'); + +np=n; +nm=n+1; +h = L/(n-1); + +bp=6; % Number of boundary points + + +%H- norm +Hm_U=[0.18373e5 / 0.136080e6 0 0 0 0 0; 0 0.228247e6 / 0.276480e6 0 0 0 0; 0 0 0.219557e6 / 0.207360e6 0 0 0; 0 0 0 0.44867e5 / 0.46080e5 0 0; 0 0 0 0 0.487757e6 / 0.483840e6 0; 0 0 0 0 0 0.2485447e7 / 0.2488320e7;]; +%H+ norm +Hp_U=[0.174529e6 / 0.552960e6 0 0 0 0 0; 0 0.769723e6 / 0.552960e6 0 0 0 0; 0 0 0.172613e6 / 0.276480e6 0 0 0; 0 0 0 0.343867e6 / 0.276480e6 0 0; 0 0 0 0 0.503237e6 / 0.552960e6 0; 0 0 0 0 0 0.560831e6 / 0.552960e6;]; +% Upper part of Qpp. Notice that the B matrix is not included here. +% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0 +Qpp_U=[-0.11e2 / 0.12288e5 0.3248629e7 / 0.4976640e7 -0.742121e6 / 0.9953280e7 -0.173089e6 / 0.1658880e7 0.100627e6 / 0.9953280e7 0.263e3 / 0.15552e5; -0.3212989e7 / 0.4976640e7 -0.341e3 / 0.20480e5 0.488429e6 / 0.995328e6 0.2108717e7 / 0.9953280e7 0.5623e4 / 0.829440e6 -0.468779e6 / 0.9953280e7; 0.635201e6 / 0.9953280e7 -0.2135641e7 / 0.4976640e7 -0.2233e4 / 0.30720e5 0.298757e6 / 0.497664e6 -0.2148901e7 / 0.9953280e7 0.99581e5 / 0.1658880e7; 0.184969e6 / 0.1658880e7 -0.2671829e7 / 0.9953280e7 -0.1044721e7 / 0.2488320e7 -0.4697e4 / 0.30720e5 0.882599e6 / 0.995328e6 -0.2114063e7 / 0.9953280e7; -0.118447e6 / 0.9953280e7 0.15761e5 / 0.829440e6 0.915757e6 / 0.9953280e7 -0.2923243e7 / 0.4976640e7 -0.4279e4 / 0.20480e5 0.4614511e7 / 0.4976640e7; -0.263e3 / 0.15552e5 0.422447e6 / 0.9953280e7 -0.5185e4 / 0.331776e6 0.424727e6 / 0.9953280e7 -0.570779e6 / 0.995328e6 -0.2761e4 / 0.12288e5;]; + +% Upper part of Qmp. Notice that the B matrix is not included here. +% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0 +Qmp_U=[-0.6660404399e10 / 0.975680535935e12 0.42802131970831759e17 / 0.60695134779444480e17 -0.27241603626152813e17 / 0.91042702169166720e17 0.148772145985039e15 / 0.1264481974571760e16 -0.386865438537449e15 / 0.30347567389722240e17 -0.739491571084877e15 / 0.182085404338333440e18; -0.40937739835891039e17 / 0.60695134779444480e17 -0.20934998251893e14 / 0.570912496455680e15 0.3069347264824655e16 / 0.3082927480860672e16 -0.531249064089227e15 / 0.1541463740430336e16 0.6343721417240047e16 / 0.107902461830123520e18 0.194756943144697e15 / 0.138731736638730240e18; 0.24212381292051773e17 / 0.91042702169166720e17 -0.13932096926615587e17 / 0.15414637404303360e17 -0.22149140293797e14 / 0.285456248227840e15 0.443549615179363e15 / 0.481707418884480e15 -0.1531771257444041e16 / 0.5994581212784640e16 0.23579222779798361e17 / 0.416195209916190720e18; -0.119651391006031e15 / 0.1264481974571760e16 0.2061363806050549e16 / 0.7707318702151680e16 -0.355212104634871e15 / 0.481707418884480e15 -0.42909037900311e14 / 0.285456248227840e15 0.25466291778687943e17 / 0.26975615457530880e17 -0.3289076301679109e16 / 0.11560978053227520e17; 0.153994229603129e15 / 0.30347567389722240e17 -0.2655886084488631e16 / 0.107902461830123520e18 0.782300684927837e15 / 0.5994581212784640e16 -0.17511430871269903e17 / 0.26975615457530880e17 -0.413193098349471e15 / 0.1998193737594880e16 0.189367309285289755e18 / 0.194224431294222336e18; 0.894580992211517e15 / 0.182085404338333440e18 -0.209441083772219e15 / 0.27746347327746048e17 -0.4946149632449393e16 / 0.416195209916190720e18 0.334964642443661e15 / 0.2890244513306880e16 -0.604469352802317407e18 / 0.971122156471111680e18 -0.128172128502407e15 / 0.570912496455680e15;]; + +% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0 +Qp_U=[-0.34660470729017653729e20 / 0.113641961250214656000e21 0.351671379135966469961e21 / 0.415604886857927884800e21 0.1819680091728191503e19 / 0.103901221714481971200e21 -0.18252344147469061739e20 / 0.346337405714939904000e21 -0.18145368485798816351e20 / 0.727308552001373798400e21 0.2627410615589536403e19 / 0.138534962285975961600e21; -0.1606450873889019037e19 / 0.7576130750014310400e19 -0.8503979509850519441e19 / 0.9235664152398397440e19 0.2208731907526094393e19 / 0.2308916038099599360e19 0.1143962309827873891e19 / 0.7696386793665331200e19 0.1263616990270014071e19 / 0.16162412266697195520e20 -0.1402288892096389187e19 / 0.27706992457195192320e20; -0.502728075208147729e18 / 0.11364196125021465600e20 0.5831273443201206481e19 / 0.41560488685792788480e20 -0.9031420599281409001e19 / 0.10390122171448197120e20 0.29977986617775158621e20 / 0.34633740571493990400e20 -0.7995649008389734663e19 / 0.72730855200137379840e20 0.728315692435313537e18 / 0.41560488685792788480e20; 0.710308100786581369e18 / 0.11364196125021465600e20 -0.2317346723533341809e19 / 0.41560488685792788480e20 -0.1357359577229545879e19 / 0.10390122171448197120e20 -0.35124499190079631261e20 / 0.34633740571493990400e20 0.81675241511291974823e20 / 0.72730855200137379840e20 0.144034831596315317e18 / 0.13853496228597596160e20; 0.13360631165154733e17 / 0.841792305557145600e18 -0.875389186128426797e18 / 0.27706992457195192320e20 0.95493318392786453e17 / 0.6926748114298798080e19 0.1714625642820850967e19 / 0.23089160380995993600e20 -0.53371483072841696197e20 / 0.48487236800091586560e20 0.30168063964639488547e20 / 0.27706992457195192320e20; -0.1943232250834614071e19 / 0.113641961250214656000e21 0.8999269402554660119e19 / 0.415604886857927884800e21 0.1242786058815312817e19 / 0.103901221714481971200e21 -0.7480218714053301461e19 / 0.346337405714939904000e21 0.28468183824757664191e20 / 0.727308552001373798400e21 -0.476082521721490837529e21 / 0.415604886857927884800e21;]; + +Hp=spdiags(ones(np,1),0,np,np); +Hp(1:bp,1:bp)=Hp_U; +Hp(np-bp+1:np,np-bp+1:np)=fliplr(flipud(Hp_U)); +Hp=Hp*h; +HIp=inv(Hp); + +Hm=spdiags(ones(nm,1),0,nm,nm); +Hm(1:bp,1:bp)=Hm_U; +Hm(nm-bp+1:nm,nm-bp+1:nm)=fliplr(flipud(Hm_U)); +Hm=Hm*h; +HIm=inv(Hm); + +Qpp=spdiags(repmat([-0.157e3 / 0.15360e5 0.537e3 / 0.5120e4 -0.3147e4 / 0.5120e4 -0.231e3 / 0.1024e4 0.999e3 / 0.1024e4 -0.1461e4 / 0.5120e4 0.949e3 / 0.15360e5 -0.33e2 / 0.5120e4;],[np,1]),-3:4,np,np); +Qpp(1:bp,1:bp)=Qpp_U; +Qpp(np-bp+1:np,np-bp+1:np)=flipud( fliplr(Qpp_U ) )'; + +Qpm=-Qpp'; + +Qmp=spdiags(repmat([-0.157e3 / 0.15360e5 0.537e3 / 0.5120e4 -0.3147e4 / 0.5120e4 -0.231e3 / 0.1024e4 0.999e3 / 0.1024e4 -0.1461e4 / 0.5120e4 0.949e3 / 0.15360e5 -0.33e2 / 0.5120e4;],[nm,1]),-3:4,nm,nm); +Qmp(1:bp,1:bp)=Qmp_U; +Qmp(nm-bp+1:nm,nm-bp+1:nm)=flipud( fliplr(Qmp_U ) )'; + +Qmm=-Qmp'; + +Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1; +Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1; + + +Dpp=HIp*(Qpp+1/2*Bpp) ; +Dpm=HIp*(Qpm+1/2*Bpp) ; + + +Dmp=HIm*(Qmp+1/2*Bmp) ; +Dmm=HIm*(Qmm+1/2*Bmp) ; + + +%%% Start with the staggered +Qp=spdiags(repmat([-0.3e1 / 0.640e3 0.25e2 / 0.384e3 -0.75e2 / 0.64e2 0.75e2 / 0.64e2 -0.25e2 / 0.384e3 0.3e1 / 0.640e3],[np,1]),-2:3,np,nm); +Qp(1:bp,1:bp)=Qp_U; +Qp(np-bp+1:np,nm-bp+1:nm)=flipud( fliplr(-Qp_U ) ); +Qm=-Qp'; + +Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1; +Bm=Bp'; + +Dp=HIp*(Qp+1/2*Bp) ; + +Dm=HIm*(Qm+1/2*Bm) ; + +% grids +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d1_staggered_upwind_8.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,82 @@ +function [xp,xm,Hp,Hm,HIp,HIm,Dp,Dm,Dpp,Dpm,Dmp,Dmm] = d1_staggered_upwind_8(n, L) + +assert(n-1 >= 12,'Not enough grid points'); + +np=n; +nm=n+1; +h = L/(n-1); + +bp=8; % Number of boundary points + + +%H- norm +Hm_U=[0.27058177e8 / 0.194594400e9 0 0 0 0 0 0 0; 0 0.378196601e9 / 0.464486400e9 0 0 0 0 0 0; 0 0 0.250683799e9 / 0.232243200e9 0 0 0 0 0; 0 0 0 0.48820787e8 / 0.51609600e8 0 0 0 0; 0 0 0 0 0.23953873e8 / 0.23224320e8 0 0 0; 0 0 0 0 0 0.55018021e8 / 0.55738368e8 0 0; 0 0 0 0 0 0 0.284772253e9 / 0.283852800e9 0; 0 0 0 0 0 0 0 0.6036097241e10 / 0.6038323200e10;]; + +%H+ norm +Hp_U=[0.273925927e9 / 0.928972800e9 0 0 0 0 0 0 0; 0 0.1417489391e10 / 0.928972800e9 0 0 0 0 0 0; 0 0 0.26528603e8 / 0.103219200e9 0 0 0 0 0; 0 0 0 0.66843235e8 / 0.37158912e8 0 0 0 0; 0 0 0 0 0.76542481e8 / 0.185794560e9 0 0 0; 0 0 0 0 0 0.132009637e9 / 0.103219200e9 0 0; 0 0 0 0 0 0 0.857580049e9 / 0.928972800e9 0; 0 0 0 0 0 0 0 0.937663193e9 / 0.928972800e9;]; + +% Upper part of Qpp. Notice that the B matrix is not included here. +% Qpp+Qpp^T=S/2,Qpp+Qpm^T=0 +Qpp_U=[-0.5053e4 / 0.43253760e8 0.299126491231e12 / 0.442810368000e12 -0.562745123e9 / 0.8856207360e10 -0.59616938177e11 / 0.398529331200e12 -0.5489464381e10 / 0.177124147200e12 0.6994759151e10 / 0.88562073600e11 0.712609553e9 / 0.181149696000e12 -0.14167e5 / 0.998400e6; -0.74652297589e11 / 0.110702592000e12 -0.995441e6 / 0.302776320e9 0.662285893e9 / 0.2236416000e10 0.2364992837e10 / 0.5535129600e10 0.8893339297e10 / 0.49816166400e11 -0.2166314077e10 / 0.8945664000e10 -0.78813191e8 / 0.2952069120e10 0.88810243397e11 / 0.1992646656000e13; 0.16939159e8 / 0.276756480e9 -0.2302477691e10 / 0.8200192000e10 -0.197067e6 / 0.9175040e7 0.46740527413e11 / 0.88562073600e11 -0.4827794723e10 / 0.11070259200e11 0.155404199e9 / 0.1230028800e10 0.15110006581e11 / 0.295206912000e12 -0.2491856531e10 / 0.88562073600e11; 0.7568509969e10 / 0.49816166400e11 -0.19762404121e11 / 0.44281036800e11 -0.2554553593e10 / 0.5535129600e10 -0.19550057e8 / 0.302776320e9 0.145135201e9 / 0.210862080e9 0.1596117533e10 / 0.11070259200e11 0.104983477e9 / 0.3558297600e10 -0.1087937837e10 / 0.25303449600e11; 0.5282544031e10 / 0.177124147200e12 -0.131784830977e12 / 0.797058662400e12 0.8310607171e10 / 0.22140518400e11 -0.56309573e8 / 0.105431040e9 -0.36477607e8 / 0.302776320e9 0.57443289107e11 / 0.88562073600e11 -0.220576549e9 / 0.691891200e9 0.39942520697e11 / 0.398529331200e12; -0.1743516779e10 / 0.22140518400e11 0.7784403199e10 / 0.32800768000e11 -0.459036521e9 / 0.4920115200e10 -0.5749179391e10 / 0.22140518400e11 -0.562460513e9 / 0.1383782400e10 -0.1500741e7 / 0.9175040e7 0.70986504791e11 / 0.73801728000e11 -0.930160589e9 / 0.3406233600e10; -0.712609553e9 / 0.181149696000e12 0.180761e6 / 0.6589440e7 -0.18016744831e11 / 0.295206912000e12 0.18651112477e11 / 0.797058662400e12 0.7155507361e10 / 0.44281036800e11 -0.49358401541e11 / 0.73801728000e11 -0.55032223e8 / 0.302776320e9 0.38275518139e11 / 0.40255488000e11; 0.14167e5 / 0.998400e6 -0.88810243397e11 / 0.1992646656000e13 0.650307179e9 / 0.22140518400e11 0.474654619e9 / 0.16102195200e11 -0.7267828861e10 / 0.199264665600e12 0.42227833e8 / 0.425779200e9 -0.71245142351e11 / 0.110702592000e12 -0.7998899e7 / 0.43253760e8;]; + +% Upper part of Qmp. Notice that the B matrix is not included here. +% Qmp+Qmp^T=S/2,Qmp+Qmm^T=0 +Qmp_U=[-0.92025012754706822244637e23 / 0.73350939131274317328275670e26 0.930302337374620084855601123690977e33 / 0.1359398284732080668181467336976000e34 -0.727851704787797291000113057457371e33 / 0.2718796569464161336362934673952000e34 0.7360201789777281105403766444579e31 / 0.90626552315472044545431155798400e32 0.13716157770579047943179985700303e32 / 0.543759313892832267272586934790400e33 -0.553550686983724599329589830113e30 / 0.21750372555713290690903477391616e32 0.39889728754596585193681401053e29 / 0.20596943708061828305779808136000e32 0.4597826214543803453588826224861e31 / 0.2718796569464161336362934673952000e34; -0.460817194517403953445488602736801e33 / 0.679699142366040334090733668488000e33 -0.2251237342833501172927156567e28 / 0.267062619272621870023659683840e30 0.190175037040902421814031988319053e33 / 0.202800676510147232549216572416000e33 -0.280264041304585424221475951435491e33 / 0.1081603608054118573595821719552000e34 -0.2948540748840639854083293613843e31 / 0.81120270604058893019686628966400e32 0.10373234521893519136055741251159e32 / 0.194688649449741343247247909519360e33 -0.391946441371315598560753809131e30 / 0.59488198442976521547770194575360e32 -0.15471439859462263133402977429037e32 / 0.6026077244872946338605292437504000e34; 0.702728186606917922841148808013121e33 / 0.2718796569464161336362934673952000e34 -0.738964213741742941634289388463587e33 / 0.811202706040588930196866289664000e33 -0.6938440904692701484464706797e28 / 0.267062619272621870023659683840e30 0.13587935221144065448123508546711e32 / 0.16900056375845602712434714368000e32 -0.398581023027193812136121833051e30 / 0.3244810824162355720787465158656e31 -0.4390487184318690815376655433561e31 / 0.486721623624353358118119773798400e33 0.39375020562052974775725982794643e32 / 0.5948819844297652154777019457536000e34 -0.337431206994896862016299739723e30 / 0.1054563517852765609255926176563200e34; -0.1626484714285090813572964936931e31 / 0.22656638078868011136357788949600e32 0.247081446636583523913555235517491e33 / 0.1081603608054118573595821719552000e34 -0.98606689532786656855690277893813e32 / 0.135200451006764821699477714944000e33 -0.54412910204947725598668603049e29 / 0.801187857817865610070979051520e30 0.44527159063104584858762191285733e32 / 0.54080180402705928679791085977600e32 -0.4170461618635408764831166558231e31 / 0.18541776138070604118785515192320e32 0.2733053052508781643445069010989e31 / 0.55081665224978260692379809792000e32 -0.62697761224771731704532896739409e32 / 0.7030423452351770728372841177088000e34; -0.16732447706243527499842078942603e32 / 0.543759313892832267272586934790400e33 0.166114324967929244010272382781e30 / 0.2897152521573531893560236748800e31 0.45101984053799299547057123965e29 / 0.811202706040588930196866289664e30 -0.17968772701532140224971787578029e32 / 0.27040090201352964339895542988800e32 -0.48635228792591105226576208151e29 / 0.400593928908932805035489525760e30 0.89480251714879473157102672858403e32 / 0.97344324724870671623623954759680e32 -0.42143863277048975845235502012773e32 / 0.148720496107441303869425486438400e33 0.4389611982769118812600028148913e31 / 0.52728175892638280462796308828160e32; 0.590475930631333482244268674627e30 / 0.21750372555713290690903477391616e32 -0.1702402353852012839500776925027e31 / 0.27812664207105906178178272788480e32 0.5538920825569170589993365115709e31 / 0.121680405906088339529529943449600e33 0.13866428236712167588433125984777e32 / 0.129792432966494228831498606346240e33 -0.16470253483258635888811378530287e32 / 0.24336081181217667905905988689920e32 -0.391812951622286986994545412081e30 / 0.2403563573453596830212937154560e31 0.345131812155051195787190432571781e33 / 0.356929190657859129286621167452160e33 -0.8126635967231661246920267342728147e34 / 0.25309524428466374622142228237516800e35; -0.88170259036818455465427409231e29 / 0.41193887416123656611559616272000e32 0.943459254925453305858100463659e30 / 0.118976396885953043095540389150720e33 -0.104249970715879362499981962934393e33 / 0.5948819844297652154777019457536000e34 0.78160374875076232344526852587e29 / 0.18360555074992753564126603264000e32 0.37536035579237342566706229296521e32 / 0.297440992214882607738850972876800e33 -0.240812609144054591374890622842541e33 / 0.356929190657859129286621167452160e33 -0.145390363516271219220036634153e30 / 0.801187857817865610070979051520e30 0.38093040928245760105862644889779897e35 / 0.38667328987934739006050626473984000e35; -0.4596580943680048141920105029111e31 / 0.2718796569464161336362934673952000e34 0.15262190991038620305994595494037e32 / 0.6026077244872946338605292437504000e34 0.1782508142749448850000523441843e31 / 0.1054563517852765609255926176563200e34 -0.33510338065544202178903991034341e32 / 0.7030423452351770728372841177088000e34 -0.8226959851063633384855147175189e31 / 0.421825407141106243702370470625280e33 0.3729744974003821244717412110773747e34 / 0.25309524428466374622142228237516800e35 -0.6554525408845613610278033864711443e34 / 0.9666832246983684751512656618496000e34 -0.49383212966905764275823531781e29 / 0.267062619272621870023659683840e30;]; + +% The staggered + operator, upper part. Notice that the B matrix is not included here.Qp+Qm^T=0 +Qp_U=[-0.6661444046602902130086192779621746771e37 / 0.23342839720855518078613888837079040000e38 0.12367666586033683088530161144944922123649e41 / 0.14797137255429936031548004199961722880000e41 0.4395151478839780503265910147432452357e37 / 0.186518536833150454179176523528929280000e39 -0.87963490365651280757669626389462038003e38 / 0.1345194295948176002868000381814702080000e40 -0.224187853266917333808369723332131873619e39 / 0.5178998039400477611041801469986603008000e40 0.55314264508663245255306547445101748003e38 / 0.2959427451085987206309600839992344576000e40 0.232303691019191510116997121080032335473e39 / 0.7398568627714968015774002099980861440000e40 -0.2816079756494683118054172509551114287e37 / 0.182680706857159704093185237036564480000e39; -0.10353150194859080046224306922028068797e38 / 0.43350988053017390717425793554575360000e38 -0.41810120772068966379630018024826564868789e41 / 0.44391411766289808094644012599885168640000e41 0.14548989645937048285445132077600344947e38 / 0.16118885899161150361163403267932160000e38 0.1353060140543056426043612771323929599939e40 / 0.6341630252327115442092001799983595520000e40 0.52584560489710238297754746193308329991e38 / 0.317081512616355772104600089999179776000e39 -0.322880558646947920788817445027309369183e39 / 0.8878282353257961618928802519977033728000e40 -0.2523050363123545986203815545725359199053e40 / 0.22195705883144904047322006299942584320000e41 0.2171073153815036601208579021528549178587e40 / 0.44391411766289808094644012599885168640000e41; -0.891108828407068615176254357157902263e36 / 0.14450329351005796905808597851525120000e38 0.9207480733237809022786174965361787199767e40 / 0.44391411766289808094644012599885168640000e41 -0.42541158717297365539935073107511004261e38 / 0.62172845611050151393058841176309760000e38 0.3752589973871193536179209598104238527889e40 / 0.4932379085143312010516001399987240960000e40 -0.503978909830206767329311700636686423011e39 / 0.2219570588314490404732200629994258432000e40 -0.25729240076595732139825342259542295873e38 / 0.328825272342887467367733426665816064000e39 0.283553624900223723890865927423934504751e39 / 0.2466189542571656005258000699993620480000e40 -0.129077863120107719239717930927898530811e39 / 0.4035582887844528008604001145444106240000e40; 0.150348887575956035991597139943595613e36 / 0.2000814833216187263881190471749632000e37 -0.20271017027971837857434511014725441049e38 / 0.422775350155141029472800119998906368000e39 -0.28593401740189393456434271848012721991e38 / 0.87041983855470211950282377646833664000e38 -0.3107901953269564253446821657614993689969e40 / 0.2959427451085987206309600839992344576000e40 0.163990240717967340033770840655582198979e39 / 0.147971372554299360315480041999617228800e39 0.1282390462809880153203810439898877805771e40 / 0.5326969411954776971357281511986220236800e40 0.55264142141446425784998699697993554089e38 / 0.1479713725542993603154800419996172288000e40 -0.11464469241412665723941542481946601159e38 / 0.328825272342887467367733426665816064000e39; 0.523088046329055388999216632374482217e36 / 0.8670197610603478143485158710915072000e37 -0.1060078188350823496391913491020744512571e40 / 0.8878282353257961618928802519977033728000e40 0.2489758378222482566046953325470732791e37 / 0.87041983855470211950282377646833664000e38 0.2223322627542429729644405110485179851147e40 / 0.8878282353257961618928802519977033728000e40 -0.405229856535104846314030690968355556777e39 / 0.443914117662898080946440125998851686400e39 0.1490365162783652350291746570529147907183e40 / 0.1775656470651592323785760503995406745600e40 -0.768972334347761129534143069144613163467e39 / 0.4439141176628980809464401259988516864000e40 0.34866718651627637034229483614924933299e38 / 0.1268326050465423088418400359996719104000e40; -0.522920230014765612739540250860177277e36 / 0.14450329351005796905808597851525120000e38 0.1453599225900263767517831305438030179313e40 / 0.44391411766289808094644012599885168640000e41 0.39466073050115777575909309244731209107e38 / 0.435209919277351059751411888234168320000e39 -0.501532616891797921922210650058733239369e39 / 0.4932379085143312010516001399987240960000e40 -0.276688880738967796756224778348832938429e39 / 0.2219570588314490404732200629994258432000e40 -0.346414321340583244146983258114429082007e39 / 0.328825272342887467367733426665816064000e39 0.12835585968364300662881260070481636919e38 / 0.10676145205937904784666669696942080000e38 -0.823464330534329517579044232163467356639e39 / 0.44391411766289808094644012599885168640000e41; -0.52802811745049309885720368328150617e35 / 0.1688999534533145092886719229399040000e37 0.306081878756402097710283710502974342821e39 / 0.4932379085143312010516001399987240960000e40 -0.43281676756908448328630499995327205907e38 / 0.1305629757832053179254235664702504960000e40 -0.136777517026276610492513752852976312597e39 / 0.4932379085143312010516001399987240960000e40 0.7290267185112688983569547959866321207e37 / 0.246618954257165600525800069999362048000e39 0.351639553804386998900032061955116226307e39 / 0.3804978151396269265255201079990157312000e40 -0.2840607921448362125113163437018746783083e40 / 0.2466189542571656005258000699993620480000e40 0.1859197240206073478058984196606793676119e40 / 0.1644126361714437336838667133329080320000e40; 0.5412884950805861225701675068383948563e37 / 0.303456916371121735021980554882027520000e39 -0.1279848124286614098666833963684894093627e40 / 0.44391411766289808094644012599885168640000e41 0.17218770500911534891873213313315117e35 / 0.39564538116122823613764717112197120000e38 0.904771800018341402297740826650985365019e39 / 0.44391411766289808094644012599885168640000e41 0.65377500469171784914136353007658339737e38 / 0.15536994118201432833125404409959809024000e41 -0.211014619053462658139013021424372225649e39 / 0.8878282353257961618928802519977033728000e40 0.195325945159072852492812627815477159003e39 / 0.3170815126163557721046000899991797760000e40 -0.52260858238454846311625894178508086247499e41 / 0.44391411766289808094644012599885168640000e41;]; + +Hp=spdiags(ones(np,1),0,np,np); +Hp(1:bp,1:bp)=Hp_U; +Hp(np-bp+1:np,np-bp+1:np)=fliplr(flipud(Hp_U)); +Hp=Hp*h; +HIp=inv(Hp); + +Hm=spdiags(ones(nm,1),0,nm,nm); +Hm(1:bp,1:bp)=Hm_U; +Hm(nm-bp+1:nm,nm-bp+1:nm)=fliplr(flipud(Hm_U)); +Hm=Hm*h; +HIm=inv(Hm); + +tt=[0.1447e4 / 0.688128e6 -0.17119e5 / 0.688128e6 0.8437e4 / 0.57344e5 -0.5543e4 / 0.8192e4 -0.15159e5 / 0.81920e5 0.16139e5 / 0.16384e5 -0.2649e4 / 0.8192e4 0.15649e5 / 0.172032e6 -0.3851e4 / 0.229376e6 0.5053e4 / 0.3440640e7;]; + +Qpp=spdiags(repmat(tt,[np,1]),-4:5,np,np); +Qpp(1:bp,1:bp)=Qpp_U; +Qpp(np-bp+1:np,np-bp+1:np)=flipud( fliplr(Qpp_U ) )'; + +Qpm=-Qpp'; + +Qmp=spdiags(repmat(tt,[nm,1]),-4:5,nm,nm); +Qmp(1:bp,1:bp)=Qmp_U; +Qmp(nm-bp+1:nm,nm-bp+1:nm)=flipud( fliplr(Qmp_U ) )'; + +Qmm=-Qmp'; + + +Bpp=spalloc(np,np,2);Bpp(1,1)=-1;Bpp(np,np)=1; +Bmp=spalloc(nm,nm,2);Bmp(1,1)=-1;Bmp(nm,nm)=1; + +Dpp=HIp*(Qpp+1/2*Bpp) ; +Dpm=HIp*(Qpm+1/2*Bpp) ; + + +Dmp=HIm*(Qmp+1/2*Bmp) ; +Dmm=HIm*(Qmm+1/2*Bmp) ; + + +%%% Start with the staggered +Qp=spdiags(repmat([0.5e1 / 0.7168e4 -0.49e2 / 0.5120e4 0.245e3 / 0.3072e4 -0.1225e4 / 0.1024e4 0.1225e4 / 0.1024e4 -0.245e3 / 0.3072e4 0.49e2 / 0.5120e4 -0.5e1 / 0.7168e4;],[np,1]),-3:4,np,nm); +Qp(1:bp,1:bp)=Qp_U; +Qp(np-bp+1:np,nm-bp+1:nm)=flipud( fliplr(-Qp_U ) ); +Qm=-Qp'; + +Bp=spalloc(np,nm,2);Bp(1,1)=-1;Bp(np,nm)=1; +Bm=Bp'; + +Dp=HIp*(Qp+1/2*Bp) ; + +Dm=HIm*(Qm+1/2*Bm) ; + +% grids +xp = h*[0:n]'; +xm = h*[0 1/2+0:n n]'; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_2.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,50 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<3) + error(['Operator requires at least ' num2str(3) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + % D1 operator + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + scheme_width = 3; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + Mm1 = -c(r-1)/2 - c(r)/2; + M0 = c(r-1)/2 + c(r) + c(r+1)/2; + Mp1 = -c(r)/2 - c(r+1)/2; + + vals = [Mm1,M0,Mp1]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_4.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,57 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_4(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<5) + error(['Operator requires at least ' num2str(5) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + S = d1_l*d1_l' + d1_r*d1_r'; + + % D1 operator + stencil = [1/12 -2/3 0 2/3 -1/12]; + diags = -2:2; + Q = stripeMatrixPeriodic(stencil, diags, m); + D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); + + + scheme_width = 5; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + % Note: these coefficients are for -M. + Mm2 = -1/8*c(r-2) + 1/6*c(r-1) - 1/8*c(r); + Mm1 = 1/6 *c(r-2) + 1/2*c(r-1) + 1/2*c(r) + 1/6*c(r+1); + M0 = -1/24*c(r-2)- 5/6*c(r-1) - 3/4*c(r) - 5/6*c(r+1) - 1/24*c(r+2); + Mp1 = 0 * c(r-2) + 1/6*c(r-1) + 1/2*c(r) + 1/2*c(r+1) + 1/6 *c(r+2); + Mp2 = 0 * c(r-2) + 0 * c(r-1) - 1/8*c(r) + 1/6*c(r+1) - 1/8 *c(r+2); + + vals = -[Mm2,Mm1,M0,Mp1,Mp2]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + + end + D2 = @D2_fun; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_6.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,58 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_6(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<7) + error(['Operator requires at least ' num2str(7) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + + % D1 operator + diags = -3:3; + stencil = [-1/60 9/60 -45/60 0 45/60 -9/60 1/60]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + % D2 operator + scheme_width = 7; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + Mm3 = c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r); + Mm2 = c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2); + Mm1 = -c(r-3)/0.40e2 - 0.3e1/0.10e2 * c(r-2) - 0.3e1/0.10e2 * c(r+1) - c(r+2)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r-1); + M0 = c(r-3)/0.180e3 + c(r-2)/0.8e1 + 0.19e2/0.20e2 * c(r-1) + 0.19e2/0.20e2 * c(r+1) + c(r+2)/0.8e1 + c(r+3)/0.180e3 + 0.101e3/0.180e3 * c(r); + Mp1 = -c(r-2)/0.40e2 - 0.3e1/0.10e2 * c(r-1) - 0.3e1/0.10e2 * c(r+2) - c(r+3)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r+1); + Mp2 = c(r-1)/0.20e2 - 0.3e1/0.10e2 * c(r+1) + c(r+3)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r+2); + Mp3 = c(r+1)/0.40e2 + c(r+2)/0.40e2 - 0.11e2/0.360e3 * c(r) - 0.11e2/0.360e3 * c(r+3); + + vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; + + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_2to2_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,17 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2 +%INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F1_ACCF2C2_STENCIL_5_BC_2_5 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F1_ACCF2C2_STENCIL_5_BC_2_5 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:07 + +stencil_F2C = [-1.0./8.0,1.0./4.0,3.0./4.0,1.0./4.0,-1.0./8.0]; +if nargout > 1 + BC_F2C = reshape([6.288191560529559e-1,-6.440957802647795e-2,6.855471317807184e-1,1.572264341096408e-1,-2.438028498638558e-1,7.469014249319279e-1,-8.431231982626708e-2,2.921561599131335e-1,1.374888185644862e-2,-1.318744409282243e-1],[2,5]); +end +if nargout > 2 + HcU = 1.0./2.0; +end +if nargout > 3 + HfU = 1.0./2.0; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,17 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1 +%INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F2_ACCF2C1_STENCIL_5_BC_1_3 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_2TO2_RATIO_2TO1_ACCC2F2_ACCF2C1_STENCIL_5_BC_1_3 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:08 + +stencil_F2C = [1.0./4.0,1.0./2.0,1.0./4.0]; +if nargout > 1 + BC_F2C = [1.0./2.0,1.0./2.0]; +end +if nargout > 2 + HcU = 1.0./2.0; +end +if nargout > 3 + HfU = 1.0./2.0; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_4to4_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3 +%INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F2_ACCF2C3_STENCIL_9_BC_3_11 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F2_ACCF2C3_STENCIL_9_BC_3_11 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:01 + +stencil_F2C = [7.0./2.56e2,-1.0./3.2e1,-7.0./6.4e1,9.0./3.2e1,8.5e1./1.28e2,9.0./3.2e1,-7.0./6.4e1,-1.0./3.2e1,7.0./2.56e2]; +if nargout > 1 + BC_F2C = reshape([7.523257802630956e-2,2.447812262221267e-1,-1.679313063616916e-1,1.290510950666589,6.315723344677289e-3,1.67178747937954e-1,1.982667903557025,-7.554383893379468e-1,7.215271362899867e-1,-2.820478831383137,1.807034411305536,-7.589683751979843e-1,-7.685973268458095e-1,3.965751544535173e-1,4.119789638051451e-1,1.574000556785898e-1,-1.113618964927466e-1,3.631000220124657e-1,1.639694219982991,-9.114074742274862e-1,4.952715520862987e-1,-4.524162151353456e-1,2.65481378213589e-1,-2.268273408674622e-1,3.455365956773523e-1,-2.009765975089827e-1,1.722179564305811e-1,-5.52933488235368e-1,3.18109976271229e-1,-2.171481232558432e-1,1.033835580108027e-1,-5.911351224351343e-2,3.960076712054991e-2],[3,11]); +end +if nargout > 2 + t2 = [1.7e1./4.8e1,5.9e1./4.8e1,4.3e1./4.8e1,4.9e1./4.8e1]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2 +%INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F3_ACCF2C2_STENCIL_9_BC_3_11 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_4TO4_RATIO_2TO1_ACCC2F3_ACCF2C2_STENCIL_9_BC_3_11 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:36:05 + +stencil_F2C = [-1.0./2.56e2,-1.0./3.2e1,1.0./6.4e1,9.0./3.2e1,6.1e1./1.28e2,9.0./3.2e1,1.0./6.4e1,-1.0./3.2e1,-1.0./2.56e2]; +if nargout > 1 + BC_F2C = reshape([1.0./2.0,0.0,0.0,1.77e2./2.72e2,3.0./8.0,-5.9e1./6.88e2,1.125919117647059e-2,3.546742584745763e-1,1.335392441860465e-2,-4.9e1./5.44e2,2.335805084745763e-1,3.204941860465116e-1,-1.194852941176471e-2,1.350635593220339e-2,5.308866279069767e-1,-9.0./5.44e2,-1.11228813559322e-2,2.943313953488372e-1,-2.803308823529412e-2,2.105402542372881e-2,-1.580668604651163e-2,-9.0./5.44e2,1.430084745762712e-2,-5.450581395348837e-2,-9.191176470588235e-4,7.944915254237288e-4,-5.450581395348837e-3,1.0./5.44e2,-1.588983050847458e-3,2.180232558139535e-3,2.297794117647059e-4,-1.986228813559322e-4,2.725290697674419e-4],[3,11]); +end +if nargout > 2 + t2 = [1.7e1./4.8e1,5.9e1./4.8e1,4.3e1./4.8e1,4.9e1./4.8e1]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_6to6_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4 +%INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F3_ACCF2C4_STENCIL_13_BC_3_17 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F3_ACCF2C4_STENCIL_13_BC_3_17 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:47 + +stencil_F2C = [-5.95703125e-3,3.0./5.12e2,3.57421875e-2,-2.5e1./5.12e2,-8.935546875e-2,7.5e1./2.56e2,3.17e2./5.12e2,7.5e1./2.56e2,-8.935546875e-2,-2.5e1./5.12e2,3.57421875e-2,3.0./5.12e2,-5.95703125e-3]; +if nargout > 1 + BC_F2C = reshape([5.233890618131365e-1,-1.594459192645467e-2,3.532688727637403e-2,8.021234957689208e-1,3.90683205173562e-1,-1.73222951632239e-1,-8.87662686483442e-2,2.90091235796637e-1,-1.600356115709148e-1,-1.044025375027475e-1,2.346179009198368e-1,6.091329306528956e-1,1.561275522703128e-1,-1.168382445709856e-1,1.040987887887311,-2.387061036980731e-1,1.363504965974361e-1,1.082611928255256e-1,-5.745310654054326e-1,3.977694249198785e-1,-9.376217911402619e-1,3.554518646054656e-2,5.609157787396987e-5,-1.564625311232018e-1,6.107907974027401e-1,-3.786608696698368e-1,7.02265869951125e-1,2.054294270642538e-1,-1.302300112378257e-1,2.478941407690889e-1,-2.657085326191479e-1,1.568761445933572e-1,-2.632906518005349e-1,-1.228047556139644e-1,7.182193248980271e-2,-1.1291238242346e-1,5.811258780405158e-2,-3.466364400805378e-2,5.683680203338252e-2,1.781337575097077e-2,-1.030572999042704e-2,1.577197067502767e-2,-1.443802115768554e-2,8.391316308374261e-3,-1.29534117369052e-2,-1.548018627738296e-3,8.794183904936319e-4,-1.298961407229805e-3,1.573818938200601e-3,-8.940753636685258e-4,1.320610764016968e-3],[3,17]); +end +if nargout > 2 + t2 = [3.159490740740741e-1,1.390393518518519,6.275462962962963e-1,1.240509259259259,9.116898148148148e-1,1.013912037037037]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3 +%INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F4_ACCF2C3_STENCIL_13_BC_4_17 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_6TO6_RATIO_2TO1_ACCC2F4_ACCF2C3_STENCIL_13_BC_4_17 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:56 + +stencil_F2C = [1.07421875e-3,3.0./5.12e2,-6.4453125e-3,-2.5e1./5.12e2,1.611328125e-2,7.5e1./2.56e2,2.45e2./5.12e2,7.5e1./2.56e2,1.611328125e-2,-2.5e1./5.12e2,-6.4453125e-3,3.0./5.12e2,1.07421875e-3]; +if nargout > 1 + BC_F2C = reshape([1.0./2.0,0.0,0.0,0.0,6.876076086160158e-1,1.5e1./3.2e1,-3.461879841386942e-1,3.502577439820862e-2,3.099721991995751e-3,2.228547003766753e-1,9.363652999354482e-3,-3.15791046603844e-3,-1.05789143206462e-1,2.35563165945226e-1,6.070385985337514e-1,-4.847496617839149e-2,-4.809232246867902e-3,5.154694068405061e-3,7.049223649022501e-1,1.016749349808733e-2,3.460117842882262e-2,-4.996621498584866e-2,5.210650763094799e-1,2.233592875303228e-1,-2.239024779745769e-3,4.254668119953384e-4,9.213872590833641e-3,3.910524351558127e-1,-9.048711215107334e-2,8.597375629526346e-2,-3.468219752858724e-1,3.250682992395969e-1,-1.309107466664224e-1,1.199249722306043e-1,-4.071552384959425e-1,1.474417123938701e-1,-3.02863877390285e-2,3.046016554982103e-2,-1.081315128181483e-1,1.120705238850532e-2,7.977722870036266e-2,-6.772545887788229e-2,2.00270418837606e-1,-5.427183680490763e-2,6.524845570188292e-2,-5.637610037043203e-2,1.711235764016968e-1,-4.203646902407166e-2,4.639548341453586e-3,-4.055884445496545e-3,1.258630924935448e-2,-2.776451559292779e-3,-1.023784366070774e-2,8.817108288936985e-3,-2.659664906860937e-2,7.14445034054861e-3,-1.435466025441424e-3,1.240274208149505e-3,-3.764718680837329e-3,1.028716295017727e-3,1.032012418492197e-3,-8.794183904936319e-4,2.597922814459609e-3,-6.571159498040679e-4,1.892022767235695e-4,-1.612267049238325e-4,4.76285849317595e-4,-1.204712574640791e-4],[4,17]); +end +if nargout > 2 + t2 = [3.159490740740741e-1,1.390393518518519,6.275462962962963e-1,1.240509259259259,9.116898148148148e-1,1.013912037037037]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,47 @@ +function [IC2F,IF2C,Hc,Hf] = IntOp_orders_8to8_ratio2to1(mc,hc,ACC) + +% ACC is a string. +% ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. +% ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. +ratio = 2; +mf = ratio*mc-1; +hf = hc/ratio; + +switch ACC + case 'F2C' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5; + case 'C2F' + [stencil_F2C,BC_F2C,HcU,HfU] = ... + sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4; +end + +stencil_width = length(stencil_F2C); +stencil_hw = (stencil_width-1)/2; +[BC_rows,BC_cols] = size(BC_F2C); + +%%% Norm matrices %%% +Hc = speye(mc,mc); +HcUm = length(HcU); +Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); +Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); +Hc = Hc*hc; + +Hf = speye(mf,mf); +HfUm = length(HfU); +Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); +Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); +Hf = Hf*hf; +%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IF2C from stencil and BC +IF2C = sparse(mc,mf); +for i = BC_rows+1 : mc-BC_rows + IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> +end +IF2C(1:BC_rows,1:BC_cols) = BC_F2C; +IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%% Create IC2F using symmetry condition %%%% +IC2F = Hf\IF2C.'*Hc; \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5 +%INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F4_ACCF2C5_STENCIL_17_BC_5_24 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F4_ACCF2C5_STENCIL_17_BC_5_24 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:02 + +stencil_F2C = [1.324608212425595e-3,-1.220703125e-3,-1.059686569940476e-2,1.1962890625e-2,3.708902994791667e-2,-5.9814453125e-2,-7.417805989583333e-2,2.99072265625e-1,5.927225748697917e-1,2.99072265625e-1,-7.417805989583333e-2,-5.9814453125e-2,3.708902994791667e-2,1.1962890625e-2,-1.059686569940476e-2,-1.220703125e-3,1.324608212425595e-3]; +if nargout > 1 + BC_F2C = reshape([2.087365925359584,-1.227221822236823,1.090916714284468e1,-1.041312149906314,1.134214373258162,-1.269686811479259,2.075368250436277,-1.520771409696474e1,1.389750473974189,-1.48485748783569,-1.040579640776707,8.899721508624742e-1,-7.177691661032869,6.882712043721976e-1,-7.599347599660409e-1,-1.26440526850463,1.160658043364166,-5.391509178445742,6.679923135342851e-1,-7.521725463922262e-1,3.752349810676949,-2.906684049793639,2.672876913566556e1,-2.493145660873,2.782825829667436,1.443588084644871e-1,-1.307230271348211e-1,2.216804748506809,1.437719804483573e-1,-1.145830952112369e-1,3.2149320801083e-1,-2.34361439272989e-1,1.855837403850803,1.286992417665136e-1,-5.758238484341108e-2,-1.301103772185027,9.973408313121463e-1,-8.837520107531879,9.534930382604288e-1,-1.429199518808459e-2,-2.384453600787036,1.818345997558567,-1.577342441239303e1,1.422846302346762,3.154159322410927e-3,3.589377915408905e-1,-2.109524979864723e-1,9.363227187483732e-1,6.301238593470628e-2,5.732390257964316e-2,2.793582225299658,-2.033528056927806,1.621324028555783e1,-1.189222718426265,6.189671387692164e-2,5.342131492864642e-1,-4.222231296997605e-1,3.787253511209165,-3.350064226195165e-1,1.245325350855226e-2,-2.184146687275183,1.551331905163865,-1.19387312118926e1,8.349589004472998e-1,-1.35383591594438e-1,-8.014367493869704e-1,5.668415377001197e-1,-4.302842643863972,2.84183319528176e-1,3.83514596570461e-2,9.45794880458861e-1,-6.732638898386457e-1,5.234390273661413,-3.83498639699941e-1,1.440453495443018e-1,5.645065377027613e-1,-4.039220855868618e-1,3.170156487564899,-2.383074337879712e-1,1.216617375214008e-1,-1.937863224145611e-1,1.358891032584724e-1,-1.023378722711107,6.830621267493578e-2,-4.169135139842575e-3,2.270306426219975e-2,-2.38399209865252e-2,2.928388757260058e-1,-4.02060181933195e-2,6.880595419697505e-2,9.638828970166005e-2,-6.98778025332229e-2,5.609317796857676e-1,-4.429526993510467e-2,2.882554468024363e-2,3.582431391759518e-2,-2.671910644265286e-2,2.251335440088556e-1,-1.966486325944329e-2,1.791756577091828e-2,-3.350790201878844e-1,2.581411034605136e-1,-2.283079071439122,2.162064132948073e-1,-2.321676317817421e-1,6.356644993195555e-2,-4.921907141164279e-2,4.384534411523264e-1,-4.198474091936761e-2,4.597203884996652e-2,2.261662512481835e-2,-1.740429407604355e-2,1.537038474929879e-1,-1.452709057733876e-2,1.556101844926456e-2,3.097679325854209e-2,-2.394872918869654e-2,2.128879105995857e-1,-2.032077838507769e-2,2.213372706946953e-2],[5,24]); +end +if nargout > 2 + t2 = [2.948906761778786e-1,1.525720623897707,2.57452876984127e-1,1.798113701499118,4.127080577601411e-1,1.278484623015873,9.232955798059965e-1,1.009333860859158]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,18 @@ +function [stencil_F2C,BC_F2C,HcU,HfU] = intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4 +%INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F5_ACCF2C4_STENCIL_17_BC_6_24 +% [STENCIL_F2C,BC_F2C,HCU,HFU] = INT_ORDERS_8TO8_RATIO_2TO1_ACCC2F5_ACCF2C4_STENCIL_17_BC_6_24 + +% This function was generated by the Symbolic Math Toolbox version 8.0. +% 21-May-2018 15:35:39 + +stencil_F2C = [-2.564929780505952e-4,-1.220703125e-3,2.051943824404762e-3,1.1962890625e-2,-7.181803385416667e-3,-5.9814453125e-2,1.436360677083333e-2,2.99072265625e-1,4.820454915364583e-1,2.99072265625e-1,1.436360677083333e-2,-5.9814453125e-2,-7.181803385416667e-3,1.1962890625e-2,2.051943824404762e-3,-1.220703125e-3,-2.564929780505952e-4]; +if nargout > 1 + BC_F2C = reshape([-1.673327822994457e1,1.665420585653232e1,-1.973927473454954e2,2.826257908914756e1,-6.156813484052936e1,3.974966126696054,2.880639373222355e1,-2.660797294399895e1,3.202303847857831e2,-4.598960974033412e1,1.003152507870138e2,-6.481221721577338,8.061874893005659,-7.706608975219419,9.234201483040626e1,-1.322147613066874e1,2.880209985918895e1,-1.859523138300886,-2.19528557898034e1,2.13763239391502e1,-2.476320412310993e2,3.572924955195535e1,-7.795290920161567e1,5.036097395109802,-1.287197329845376,1.244100160910233,-1.394684346537857e1,2.112320869857597,-4.60441856527933,2.978264879234002e-1,1.935199269380359,-1.88580031845933,2.330807783420723e1,-2.91749306966336,6.644427567385864,-4.302395355488564e-1,-1.148984264037686,1.109982778464842,-1.314799057941425e1,2.137017065505111,-4.084082382904703,2.606486561845403e-1,1.565462612018858,-1.508309105065004,1.772442871757061e1,-2.40060893717649,6.31422909250156,-4.052444871368164e-1,3.179380137695096e-1,-3.071741522636364e-1,3.636264382483117,-5.183996466001354e-1,2.322192602480252,-6.467931586900172e-2,6.277081831671831e-1,-6.107890751034906e-1,7.335850436798604,-1.091016872160133,3.089218691662174,6.988892111841499e-2,1.639756762532723,-1.583883992641009,1.876236237708381e1,-2.685390559486947,5.859855615166555,3.919312696253764e-3,1.611803934956753,-1.540778869185931,1.795955382224485e1,-2.495154500270907,5.003492494640851,-4.157919873785453e-2,-6.213413520816026e-1,6.242581523412814e-1,-7.823535324240331,1.222418137438074,-3.120142966701915,2.975831327541895e-1,2.12230710100161,-2.0451058741188,2.412395225004063e1,-3.423544496956294,7.325914828631551,-4.793445914188063e-1,-4.31809447483846,4.158975059021769,-4.906211424150023e1,6.975058620124194,-1.501259943091044e1,9.418983630154896e-1,-6.138797742037704e-1,5.810969368378602e-1,-6.6821474226667,9.113864351863538e-1,-1.814400062103318,1.002574935832871e-1,-9.860044240994036e-2,9.448685284232805e-2,-1.106188847377502,1.552858946803229e-1,-3.265426230113952e-1,2.062042480203908e-2,-2.745034502159738,2.655226942873441,-3.151193386341271e1,4.520932327069998,-9.883338555855938,6.423439846307804e-1,2.139968493382333,-2.067738462547628,2.450223844197562e1,-3.50699574523559,7.635005983192545,-4.923960464980415e-1,1.036406872394752,-1.002008838408551,1.188339891704641e1,-1.70302259679446,3.715789696036691,-2.407301455417931e-1,1.606428513214277,-1.552524985439685,1.840245870449254e1,-2.635131753653705,5.741508438443495,-3.708346905830083e-1,5.659562678739624e-2,-5.465656201303744e-2,6.471932864664631e-1,-9.253169145609999e-2,2.010909570941804e-1,-1.292046400706276e-2,5.435391241988039e-1,-5.252672844681129e-1,6.225561585258754,-8.91344337215917e-1,1.941632399616401,-1.253426955105431e-1,-1.552116972709218,1.499962759958308,-1.777819805127299e1,2.545472086708348,-5.545140384142844,3.580057322157566e-1],[6,24]); +end +if nargout > 2 + t2 = [2.948906761778786e-1,1.525720623897707,2.57452876984127e-1,1.798113701499118,4.127080577601411e-1,1.278484623015873,9.232955798059965e-1,1.009333860859158]; + HcU = t2; +end +if nargout > 3 + HfU = t2; +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_2to2_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,35 @@ +function [IC2F,IF2C] = intOpMC_orders_2to2_ratio2to1(mc) + +mf = 2*(mc-1) + 1; + +stencil_F2C = [1.0./4.0,1.0./2.0,1.0./4.0]; +stencil_width = length(stencil_F2C); +stencil_halfwidth = (stencil_width-1)/2; + +BC_F2C = [1.0./2.0,1.0./2.0]; + +Hc = speye(mc,mc); +Hc(1,1) = 1/2; +Hc(end,end) = 1/2; + +Hf = 1/2*speye(mf,mf); +Hf(1,1) = 1/4; +Hf(end,end) = 1/4; + +IF2C = sparse(mc,mf); +[BCrows, BCcols] = size(BC_F2C); +IF2C(1:BCrows, 1:BCcols) = BC_F2C; +IF2C(mc-BCrows+1:mc, mf-BCcols+1:mf) = rot90(BC_F2C,2); + +for i = BCrows+1 : mc-BCrows + IF2C(i,(2*i-stencil_halfwidth-1) :(2*i+stencil_halfwidth-1))... + = stencil_F2C; +end + +IC2F = Hf\(IF2C'*Hc); + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_4to4_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,66 @@ +% Marks New interpolation operators +% 4th order to 2nd order accurate (diagonal norm) +% M=9 is the minimum amount of points on the coarse mesh + +function [IC2F,IF2C] = intOpMC_orders_4to4_ratio2to1(M_C) + +M_F=M_C*2-1; + +% Coarse to fine +I1=zeros(M_F,M_C); +t1=[ 2047/2176 , 129/1088 , -129/2176 , 0 , 0 , 0 , 0 ,... + 0 ;... + 429/944 , 279/472 , -43/944 , 0 , 0 , 0 , 0 , 0 ;... + 111/2752 , 4913/5504 , 3/32 , -147/5504 ,... + 0 , 0 , 0 , 0 ;... + -103/784 , 549/784 , 387/784 , -1/16 , 0 ,... + 0 , 0 , 0 ;... + -335/3072 , 205/768 , 2365/3072 , 49/512 ,... + -3/128 , 0 , 0 , 0 ;... + -9/256 , 5/256 , 129/256 , 147/256 ,... + -1/16 , 0 , 0 , 0 ;... + 5/192 , -59/1024 , 43/512 , 2695/3072 ,... + 3/32 , -3/128 , 0 , 0 ;... + 23/768 , -37/768 , -43/768 , 147/256 ,... + 9/16 , -1/16 , 0 , 0 ;... + 13/2048 , -11/1024 , -43/2048 , 49/512 ,... + 55/64 , 3/32 , -3/128 , 0 ;... + -1/384 , 1/256 , 0 , -49/768 , 9/16 ,... + 9/16 , -1/16 , 0 ;... + -1/1024 , 3/2048 , 0 , -49/2048 , 3/32 ,... + 55/64 , 3/32 , -3/128]; + + +t2=[-3/128 3/32 55/64 3/32 -3/128]; +t3=[-1/16 9/16 9/16 -1/16]; +I1(1:11,1:8)=t1; +I1(M_F-10:M_F,M_C-7:M_C)=fliplr(flipud(t1)); +I1(12,5:8)=t3; +for i=13:2:M_F-12 + j=(i-1)/2; + I1(i,j-1:j+3)=t2; + I1(i+1,j:j+3)=t3; +end + +% Fine to coarse +I2=zeros(M_C,M_F); + +t1=[ 2047/4352 , 429/544 , 111/2176 , -103/544 ... + , -335/2176 , -27/544 , 5/136 , 23/544 ,... + 39/4352 , -1/272 , -3/2176 ;... + 129/7552 , 279/944 , 4913/15104 , 549/1888 ... + , 205/1888 , 15/1888 , -3/128 , -37/1888 ... + , -33/7552 , 3/1888 , 9/15104 ]; +t2=[-3/256 -1/32 3/64 9/32 55/128 9/32 3/64 -1/32 -3/256]; + +I2(1:2,1:11)=t1; + +I2(M_C-1:M_C,M_F-10:M_F)=fliplr(flipud(t1)); + +for i=3:M_C-2 + j=2*(i-3)+1; + I2(i,j:j+8)=t2; +end + +IC2F = sparse(I1); +IF2C = sparse(I2);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_6to6_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,115 @@ +% Marks New interpolation operators +% 6th order accurate (diagonal norm) +% M=19 is the minimum amount of points on the coarse mesh + +function [IC2F,IF2C] = intOpMC_orders_6to6_ratio2to1(M_C) + +M_F=M_C*2-1; + +% Coarse to fine +I1=zeros(M_F,M_C); + +t1= [6854313/6988288 , 401925/6988288 , -401925/6988288 ... + , 133975/6988288 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + 560547/1537664 , 1201479/1537664 , -240439/1537664 ... + , 16077/1537664 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + 203385/5552128 , 1225647/1388032 , 364155/2776064 ... + , -80385/1388032 , 39385/5552128 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -145919/2743808 , 721527/1371904 , 105687/171488 , ... + -25/256 , 23631/2743808 , 0 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -178863/4033024 , 1178085/8066048 , ... + 1658587/2016512 , 401925/4033024 , -15/512 , ... + 43801/8066048 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -1668147/11213056 , 4193225/11213056 , ... + 375675/2803264 , 2009625/2803264 , ... + -984625/11213056 , 3/256 , 0 , 0 , 0 , 0 , 0 , 0 ; ... + -561187/2949120 , 831521/1474560 , -788801/1474560 ... + , 412643/368640 , 39385/589824 , -43801/1474560 ... + , 5/1024 , 0 , 0 , 0 , 0 , 0 ; ... + 23/1024 , 23/147456 , -43435/221184 , ... + 26795/36864 , 39385/73728 , -43801/442368 , ... + 3/256 , 0 , 0 , 0 , 0 , 0 ; ... + 79379/368640 , -1664707/2949120 , 284431/737280 , ... + 26795/294912 , 606529/737280 , 43801/589824 , ... + -15/512 , 5/1024 , 0 , 0 , 0 , 0 ; ... + 3589/27648 , -2225/6144 , 22939/73728 , ... + -26795/221184 , 39385/73728 , 43801/73728 , ... + -25/256 , 3/256 , 0 , 0 , 0 , 0 ; ... + -720623/14745600 , 10637/92160 , -89513/1474560 , ... + -5359/147456 , 39385/589824 , 3372677/3686400 , ... + 75/1024 , -15/512 , 5/1024 , 0 , 0 , 0 ; ... + -6357/81920 , 55219/276480 , -8707/61440 , ... + 5359/368640 , -39385/442368 , 43801/73728 , ... + 75/128 , -25/256 , 3/256 , 0 , 0 , 0 ; ... + -13315/884736 , 2589/65536 , -479/16384 , ... + 5359/884736 , -7877/294912 , 43801/589824 , ... + 231/256 , 75/1024 , -15/512 , 5/1024 , 0 , 0 ; ... + 8299/737280 , -7043/245760 , 5473/276480 , 0 , ... + 7877/737280 , -43801/442368 , 75/128 , ... + 75/128 , -25/256 , 3/256 , 0 , 0 ; ... + 11027/2949120 , -8461/884736 , 655/98304 , 0 , ... + 7877/1769472 , -43801/1474560 , 75/1024 , ... + 231/256 , 75/1024 , -15/512 , 5/1024 , 0 ; ... + -601/614400 , 601/245760 , -601/368640 , 0 , 0 , ... + 43801/3686400 , -25/256 , 75/128 , ... + 75/128 , -25/256 , 3/256 , 0 ; ... + -601/1474560 , 601/589824 , -601/884736 , 0 , 0 , ... + 43801/8847360 , -15/512 , 75/1024 , ... + 231/256 , 75/1024 , -15/512 , 5/1024 ] ; + + t2= [5/1024 , -15/512 , 75/1024, 231/256 , 75/1024 , -15/512 , 5/1024]; + + t3= [3/256 , -25/256 , 75/128 , 75/128 , -25/256 , 3/256]; + +I1(1:17,1:12)=t1; +I1(M_F-16:M_F,M_C-11:M_C)=fliplr(flipud(t1)); +I1(18,7:12)=t3; +for i=19:2:M_F-18 + j=(i-3)/2; + I1(i,j-1:j+5)=t2; + I1(i+1,j:j+5)=t3; +end + +% Fine to coarse +I2=zeros(M_C,M_F); + +t1=[6854313/13976576 , 2802735/3494144 , ... + 1016925/27953152 , -729595/6988288 , ... + -894315/13976576 , -1668147/6988288 , ... + -8417805/27953152 , 15525/436768 , ... + 1190685/3494144 , 89725/436768 , ... + -2161869/27953152 , -858195/6988288 , ... + -332875/13976576 , 124485/6988288 , ... + 165405/27953152 , -5409/3494144 , -9015/13976576; ... + 80385/12301312 , 1201479/3075328 , 1225647/6150656 ... + , 721527/3075328 , 1178085/24602624 , ... + 838645/6150656 , 60843/300032 , 345/6150656 , ... + -4994121/24602624 , -100125/768832 , ... + 31911/768832 , 55219/768832 , 349515/24602624 , ... + -63387/6150656 , -42305/12301312 , 5409/6150656 ... + , 9015/24602624 ; ... + -80385/5552128 , -240439/1388032 , 364155/5552128 ... + , 105687/173504 , 1658587/2776064 , 75135/694016 ... + , -2366403/5552128 , -217175/1388032 , ... + 853293/2776064 , 344085/1388032 , ... + -268539/5552128 , -78363/694016 , -64665/2776064 ... + , 5473/347008 , 29475/5552128 , -1803/1388032 , ... + -3005/5552128]; + + + + t2=[ 5/2048 , 3/512 , -15/1024 , -25/512 , ... + 75/2048 , 75/256 , 231/512 , 75/256 , ... + 75/2048 , -25/512 , -15/1024 , 3/512 , 5/2048]; + +I2(1:3,1:17)=t1; + +I2(M_C-2:M_C,M_F-16:M_F)=fliplr(flipud(t1)); + +for i=4:M_C-3 + j=2*(i-4)+1; + I2(i,j:j+12)=t2; +end +IC2F = sparse(I1); +IF2C = sparse(I2); +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/intOpMC_orders_8to8_ratio2to1.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,53 @@ +% Marks New interpolation operators +% 8th order accurate (diagonal norm) +% M=19 is the minimum amount of points on the coarse mesh + +function [IC2F,IF2C] = intOpMC_orders_8to8_ratio2to1(M_C) + +M_F=M_C*2-1; + +% Coarse to fine +I1=zeros(M_F,M_C); + +t1= [0.99340647972821530406e0 0.26374081087138783775e-1 -0.39561121630708175663e-1 0.26374081087138783775e-1 -0.65935202717846959438e-2 0 0 0 0 0 0 0 0 0 0 0; 0.31183959860287729600e0 0.94014160558849081601e0 -0.31646240838273622401e0 0.65141605588490816007e-1 -0.66040139712270400169e-3 0 0 0 0 0 0 0 0 0 0 0; -0.33163591467432411328e-1 0.11092588191512759415e1 -0.10539936193077965253e0 -0.77189144409925778387e-2 0.60418595406382404105e-1 -0.23395546718453703858e-1 0 0 0 0 0 0 0 0 0 0; -0.63951987538041216890e-1 0.56657207530367360435e0 0.56073157416571775151e0 -0.67107298938782711711e-1 0.54915118559238359570e-2 -0.17358748484912632117e-2 0 0 0 0 0 0 0 0 0 0; 0.22970968995770386894e0 -0.84424237330911973043e0 0.20693327723706358111e1 -0.42910115554542331920e0 -0.13191478386190563918e0 0.11675567167691342983e0 -0.10539821288804421121e-1 0 0 0 0 0 0 0 0 0; 0.10195433776615307267e0 -0.45344410547614678949e0 0.11049699103276728065e1 0.26297465812771521534e0 -0.38617448080010680341e-1 0.23925781250000000000e-1 -0.17631339153836246986e-2 0 0 0 0 0 0 0 0 0; -0.33882356049961665258e0 0.12718893449497745294e1 -0.16822328941798626751e1 0.17813590728082409984e1 0.11793036905675500906e0 -0.18266200597575250398e0 0.37689938246258754051e-1 -0.51502644057974593120e-2 0 0 0 0 0 0 0 0; -0.50400650482311136546e0 0.19901277318227816880e1 -0.29708252195230746436e1 0.23722122500767090037e1 0.24457622727717445252e0 -0.15152936311741758892e0 0.21886284536938453771e-1 -0.24414062500000000000e-2 0 0 0 0 0 0 0 0; 0.37882732714731866331e-2 0.12115606818363056270e0 -0.53924884156527999826e0 0.88886598075786086512e0 0.27660232216640139169e0 0.33730204543181760125e0 -0.12179633685076087365e0 0.38041730885639608773e-1 -0.47112422807823442564e-2 0 0 0 0 0 0 0; 0.41123781086486062886e0 -0.14418811327145423418e1 0.16793759240044487847e1 -0.57156511000645013503e0 0.24685906775203751929e0 0.76471858554416232639e0 -0.11045284035765094522e0 0.24149101163134162809e-1 -0.24414062500000000000e-2 0 0 0 0 0 0 0; 0.36463669929508579110e0 -0.13698743320477484358e1 0.18268133133896437122e1 -0.93074264690533336964e0 0.10888458847499176143e0 0.85685706622610101431e0 0.24359267370152174731e0 -0.13314605809973863070e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0 0 0; 0.25541802394537278164e0 -0.10497853549910180363e1 0.15921730465359157986e1 -0.96615555845660927855e0 -0.49371813550407503858e-1 0.76471858554416232639e0 0.55226420178825472608e0 -0.12074550581567081404e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0 0 0; -0.93469656691661424206e-1 0.26634337856674539612e0 -0.17694629839462736800e0 -0.64947940656920645005e-1 -0.54442294237495880713e-1 0.33730204543181760125e0 0.61880473767909428853e0 0.26629211619947726141e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0 0; -0.49445400002561969522e0 0.18168098496802059221e1 -0.23462477366129557292e1 0.11091140417405116705e1 0.98743627100815007716e-2 -0.15294371710883246528e0 0.55226420178825472608e0 0.60372752907835407022e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0 0; -0.24633538738293361369e0 0.93021448723433316305e0 -0.12527182680719820488e1 0.63698038058706249354e0 0.15554941210713108775e-1 -0.16865102271590880063e0 0.24359267370152174731e0 0.67646871560981190145e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0 0; 0.23150628239599695492e0 -0.82682792555928440531e0 0.10237569354829334077e1 -0.45129113643047460593e0 -0.10075880316409694665e-2 0.30588743421766493056e-1 -0.11045284035765094522e0 0.60372752907835407022e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0 0; 0.19518658723021413499e0 -0.70515100787561674409e0 0.88870680011413432419e0 -0.40458631782898657259e0 -0.19443676513391385969e-2 0.48186006490259657322e-1 -0.12179633685076087365e0 0.26629211619947726141e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0 0; -0.43403699494753775353e-1 0.15442805550926787092e0 -0.18997683853068679729e0 0.82584189359473172950e-1 0 -0.31213003491598462302e-2 0.22090568071530189043e-1 -0.12074550581567081404e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0 0; -0.58783367481931040778e-1 0.20994477957758583150e0 -0.25976152530269196017e0 0.11403438083569734975e0 0 -0.60232508112824571652e-2 0.34798953385931678187e-1 -0.13314605809973863070e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0 0; 0.63390647713259204145e-2 -0.22373993350504987875e-1 0.27185871992161665013e-1 -0.11561529976981026786e-1 0 0 -0.22541395991357335758e-2 0.24149101163134162809e-1 -0.11962890625000000000e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0 0; 0.10649583863025939259e-1 -0.37634916628131671889e-1 0.45812371547331598336e-1 -0.19540204529147604911e-1 0 0 -0.43498691732414597734e-2 0.38041730885639608773e-1 -0.13191478386190563918e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2 0; -0.45575492476359756866e-3 0.15951422366725914903e-2 -0.19141706840071097884e-2 0.79757111833629574515e-3 0 0 0 -0.24641939962381798784e-2 0.23925781250000000000e-1 -0.11962890625000000000e0 0.59814453125000000000e0 0.59814453125000000000e0 -0.11962890625000000000e0 0.23925781250000000000e-1 -0.24414062500000000000e-2 0; -0.87948159845213680356e-3 0.30781855945824788125e-2 -0.36938227134989745750e-2 0.15390927972912394062e-2 0 0 0 -0.47552163607049510966e-2 0.37689938246258754051e-1 -0.13191478386190563918e0 0.26382956772381127836e0 0.67021304034523590205e0 0.26382956772381127836e0 -0.13191478386190563918e0 0.37689938246258754051e-1 -0.47112422807823442564e-2;]; + + + t2= [-0.456778318652801649905139026187255979136056340856723240855601037075101451671715366e81 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 -0.1827113274611206599620556104749023916544225363426892963422404148300405806686861464e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.1856585148354920384923865861096125662293072684152233190798249652677391616531107981e82 / 0.2770141785660523684734143913470637620565185365865679672509451726827594519874538713e82 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 -0.1827113274611206599620556104749023916544225363426892963422404148300405806686861464e82 / 0.13850708928302618423670719567353188102825926829328398362547258634137972599372693565e83 0.3654226549222413199241112209498047833088450726853785926844808296600811613373722928e82 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83 -0.456778318652801649905139026187255979136056340856723240855601037075101451671715366e81 / 0.96954962498118328965695036971472316719781487805298788537830810438965808195608854955e83;]; + + t3= [-0.5e1 / 0.2048e4 0.49e2 / 0.2048e4 -0.245e3 / 0.2048e4 0.1225e4 / 0.2048e4 0.1225e4 / 0.2048e4 -0.245e3 / 0.2048e4 0.49e2 / 0.2048e4 -0.5e1 / 0.2048e4;]; +I1(1:23,1:16)=t1; +I1(M_F-22:M_F,M_C-15:M_C)=fliplr(flipud(t1)); +I1(24,9:16)=t3; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DEBUGGING + +for i=25:2:M_F-24 + j=(i-3)/2; + I1(i,j-2:j+6)=t2; + I1(i+1,j-1:j+6)=t3; +end + +% Fine to coarse: A +I2=zeros(M_C,M_F); + +t1=[0.49670323986410765203e0 0.80670591743192512511e0 -0.14476656476698073533e-1 -0.19497555250083395132e0 0.16074268813765884450e0 0.22100911221277072755e0 -0.53042418939472250452e0 -0.86254139670456354517e0 0.64231825172866668299e-2 0.69727164011248914487e0 0.61825742343094006838e0 0.43307239695721032895e0 -0.15848187861199173328e0 -0.83836831742920925562e0 -0.41767239062238727391e0 0.39252899650233765024e0 0.33094736964907844809e0 -0.73592865087013788440e-1 -0.99669762780958210515e-1 0.10748160731101219580e-1 0.18056833808814782786e-1 -0.77275234787125894193e-3 -0.14911993994710636483e-2; 0.25487859584304862664e-2 0.47007080279424540800e0 0.93589176759289320118e-1 0.33386224041716468423e0 -0.11418395501435496457e0 -0.18998278818813064037e0 0.38484431284491986603e0 0.65828018436035862232e0 0.39704539050575728856e-1 -0.47252462545568042557e0 -0.44892698918501097924e0 -0.34402935194949605213e0 0.87284452472801643396e-1 0.59539401290875351190e0 0.30484430526276345964e0 -0.27096308216867871783e0 -0.23108785344796321535e0 0.50608234918774219796e-1 0.68801842319351676214e-1 -0.73322707316320135247e-2 -0.12333488857215226757e-1 0.52275043402033040521e-3 0.10087644967132781708e-2; -0.22656973277393401539e-1 -0.93771184228725468159e0 -0.52699680965389826267e-1 0.19581430555012001670e1 0.16586148101136731602e1 0.27435837109255836264e1 -0.30164708462284474624e1 -0.58235016129647971089e1 -0.10472767830023645070e1 0.32615209891555982371e1 0.35478595826728208891e1 0.30921640208240511054e1 -0.34364793368678654581e0 -0.45566547247355317663e1 -0.24329078834671892590e1 0.19882413967858906122e1 0.17259601262271516763e1 -0.36895458453625989435e0 -0.50448363278283993228e0 0.52797763052066775846e-1 0.88972343374038343284e-1 -0.37175165926095403240e-2 -0.71737841052106668688e-2; 0.21626748627943672418e-2 0.27636709246281031633e-1 -0.55259484658035070394e-3 -0.33553649469391355855e-1 -0.49244245327794891233e-1 0.93489376222103426899e-1 0.45734620580442859300e0 0.66579609152388478842e0 0.24716623315221892947e0 -0.15893464065423852815e0 -0.25881084331023040874e0 -0.26865808253702445366e0 -0.18060020510041282529e-1 0.30841043055726240020e0 0.17712461121229459766e0 -0.12549015561536110372e0 -0.11250298507032009025e0 0.22964117700293735857e-1 0.31709446610808019222e-1 -0.32149051440245356510e-2 -0.54335286230388551025e-2 0.22177994574852164638e-3 0.42797426992744435493e-3;]; + + + + t2=[-0.23556211403911721282e-2 -0.12207031250000000000e-2 0.18844969123129377026e-1 0.11962890625000000000e-1 -0.65957391930952819589e-1 -0.59814453125000000000e-1 0.13191478386190563918e0 0.29907226562500000000e0 0.33510652017261795103e0 0.29907226562500000000e0 0.13191478386190563918e0 -0.59814453125000000000e-1 -0.65957391930952819589e-1 0.11962890625000000000e-1 0.18844969123129377026e-1 -0.12207031250000000000e-2 -0.23556211403911721282e-2;]; + + +I2(1:4,1:23)=t1; + +I2(M_C-3:M_C,M_F-22:M_F)=fliplr(flipud(t1)); + +for i=5:M_C-4 + j=2*(i-5)+1; + I2(i,j:j+16)=t2; +end + +IC2F = sparse(I1); +IF2C = sparse(I2); + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D1Staggered.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,86 @@ +classdef D1Staggered < sbp.OpSet + properties + % x_primal: "primal" grid with m points. Equidistant. Called Plus grid in Ossian's paper. + % x_dual: "dual" grid with m+1 points. Called Minus grid in Ossian's paper. + + % D1_primal takes FROM dual grid TO primal grid + % D1_dual takes FROM primal grid TO dual grid + + D1_primal % SBP operator approximating first derivative + D1_dual % SBP operator approximating first derivative + H_primal % Norm matrix + H_dual % Norm matrix + H_primalI % H^-1 + H_dualI % H^-1 + e_primal_l % Left boundary operator + e_dual_l % Left boundary operator + e_primal_r % Right boundary operator + e_dual_r % Right boundary operator + m % Number of grid points. + m_primal % Number of grid points. + m_dual % Number of grid points. + h % Step size + x_primal % grid + x_dual % grid + x + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D1Staggered(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + + m_primal = m; + m_dual = m+1; + + switch order + case 2 + [x_primal, x_dual, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_2(m, L); + case 4 + [x_primal, x_dual, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_4(m, L); + case 6 + [x_primal, x_dual, Pp, Pm, Qp, Qm] = sbp.implementations.d1_staggered_6(m, L); + otherwise + error('Invalid operator order %d.',order); + end + + obj.m = m; + obj.m_primal = m_primal; + obj.m_dual = m_dual; + obj.x_primal = x_l + x_primal'; + obj.x_dual = x_l + x_dual'; + + D1_primal = Pp\Qp; + D1_dual = Pm\Qm; + + obj.D1_primal = D1_primal; + obj.D1_dual = D1_dual; + obj.H_primal = Pp; + obj.H_dual = Pm; + + obj.e_primal_l = sparse(m_primal,1); + obj.e_primal_r = sparse(m_primal,1); + obj.e_primal_l(1) = 1; + obj.e_primal_r(m_primal) = 1; + + obj.e_dual_l = sparse(m_dual,1); + obj.e_dual_r = sparse(m_dual,1); + obj.e_dual_l(1) = 1; + obj.e_dual_r(m_dual) = 1; + + obj.H_primalI = inv(obj.H_primal); + obj.H_dualI = inv(obj.H_dual); + + obj.borrowing = []; + obj.x = []; + + end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D1StaggeredUpwind.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,95 @@ +classdef D1StaggeredUpwind < sbp.OpSet + % Compatible staggered and upwind operators by Ken Mattsson and Ossian O'reilly + properties + % x_primal: "primal" grid with m points. Equidistant. Called Plus grid in Ossian's paper. + % x_dual: "dual" grid with m+1 points. Called Minus grid in Ossian's paper. + + % D1_primal takes FROM dual grid TO primal grid + % D1_dual takes FROM primal grid TO dual grid + + D1_primal % SBP operator approximating first derivative + D1_dual % SBP operator approximating first derivative + + Dplus_primal % Upwind operator on primal grid + Dminus_primal % Upwind operator on primal grid + Dplus_dual % Upwind operator on dual grid + Dminus_dual % Upwind operator on dual grid + + H_primal % Norm matrix + H_dual % Norm matrix + H_primalI % H^-1 + H_dualI % H^-1 + e_primal_l % Left boundary operator + e_dual_l % Left boundary operator + e_primal_r % Right boundary operator + e_dual_r % Right boundary operator + m % Number of grid points. + m_primal % Number of grid points. + m_dual % Number of grid points. + h % Step size + x_primal % grid + x_dual % grid + x + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D1StaggeredUpwind(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + + m_primal = m; + m_dual = m+1; + + switch order + case 2 + [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... + obj.H_primalI, obj.H_dualI,... + obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... + obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_2(m, L); + case 4 + [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... + obj.H_primalI, obj.H_dualI,... + obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... + obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_4(m, L); + case 6 + [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... + obj.H_primalI, obj.H_dualI,... + obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... + obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_6(m, L); + case 8 + [obj.x_primal, obj.x_dual, obj.H_primal, obj.H_dual,... + obj.H_primalI, obj.H_dualI,... + obj.D1_primal, obj.D1_dual, obj.Dplus_primal, obj.Dminus_primal,... + obj.Dplus_dual, obj.Dminus_dual] = sbp.implementations.d1_staggered_upwind_8(m, L); + otherwise + error('Invalid operator order %d.',order); + end + + obj.m = m; + obj.m_primal = m_primal; + obj.m_dual = m_dual; + obj.h = L/(m-1); + + obj.e_primal_l = sparse(m_primal,1); + obj.e_primal_r = sparse(m_primal,1); + obj.e_primal_l(1) = 1; + obj.e_primal_r(m_primal) = 1; + + obj.e_dual_l = sparse(m_dual,1); + obj.e_dual_r = sparse(m_dual,1); + obj.e_dual_l(1) = 1; + obj.e_dual_r(m_dual) = 1; + + obj.borrowing = []; + obj.x = []; + + end + + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end +end
--- a/+sbp/D2Variable.m Mon Jun 04 17:48:19 2018 +0200 +++ b/+sbp/D2Variable.m Sun Jun 17 15:29:46 2018 -0700 @@ -26,22 +26,39 @@ obj.x = linspace(x_l,x_r,m)'; switch order + + case 6 + + [obj.H, obj.HI, obj.D1, obj.D2, ... + ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... + obj.d1_l, obj.d1_r] = ... + sbp.implementations.d4_variable_6(m, obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + case 4 [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_4(m,obj.h); obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R case 2 [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_2(m,obj.h); obj.borrowing.M.d1 = 0.3636363636; % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R otherwise error('Invalid operator order %d.',order); end - + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, obj.m = m; obj.M = []; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2VariablePeriodic.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,71 @@ +classdef D2VariablePeriodic < sbp.OpSet + 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 + D2 % SBP operator for second derivative + M % Norm matrix, second derivative + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D2VariablePeriodic(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + obj.h = L/m; + x = linspace(x_l,x_r,m+1)'; + obj.x = x(1:end-1); + + switch order + + case 6 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_6(m,obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + + case 4 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_4(m,obj.h); + obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R + case 2 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_2(m,obj.h); + obj.borrowing.M.d1 = 0.3636363636; + % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R + + otherwise + error('Invalid operator order %d.',order); + end + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, + + obj.m = m; + obj.M = []; + end + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end + + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpAWW.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,67 @@ +classdef InterpAWW < sbp.InterpOps + properties + + % Interpolation operators + IC2F + IF2C + + % Orders used on coarse and fine sides + order_C + order_F + + % Grid points, refinement ratio. + ratio + m_C + m_F + + % Boundary accuracy of IC2F and IF2C. + acc_C2F + acc_F2C + end + + methods + % accOp : String, 'C2F' or 'F2C'. Specifies which of the operators + % should have higher accuracy. + function obj = InterpAWW(m_C,m_F,order_C,order_F,accOp) + + ratio = (m_F-1)/(m_C-1); + h_C = 1; + + assert(order_C == order_F,... + 'Error: Different orders of accuracy not available'); + + switch ratio + case 2 + switch order_C + case 2 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_2to2_ratio2to1(m_C, h_C, accOp); + case 4 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_4to4_ratio2to1(m_C, h_C, accOp); + case 6 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_6to6_ratio2to1(m_C, h_C, accOp); + case 8 + [IC2F,IF2C] = sbp.implementations.intOpAWW_orders_8to8_ratio2to1(m_C, h_C, accOp); + otherwise + error(['Order ' num2str(order_C) ' not available.']); + end + otherwise + error(['Grid ratio ' num2str(ratio) ' not available']); + end + + obj.IC2F = IC2F; + obj.IF2C = IF2C; + obj.order_C = order_C; + obj.order_F = order_F; + obj.ratio = ratio; + obj.m_C = m_C; + obj.m_F = m_F; + + end + + function str = string(obj) + str = [class(obj) '_orders' num2str(obj.order_F) 'to'... + num2str(obj.order_C) '_ratio' num2str(obj.ratio) 'to1']; + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpMC.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,61 @@ +classdef InterpMC < sbp.InterpOps + properties + + % Interpolation operators + IC2F + IF2C + + % Orders used on coarse and fine sides + order_C + order_F + + % Grid points, refinement ratio. + ratio + m_C + m_F + end + + methods + function obj = InterpMC(m_C,m_F,order_C,order_F) + + ratio = (m_F-1)/(m_C-1); + + assert(order_C == order_F,... + 'Error: Different orders of accuracy not available'); + + switch ratio + case 2 + switch order_C + case 2 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_2to2_ratio2to1(m_C); + case 4 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_4to4_ratio2to1(m_C); + case 6 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_6to6_ratio2to1(m_C); + case 8 + [IC2F,IF2C] = sbp.implementations.intOpMC_orders_8to8_ratio2to1(m_C); + otherwise + error(['Order ' num2str(order_C) ' not available.']); + end + otherwise + error(['Grid ratio ' num2str(ratio) ' not available']); + end + + obj.IC2F = IC2F; + obj.IF2C = IF2C; + obj.order_C = order_C; + obj.order_F = order_F; + obj.ratio = ratio; + obj.m_C = m_C; + obj.m_F = m_F; + + + end + + function str = string(obj) + str = [class(obj) '_orders' num2str(obj.order_F) 'to'... + num2str(obj.order_C) '_ratio' num2str(obj.ratio) 'to1']; + end + + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/InterpOps.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,14 @@ +classdef (Abstract) InterpOps + properties (Abstract) + % C and F may refer to coarse and fine, but it's not a must. + IC2F % Interpolation operator from "C" to "F" + IF2C % Interpolation operator from "F" to "C" + + end + + methods (Abstract) + % Returns a string representation of the type of operator. + str = string(obj) + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Elastic2dVariable.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,420 @@ +classdef Elastic2dVariable < scheme.Scheme + +% Discretizes the elastic wave equation: +% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrices for varible coefficients + LAMBDA % Variable coefficient, related to dilation + MU % Shear modulus, variable coefficient + RHO, RHOi % Density, variable + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_lambda + D2_mu + + % Traction operators used for BC + T_l, T_r + tau_l, tau_r + + H, Hi % Inner products + phi % Borrowing constant for (d1 - e^T*D1) from R + gamma % Borrowing constant for d1 from M + H11 % First element of H + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + E % E{i}^T picks out component i + + H_boundary % Boundary inner products + + % Kroneckered norms and coefficients + RHOi_kron + Hi_kron + end + + methods + + function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('lambda_fun', @(x,y) 0*x+1); + default_arg('mu_fun', @(x,y) 0*x+1); + default_arg('rho_fun', @(x,y) 0*x+1); + dim = 2; + + assert(isa(g, 'grid.Cartesian')) + + lambda = grid.evalOn(g, lambda_fun); + mu = grid.evalOn(g, mu_fun); + rho = grid.evalOn(g, rho_fun); + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + lim = g.lim; + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, order); + end + + % Borrowing constants + for i = 1:dim + beta = ops{i}.borrowing.R.delta_D; + obj.H11{i} = ops{i}.borrowing.H11; + obj.phi{i} = beta/obj.H11{i}; + obj.gamma{i} = ops{i}.borrowing.M.d1; + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + %====== Assemble full operators ======== + LAMBDA = spdiag(lambda); + obj.LAMBDA = LAMBDA; + MU = spdiag(mu); + obj.MU = MU; + RHO = spdiag(rho); + obj.RHO = RHO; + obj.RHOi = inv(RHO); + + obj.D1 = cell(dim,1); + obj.D2_lambda = cell(dim,1); + obj.D2_mu = cell(dim,1); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 + for i = 1:dim + obj.D2_lambda{i} = sparse(m_tot); + obj.D2_mu{i} = sparse(m_tot); + end + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m(2) + D_lambda = D2{1}(lambda(ind(:,i))); + D_mu = D2{1}(mu(ind(:,i))); + + p = ind(:,i); + obj.D2_lambda{1}(p,p) = D_lambda; + obj.D2_mu{1}(p,p) = D_mu; + end + + for i = 1:m(1) + D_lambda = D2{2}(lambda(ind(i,:))); + D_mu = D2{2}(mu(ind(i,:))); + + p = ind(i,:); + obj.D2_lambda{2}(p,p) = D_lambda; + obj.D2_mu{2}(p,p) = D_mu; + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary = cell(dim,1); + obj.H_boundary{1} = H{2}; + obj.H_boundary{2} = H{1}; + + % E{i}^T picks out component i. + E = cell(dim,1); + I = speye(m_tot,m_tot); + for i = 1:dim + e = sparse(dim,1); + e(i) = 1; + E{i} = kron(I,e); + end + obj.E = E; + + % Differentiation matrix D (without SAT) + D2_lambda = obj.D2_lambda; + D2_mu = obj.D2_mu; + D1 = obj.D1; + D = sparse(dim*m_tot,dim*m_tot); + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + for i = 1:dim + for j = 1:dim + D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +... + db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ... + ); + D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +... + db(i,j)*D1{j}*MU*D1{i}*E{j}' + ... + D2_mu{j}*E{i}' ... + ); + end + end + obj.D = D; + %=========================================% + + % Numerical traction operators for BC. + % Because d1 =/= e0^T*D1, the numerical tractions are different + % at every boundary. + T_l = cell(dim,1); + T_r = cell(dim,1); + tau_l = cell(dim,1); + tau_r = cell(dim,1); + % tau^{j}_i = sum_k T^{j}_{ik} u_k + + d1_l = obj.d1_l; + d1_r = obj.d1_r; + e_l = obj.e_l; + e_r = obj.e_r; + D1 = obj.D1; + + % Loop over boundaries + for j = 1:dim + T_l{j} = cell(dim,dim); + T_r{j} = cell(dim,dim); + tau_l{j} = cell(dim,1); + tau_r{j} = cell(dim,1); + + % Loop over components + for i = 1:dim + tau_l{j}{i} = sparse(m_tot,dim*m_tot); + tau_r{j}{i} = sparse(m_tot,dim*m_tot); + for k = 1:dim + T_l{j}{i,k} = ... + -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... + -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... + -d(i,k)*MU*e_l{j}*d1_l{j}'; + + T_r{j}{i,k} = ... + d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... + +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... + +d(i,k)*MU*e_r{j}*d1_r{j}'; + + tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; + tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; + end + + end + end + obj.T_l = T_l; + obj.T_r = T_r; + obj.tau_l = tau_l; + obj.tau_r = tau_r; + + % Kroneckered norms and coefficients + I_dim = speye(dim); + obj.RHOi_kron = kron(obj.RHOi, I_dim); + obj.Hi_kron = kron(obj.Hi, I_dim); + + % Misc. + obj.m = m; + obj.h = h; + obj.order = order; + obj.grid = g; + obj.dim = dim; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a cell array of strings specifying the type of boundary condition for each component. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) + default_arg('type',{'free','free'}); + default_arg('parameter', []); + + % j is the coordinate direction of the boundary + % nj: outward unit normal component. + % nj = -1 for west, south, bottom boundaries + % nj = 1 for east, north, top boundaries + [j, nj] = obj.get_boundary_number(boundary); + switch nj + case 1 + e = obj.e_r; + d = obj.d1_r; + tau = obj.tau_r{j}; + T = obj.T_r{j}; + case -1 + e = obj.e_l; + d = obj.d1_l; + tau = obj.tau_l{j}; + T = obj.T_l{j}; + end + + E = obj.E; + Hi = obj.Hi; + H_gamma = obj.H_boundary{j}; + LAMBDA = obj.LAMBDA; + MU = obj.MU; + RHOi = obj.RHOi; + + dim = obj.dim; + m_tot = obj.grid.N(); + + RHOi_kron = obj.RHOi_kron; + Hi_kron = obj.Hi_kron; + + % Preallocate + closure = sparse(dim*m_tot, dim*m_tot); + penalty = cell(dim,1); + for k = 1:dim + penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j)); + end + + % Loop over components that we (potentially) have different BC on + for k = 1:dim + switch type{k} + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + + tuning = 1.2; + phi = obj.phi{j}; + h = obj.h(j); + h11 = obj.H11{j}*h; + gamma = obj.gamma{j}; + + a_lambda = dim/h11 + 1/(h11*phi); + a_mu_i = 2/(gamma*h); + a_mu_ij = 2/h11 + 1/(h11*phi); + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... + + d(i,j)* a_mu_i*MU ... + + db(i,j)*a_mu_ij*MU ); + + % Loop over components that Dirichlet penalties end up on + for i = 1:dim + C = T{k,i}; + A = -d(i,k)*alpha(i,j); + B = A + C; + closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); + penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; + end + + % Free boundary condition + case {'F','f','Free','free','traction','Traction','t','T'} + closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); + penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + tuning = 1.2; + % tuning = 20.2; + error('Interface not implemented'); + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [j, nj] = get_boundary_number(obj, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + nj = 1; + end + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operatr = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Heat2dVariable.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,262 @@ +classdef Heat2dVariable < scheme.Scheme + +% Discretizes the Laplacian with variable coefficent, +% In the Heat equation way (i.e., the discretization matrix is not necessarily +% symmetric) +% u_t = div * (kappa * grad u ) +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrix for variable coefficients + KAPPA % Variable coefficient + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_kappa + + H, Hi % Inner products + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + + H_boundary % Boundary inner products + + end + + methods + + function obj = Heat2dVariable(g ,order, kappa_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('kappa_fun', @(x,y) 0*x+1); + dim = 2; + + assert(isa(g, 'grid.Cartesian')) + + kappa = grid.evalOn(g, kappa_fun); + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + lim = g.lim; + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, order); + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + %====== Assemble full operators ======== + KAPPA = spdiag(kappa); + obj.KAPPA = KAPPA; + + obj.D1 = cell(dim,1); + obj.D2_kappa = cell(dim,1); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 + for i = 1:dim + obj.D2_kappa{i} = sparse(m_tot); + end + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m(2) + D_kappa = D2{1}(kappa(ind(:,i))); + p = ind(:,i); + obj.D2_kappa{1}(p,p) = D_kappa; + end + + for i = 1:m(1) + D_kappa = D2{2}(kappa(ind(i,:))); + p = ind(i,:); + obj.D2_kappa{2}(p,p) = D_kappa; + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary = cell(dim,1); + obj.H_boundary{1} = H{2}; + obj.H_boundary{2} = H{1}; + + % Differentiation matrix D (without SAT) + D2_kappa = obj.D2_kappa; + D1 = obj.D1; + D = sparse(m_tot,m_tot); + for i = 1:dim + D = D + D2_kappa{i}; + end + obj.D = D; + %=========================================% + + % Misc. + obj.m = m; + obj.h = h; + obj.order = order; + obj.grid = g; + obj.dim = dim; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) + default_arg('type','Neumann'); + default_arg('parameter', []); + + % j is the coordinate direction of the boundary + % nj: outward unit normal component. + % nj = -1 for west, south, bottom boundaries + % nj = 1 for east, north, top boundaries + [j, nj] = obj.get_boundary_number(boundary); + switch nj + case 1 + e = obj.e_r; + d = obj.d1_r; + case -1 + e = obj.e_l; + d = obj.d1_l; + end + + Hi = obj.Hi; + H_gamma = obj.H_boundary{j}; + KAPPA = obj.KAPPA; + kappa_gamma = e{j}'*KAPPA*e{j}; + + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); + penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; + + % Free boundary condition + case {'N','n','neumann','Neumann'} + closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); + penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + error('Interface not implemented'); + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [j, nj] = get_boundary_number(obj, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + nj = 1; + end + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operatr = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
--- a/+scheme/LaplaceCurvilinear.m Mon Jun 04 17:48:19 2018 +0200 +++ b/+scheme/LaplaceCurvilinear.m Sun Jun 17 15:29:46 2018 -0700 @@ -38,22 +38,29 @@ du_n, dv_n gamm_u, gamm_v lambda + + interpolation_type end methods % Implements a*div(b*grad(u)) as a SBP scheme % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) - function obj = LaplaceCurvilinear(g ,order, a, b, opSet) + function obj = LaplaceCurvilinear(g ,order, a, b, opSet, interpolation_type) default_arg('opSet',@sbp.D2Variable); default_arg('a', 1); default_arg('b', 1); + default_arg('interpolation_type','AWW'); if b ~=1 error('Not implemented yet') end - assert(isa(g, 'grid.Curvilinear')) + % assert(isa(g, 'grid.Curvilinear')) + if isa(a, 'function_handle') + a = grid.evalOn(g, a); + a = spdiag(a); + end m = g.size(); m_u = m(1); @@ -209,6 +216,7 @@ obj.h = [h_u h_v]; obj.order = order; obj.grid = g; + obj.interpolation_type = interpolation_type; obj.a = a; obj.b = b; @@ -269,13 +277,70 @@ end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); + [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + Hi = obj.Hi; + a = obj.a; + + m_u = length(H_b_u); + m_v = length(H_b_v); + + grid_ratio = m_u/m_v; + if grid_ratio ~= 1 + + [ms, index] = sort([m_u, m_v]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_v2u_good = interpOpSet.IF2C; + I_v2u_bad = interpOpSet.IF2C; + I_u2v_good = interpOpSet.IC2F; + I_u2v_bad = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_v2u_good = interpOpSet.IC2F; + I_v2u_bad = interpOpSet.IC2F; + I_u2v_good = interpOpSet.IF2C; + I_u2v_bad = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_v2u_good = interpOpSetF2C.IF2C; + I_v2u_bad = interpOpSetC2F.IF2C; + I_u2v_good = interpOpSetC2F.IC2F; + I_u2v_bad = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_v2u_good = interpOpSetC2F.IC2F; + I_v2u_bad = interpOpSetF2C.IC2F; + I_u2v_good = interpOpSetF2C.IF2C; + I_u2v_bad = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_v2u_good = speye(m_u,m_u); + I_v2u_bad = speye(m_u,m_u); + I_u2v_good = speye(m_u,m_u); + I_u2v_bad = speye(m_u,m_u); + end + % u denotes the solution in the own domain % v denotes the solution in the neighbour domain tuning = 1.2; % tuning = 20.2; - [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); - [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - u = obj; v = neighbour_scheme; @@ -284,18 +349,24 @@ b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; - tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); - tau1 = tuning * spdiag(tau1); - tau2 = 1/2; + tau_u = -1./(4*b1_u) -1./(4*b2_u); + tau_v = -1./(4*b1_v) -1./(4*b2_v); + + tau_u = tuning * spdiag(tau_u); + tau_v = tuning * spdiag(tau_v); + beta_u = tau_v; - sig1 = -1/2; - sig2 = 0; + closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... + a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*I_u2v_good*e_u' + ... + a*1/2*Hi*d_u*H_b_u*e_u' + ... + -a*1/2*Hi*e_u*H_b_u*d_u'; - tau = (e_u*tau1 + tau2*d_u)*H_b_u; - sig = (sig1*e_u + sig2*d_u)*H_b_u; + penalty = -a*Hi*e_u*tau_u*H_b_u*I_v2u_good*e_v' + ... + -a*Hi*e_u*H_b_u*I_v2u_bad*beta_u*e_v' + ... + -a*1/2*Hi*d_u*H_b_u*I_v2u_good*e_v' + ... + -a*1/2*Hi*e_u*H_b_u*I_v2u_bad*d_v'; + - closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u'); - penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v'); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Schrodinger2d.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,374 @@ +classdef Schrodinger2d < scheme.Scheme + +% Discretizes the Laplacian with constant coefficent, +% in the Schrödinger equation way (i.e., the discretization matrix is not necessarily +% definite) +% u_t = a*i*Laplace u +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrix for variable coefficients + a % Constant coefficient + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2 + + H, Hi % Inner products + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + e_w, e_e, e_s, e_n + d_w, d_e, d_s, d_n + + H_boundary % Boundary inner products + + interpolation_type % MC or AWW + + end + + methods + + function obj = Schrodinger2d(g ,order, a, opSet, interpolation_type) + default_arg('interpolation_type','AWW'); + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('a',1); + dim = 2; + + assert(isa(g, 'grid.Cartesian')) + if isa(a, 'function_handle') + a = grid.evalOn(g, a); + a = spdiag(a); + end + + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + lim = {xlim, ylim}; + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, order); + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + % Constant coeff D2 + for i = 1:dim + D2{i} = D2{i}(ones(m(i),1)); + end + + %====== Assemble full operators ======== + obj.D1 = cell(dim,1); + obj.D2 = cell(dim,1); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 + obj.D2{1} = kron(D2{1},I{2}); + obj.D2{2} = kron(I{1},D2{2}); + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary = cell(dim,1); + obj.H_boundary{1} = H{2}; + obj.H_boundary{2} = H{1}; + + % Differentiation matrix D (without SAT) + D2 = obj.D2; + D = sparse(m_tot,m_tot); + for j = 1:dim + D = D + a*1i*D2{j}; + end + obj.D = D; + %=========================================% + + % Misc. + obj.m = m; + obj.h = h; + obj.order = order; + obj.grid = g; + obj.dim = dim; + obj.a = a; + obj.e_w = obj.e_l{1}; + obj.e_e = obj.e_r{1}; + obj.e_s = obj.e_l{2}; + obj.e_n = obj.e_r{2}; + obj.d_w = obj.d1_l{1}; + obj.d_e = obj.d1_r{1}; + obj.d_s = obj.d1_l{2}; + obj.d_n = obj.d1_r{2}; + obj.interpolation_type = interpolation_type; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) + default_arg('type','Neumann'); + default_arg('parameter', []); + + % j is the coordinate direction of the boundary + % nj: outward unit normal component. + % nj = -1 for west, south, bottom boundaries + % nj = 1 for east, north, top boundaries + [j, nj] = obj.get_boundary_number(boundary); + switch nj + case 1 + e = obj.e_r; + d = obj.d1_r; + case -1 + e = obj.e_l; + d = obj.d1_l; + end + + Hi = obj.Hi; + H_gamma = obj.H_boundary{j}; + a = e{j}'*obj.a*e{j}; + + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); + penalty = -nj*Hi*d{j}*a*1i*H_gamma; + + % Free boundary condition + case {'N','n','neumann','Neumann'} + closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); + penalty = nj*Hi*e{j}*a*1i*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + % Get neighbour boundary operator + + [coord_nei, n_nei] = get_boundary_number(obj, neighbour_boundary); + [coord, n] = get_boundary_number(obj, boundary); + switch n_nei + case 1 + % North or east boundary + e_neighbour = neighbour_scheme.e_r; + d_neighbour = neighbour_scheme.d1_r; + case -1 + % South or west boundary + e_neighbour = neighbour_scheme.e_l; + d_neighbour = neighbour_scheme.d1_l; + end + + e_neighbour = e_neighbour{coord_nei}; + d_neighbour = d_neighbour{coord_nei}; + H_gamma = obj.H_boundary{coord}; + Hi = obj.Hi; + a = obj.a; + + switch coord_nei + case 1 + m_neighbour = neighbour_scheme.m(2); + case 2 + m_neighbour = neighbour_scheme.m(1); + end + + switch coord + case 1 + m = obj.m(2); + case 2 + m = obj.m(1); + end + + switch n + case 1 + % North or east boundary + e = obj.e_r; + d = obj.d1_r; + case -1 + % South or west boundary + e = obj.e_l; + d = obj.d1_l; + end + e = e{coord}; + d = d{coord}; + + Hi = obj.Hi; + sigma = -n*1i*a/2; + tau = -n*(1i*a)'/2; + + grid_ratio = m/m_neighbour; + if grid_ratio ~= 1 + + [ms, index] = sort([m, m_neighbour]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_neighbour2local_e = interpOpSet.IF2C; + I_neighbour2local_d = interpOpSet.IF2C; + I_local2neighbour_e = interpOpSet.IC2F; + I_local2neighbour_d = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_neighbour2local_e = interpOpSet.IC2F; + I_neighbour2local_d = interpOpSet.IC2F; + I_local2neighbour_e = interpOpSet.IF2C; + I_local2neighbour_d = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_neighbour2local_e = interpOpSetF2C.IF2C; + I_neighbour2local_d = interpOpSetC2F.IF2C; + I_local2neighbour_e = interpOpSetC2F.IC2F; + I_local2neighbour_d = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_neighbour2local_e = interpOpSetC2F.IC2F; + I_neighbour2local_d = interpOpSetF2C.IC2F; + I_local2neighbour_e = interpOpSetF2C.IF2C; + I_local2neighbour_d = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_neighbour2local_e = speye(m,m); + I_neighbour2local_d = speye(m,m); + I_local2neighbour_e = speye(m,m); + I_local2neighbour_d = speye(m,m); + end + + closure = tau*Hi*d*H_gamma*e' + sigma*Hi*e*H_gamma*d'; + penalty = -tau*Hi*d*H_gamma*I_neighbour2local_e*e_neighbour' ... + -sigma*Hi*e*H_gamma*I_neighbour2local_d*d_neighbour'; + + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [j, nj] = get_boundary_number(obj, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + nj = 1; + end + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operator = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Staggered1DAcoustics.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,231 @@ +classdef Staggered1DAcoustics < scheme.Scheme + properties + m % Number of points of primal grid in each direction, possibly a vector + h % Grid spacing + + % Grids + grid % Total grid object + grid_primal + grid_dual + + order % Order accuracy for the approximation + + H % Combined norm + Hi % Inverse + D % Semi-discrete approximation matrix + + D1_primal + D1_dual + + % Pick out left or right boundary + e_l + e_r + + % Initial data + v0 + + % Pick out primal or dual component + e_primal + e_dual + + % System matrices + A + B + end + + + methods + % Scheme for A*u_t + B u_x = 0, + % u = [p, v]; + % A: Diagonal and A > 0, + % A = [a1, 0; + % 0, a2] + % + % B = B^T and with diagonal entries = 0. + % B = [0 b + % b 0] + % Here we store p on the primal grid and v on the dual + function obj = Staggered1DAcoustics(g, order, A, B) + default_arg('A',[1, 0; 0, 1]); + default_arg('B',[0, 1; 1, 0]); + + obj.order = order; + obj.A = A; + obj.B = B; + + % Grids + obj.m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + + obj.grid = g; + obj.grid_primal = g.grids{1}; + obj.grid_dual = g.grids{2}; + + % Get operators + ops = sbp.D1StaggeredUpwind(obj.m, xlim, order); + obj.h = ops.h; + + % Build combined operators + H_primal = ops.H_primal; + H_dual = ops.H_dual; + obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } ); + obj.Hi = inv(obj.H); + + D1_primal = ops.D1_primal; + D1_dual = ops.D1_dual; + D = {[], -1/A(1,1)*B(1,2)*D1_primal;... + -1/A(2,2)*B(2,1)*D1_dual, []}; + obj.D = blockmatrix.toMatrix(D); + obj.D1_primal = D1_primal; + obj.D1_dual = D1_dual; + + % Combined boundary operators + e_primal_l = ops.e_primal_l; + e_primal_r = ops.e_primal_r; + e_dual_l = ops.e_dual_l; + e_dual_r = ops.e_dual_r; + e_l = {e_primal_l, [];... + [] , e_dual_l}; + e_r = {e_primal_r, [];... + [] , e_dual_r}; + obj.e_l = blockmatrix.toMatrix(e_l); + obj.e_r = blockmatrix.toMatrix(e_r); + + % Pick out first or second component of solution + N_primal = obj.grid_primal.N(); + N_dual = obj.grid_dual.N(); + obj.e_primal = [speye(N_primal, N_primal); sparse(N_dual, N_primal)]; + obj.e_dual = [sparse(N_primal, N_dual); speye(N_dual, N_dual)]; + + + end + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary, type) + + default_arg('type','p'); + + % type = 'p' => boundary condition for p + % type = 'v' => boundary condition for v + % type = 'characteristic' => incoming characteristc = data + % No other types implemented yet + + % BC on the form Lu - g = 0; + + % Need a transformation T such that w = T^{−1}*u is a + % meaningful change of variables and T^T*B*T is block-diagonal, + % For linear acoustics, T = T_C meets these criteria. + + A = obj.A; + B = obj.B; + C = inv(A)*B; + + % Diagonalize C and use T_C to diagonalize B + [T, ~] = eig(C); + Lambda = T'*B*T; + lambda = diag(Lambda); + + + % Identify in- and outgoing variables + Iplus = lambda > 0; + Iminus = lambda < 0; + + switch boundary + case 'l' + Iin = Iplus; + case 'r' + Iin = Iminus; + end + Tin = T(:,Iin); + + switch type + case 'p' + L = [1, 0]; + case 'v' + L = [0, 1]; + case 'characteristic' + % Diagonalize C + [T_C, Lambda_C] = eig(C); + lambda_C = diag(Lambda_C); + % Identify in- and outgoing characteristic variables + Iplus = lambda_C > 0; + Iminus = lambda_C < 0; + + switch boundary + case 'l' + Iin_C = Iplus; + case 'r' + Iin_C = Iminus; + end + T_C_inv = inv(T_C); + L = T_C_inv(Iin_C,:); + otherwise + error('Boundary condition not implemented.'); + end + + % Penalty parameters + sigma = [0; 0]; + sigma(Iin) = lambda(Iin); + + % Sparsify + A = sparse(A); + T = sparse(T); + sigma = sparse(sigma); + L = sparse(L); + Tin = sparse(Tin); + + switch boundary + case 'l' + tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin); + closure = obj.Hi*tau*L*obj.e_l'; + + case 'r' + tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin); + closure = obj.Hi*tau*L*obj.e_r'; + + end + + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + error('Staggered1DAcoustics, interface not implemented'); + + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; + penalty = -obj.Hi*tau*neighbour_scheme.e_r'; + case {'r','right'} + tau = 0*obj.e_r; + closure = obj.Hi*tau*obj.e_r'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Staggered1DAcousticsVariable.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,254 @@ +classdef Staggered1DAcousticsVariable < scheme.Scheme + properties + m % Number of points of primal grid in each direction, possibly a vector + h % Grid spacing + + % Grids + grid % Total grid object + grid_primal + grid_dual + + order % Order accuracy for the approximation + + H % Combined norm + Hi % Inverse + D % Semi-discrete approximation matrix + + D1_primal + D1_dual + + % Pick out left or right boundary + e_l + e_r + + % Initial data + v0 + + % Pick out primal or dual component + e_primal + e_dual + + % System matrices, cell matrices of function handles + A + B + + % Variable coefficient matrices evaluated at boundaries. + A_l, A_r + B_l, B_r + end + + + methods + % Scheme for A*u_t + B u_x = 0, + % u = [p, v]; + % A: Diagonal and A > 0, + % A = [a1, 0; + % 0, a2] + % + % B = B^T and with diagonal entries = 0. + % B = [0 b + % b 0] + % Here we store p on the primal grid and v on the dual + % A and B are cell matrices of function handles + function obj = Staggered1DAcousticsVariable(g, order, A, B) + default_arg('B',{@(x)0*x, @(x)0*x+1; @(x)0*x+1, @(x)0*x}); + default_arg('A',{@(x)0*x+1, @(x)0*x; @(x)0*x, @(x)0*x+1}); + + obj.order = order; + obj.A = A; + obj.B = B; + + % Grids + obj.m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + obj.grid = g; + obj.grid_primal = g.grids{1}; + obj.grid_dual = g.grids{2}; + x_primal = obj.grid_primal.points(); + x_dual = obj.grid_dual.points(); + + % Boundary matrices + obj.A_l = [A{1,1}(xl), A{1,2}(xl);.... + A{2,1}(xl), A{2,2}(xl)]; + obj.A_r = [A{1,1}(xr), A{1,2}(xr);.... + A{2,1}(xr), A{2,2}(xr)]; + obj.B_l = [B{1,1}(xl), B{1,2}(xl);.... + B{2,1}(xl), B{2,2}(xl)]; + obj.B_r = [B{1,1}(xr), B{1,2}(xr);.... + B{2,1}(xr), B{2,2}(xr)]; + + % Get operators + ops = sbp.D1StaggeredUpwind(obj.m, xlim, order); + obj.h = ops.h; + + % Build combined operators + H_primal = ops.H_primal; + H_dual = ops.H_dual; + obj.H = blockmatrix.toMatrix( {H_primal, []; [], H_dual } ); + obj.Hi = inv(obj.H); + + D1_primal = ops.D1_primal; + D1_dual = ops.D1_dual; + A11_B12 = spdiag(-1./A{1,1}(x_primal).*B{1,2}(x_primal), 0); + A22_B21 = spdiag(-1./A{2,2}(x_dual).*B{2,1}(x_dual), 0); + D = {[], A11_B12*D1_primal;... + A22_B21*D1_dual, []}; + obj.D = blockmatrix.toMatrix(D); + obj.D1_primal = D1_primal; + obj.D1_dual = D1_dual; + + % Combined boundary operators + e_primal_l = ops.e_primal_l; + e_primal_r = ops.e_primal_r; + e_dual_l = ops.e_dual_l; + e_dual_r = ops.e_dual_r; + e_l = {e_primal_l, [];... + [] , e_dual_l}; + e_r = {e_primal_r, [];... + [] , e_dual_r}; + obj.e_l = blockmatrix.toMatrix(e_l); + obj.e_r = blockmatrix.toMatrix(e_r); + + % Pick out first or second component of solution + N_primal = obj.grid_primal.N(); + N_dual = obj.grid_dual.N(); + obj.e_primal = [speye(N_primal, N_primal); sparse(N_dual, N_primal)]; + obj.e_dual = [sparse(N_primal, N_dual); speye(N_dual, N_dual)]; + + + end + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other domain. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary, type) + + default_arg('type','p'); + + % type = 'p' => boundary condition for p + % type = 'v' => boundary condition for v + % No other types implemented yet + + % BC on the form Lu - g = 0; + + % Need a transformation T such that w = T^{−1}*u is a + % meaningful change of variables and T^T*B*T is block-diagonal, + % For linear acoustics, T = T_C meets these criteria. + + % Get boundary matrices + switch boundary + case 'l' + A = obj.A_l; + B = obj.B_l; + case 'r' + A = obj.A_r; + B = obj.B_r; + end + C = inv(A)*B; + + % Diagonalize C and use T_C to diagonalize B + [T, ~] = eig(C); + Lambda = T'*B*T; + lambda = diag(Lambda); + + % Identify in- and outgoing characteristic variables + Iplus = lambda > 0; + Iminus = lambda < 0; + + switch boundary + case 'l' + Iin = Iplus; + case 'r' + Iin = Iminus; + end + Tin = T(:,Iin); + + switch type + case 'p' + L = [1, 0]; + case 'v' + L = [0, 1]; + case 'characteristic' + % Diagonalize C + [T_C, Lambda_C] = eig(C); + lambda_C = diag(Lambda_C); + % Identify in- and outgoing characteristic variables + Iplus = lambda_C > 0; + Iminus = lambda_C < 0; + + switch boundary + case 'l' + Iin_C = Iplus; + case 'r' + Iin_C = Iminus; + end + T_C_inv = inv(T_C); + L = T_C_inv(Iin_C,:); + otherwise + error('Boundary condition not implemented.'); + end + + % Penalty parameters + sigma = [0; 0]; + sigma(Iin) = lambda(Iin); + + % Sparsify + A = sparse(A); + T = sparse(T); + sigma = sparse(sigma); + L = sparse(L); + Tin = sparse(Tin); + + switch boundary + case 'l' + tau = -1*obj.e_l * inv(A) * inv(T)' * sigma * inv(L*Tin); + closure = obj.Hi*tau*L*obj.e_l'; + + case 'r' + tau = 1*obj.e_r * inv(A) * inv(T)' * sigma * inv(L*Tin); + closure = obj.Hi*tau*L*obj.e_r'; + + end + + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + error('Staggered1DAcoustics, interface not implemented'); + + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; + penalty = -obj.Hi*tau*neighbour_scheme.e_r'; + case {'r','right'} + tau = 0*obj.e_r; + closure = obj.Hi*tau*obj.e_r'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
--- a/+scheme/Utux.m Mon Jun 04 17:48:19 2018 +0200 +++ b/+scheme/Utux.m Sun Jun 17 15:29:46 2018 -0700 @@ -2,7 +2,7 @@ properties m % Number of points in each direction, possibly a vector h % Grid spacing - x % Grid + grid % Grid order % Order accuracy for the approximation H % Discrete norm @@ -17,41 +17,40 @@ methods - function obj = Utux(m,xlim,order,operator) - default_arg('a',1); - - %Old operators - % [x, h] = util.get_grid(xlim{:},m); - %ops = sbp.Ordinary(m,h,order); - - + function obj = Utux(g ,order, operator) + default_arg('operator','Standard'); + + m = g.size(); + xl = g.getBoundary('l'); + xr = g.getBoundary('r'); + xlim = {xl, xr}; + switch operator - case 'NonEquidistant' - ops = sbp.D1Nonequidistant(m,xlim,order); - obj.D1 = ops.D1; +% case 'NonEquidistant' +% ops = sbp.D1Nonequidistant(m,xlim,order); +% obj.D1 = ops.D1; case 'Standard' ops = sbp.D2Standard(m,xlim,order); obj.D1 = ops.D1; - case 'Upwind' - ops = sbp.D1Upwind(m,xlim,order); - obj.D1 = ops.Dm; +% case 'Upwind' +% ops = sbp.D1Upwind(m,xlim,order); +% obj.D1 = ops.Dm; otherwise error('Unvalid operator') end - obj.x=ops.x; + + obj.grid = g; - obj.H = ops.H; obj.Hi = ops.HI; obj.e_l = ops.e_l; obj.e_r = ops.e_r; - obj.D=obj.D1; + obj.D = -obj.D1; obj.m = m; obj.h = ops.h; obj.order = order; - obj.x = ops.x; end % Closure functions return the opertors applied to the own doamin to close the boundary @@ -61,17 +60,27 @@ % data is a function returning the data that should be applied at the boundary. % neighbour_scheme is an instance of Scheme that should be interfaced to. % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj,boundary,type,data) - default_arg('type','neumann'); - default_arg('data',0); + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); tau =-1*obj.e_l; closure = obj.Hi*tau*obj.e_l'; - penalty = 0*obj.e_l; + penalty = -obj.Hi*tau; end function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) - error('An interface function does not exist yet'); + switch boundary + % Upwind coupling + case {'l','left'} + tau = -1*obj.e_l; + closure = obj.Hi*tau*obj.e_l'; + penalty = -obj.Hi*tau*neighbour_scheme.e_r'; + case {'r','right'} + tau = 0*obj.e_r; + closure = obj.Hi*tau*obj.e_r'; + penalty = -obj.Hi*tau*neighbour_scheme.e_l'; + end + end function N = size(obj) @@ -81,9 +90,9 @@ end methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Utux2D.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,286 @@ +classdef Utux2D < scheme.Scheme + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + grid % Grid + order % Order accuracy for the approximation + v0 % Initial data + + a % Wave speed a = [a1, a2]; + % Can either be a constant vector or a cell array of function handles. + + H % Discrete norm + H_x, H_y % Norms in the x and y directions + Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms + + % Derivatives + Dx, Dy + + % Boundary operators + e_w, e_e, e_s, e_n + + D % Total discrete operator + + % String, type of interface coupling + % Default: 'upwind' + % Other: 'centered' + coupling_type + + % String, type of interpolation operators + % Default: 'AWW' (Almquist Wang Werpers) + % Other: 'MC' (Mattsson Carpenter) + interpolation_type + + + % Cell array, damping on upwstream and downstream sides. + interpolation_damping + + end + + + methods + function obj = Utux2D(g ,order, opSet, a, coupling_type, interpolation_type, interpolation_damping) + + default_arg('interpolation_damping',{0,0}); + default_arg('interpolation_type','AWW'); + default_arg('coupling_type','upwind'); + default_arg('a',1/sqrt(2)*[1, 1]); + default_arg('opSet',@sbp.D2Standard); + + assert(isa(g, 'grid.Cartesian')) + if iscell(a) + a1 = grid.evalOn(g, a{1}); + a2 = grid.evalOn(g, a{2}); + a = {spdiag(a1), spdiag(a2)}; + else + a = {a(1), a(2)}; + end + + m = g.size(); + m_x = m(1); + m_y = m(2); + m_tot = g.N(); + + xlim = {g.x{1}(1), g.x{1}(end)}; + ylim = {g.x{2}(1), g.x{2}(end)}; + obj.grid = g; + + % Operator sets + ops_x = opSet(m_x, xlim, order); + ops_y = opSet(m_y, ylim, order); + Ix = speye(m_x); + Iy = speye(m_y); + + % Norms + Hx = ops_x.H; + Hy = ops_y.H; + Hxi = ops_x.HI; + Hyi = ops_y.HI; + + obj.H_x = Hx; + obj.H_y = Hy; + obj.H = kron(Hx,Hy); + obj.Hi = kron(Hxi,Hyi); + obj.Hx = kron(Hx,Iy); + obj.Hy = kron(Ix,Hy); + obj.Hxi = kron(Hxi,Iy); + obj.Hyi = kron(Ix,Hyi); + + % Derivatives + Dx = ops_x.D1; + Dy = ops_y.D1; + obj.Dx = kron(Dx,Iy); + obj.Dy = kron(Ix,Dy); + + % Boundary operators + obj.e_w = kr(ops_x.e_l, Iy); + obj.e_e = kr(ops_x.e_r, Iy); + obj.e_s = kr(Ix, ops_y.e_l); + obj.e_n = kr(Ix, ops_y.e_r); + + obj.m = m; + obj.h = [ops_x.h ops_y.h]; + obj.order = order; + obj.a = a; + obj.coupling_type = coupling_type; + obj.interpolation_type = interpolation_type; + obj.interpolation_damping = interpolation_damping; + obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); + + end + % Closure functions return the opertors applied to the own domain to close the boundary + % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition if there are several. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj,boundary,type) + default_arg('type','dirichlet'); + + sigma = -1; % Scalar penalty parameter + switch boundary + case {'w','W','west','West'} + tau = sigma*obj.a{1}*obj.e_w*obj.H_y; + closure = obj.Hi*tau*obj.e_w'; + + case {'s','S','south','South'} + tau = sigma*obj.a{2}*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + end + penalty = -obj.Hi*tau; + + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + + % Get neighbour boundary operator + switch neighbour_boundary + case {'e','E','east','East'} + e_neighbour = neighbour_scheme.e_e; + m_neighbour = neighbour_scheme.m(2); + case {'w','W','west','West'} + e_neighbour = neighbour_scheme.e_w; + m_neighbour = neighbour_scheme.m(2); + case {'n','N','north','North'} + e_neighbour = neighbour_scheme.e_n; + m_neighbour = neighbour_scheme.m(1); + case {'s','S','south','South'} + e_neighbour = neighbour_scheme.e_s; + m_neighbour = neighbour_scheme.m(1); + end + + switch obj.coupling_type + + % Upwind coupling (energy dissipation) + case 'upwind' + sigma_ds = -1; %"Downstream" penalty + sigma_us = 0; %"Upstream" penalty + + % Energy-preserving coupling (no energy dissipation) + case 'centered' + sigma_ds = -1/2; %"Downstream" penalty + sigma_us = 1/2; %"Upstream" penalty + + otherwise + error(['Interface coupling type ' coupling_type ' is not available.']) + end + + % Check grid ratio for interpolation + switch boundary + case {'w','W','west','West','e','E','east','East'} + m = obj.m(2); + case {'s','S','south','South','n','N','north','North'} + m = obj.m(1); + end + grid_ratio = m/m_neighbour; + if grid_ratio ~= 1 + + [ms, index] = sort([m, m_neighbour]); + orders = [obj.order, neighbour_scheme.order]; + orders = orders(index); + + switch obj.interpolation_type + case 'MC' + interpOpSet = sbp.InterpMC(ms(1),ms(2),orders(1),orders(2)); + if grid_ratio < 1 + I_neighbour2local_us = interpOpSet.IF2C; + I_neighbour2local_ds = interpOpSet.IF2C; + I_local2neighbour_us = interpOpSet.IC2F; + I_local2neighbour_ds = interpOpSet.IC2F; + elseif grid_ratio > 1 + I_neighbour2local_us = interpOpSet.IC2F; + I_neighbour2local_ds = interpOpSet.IC2F; + I_local2neighbour_us = interpOpSet.IF2C; + I_local2neighbour_ds = interpOpSet.IF2C; + end + case 'AWW' + %String 'C2F' indicates that ICF2 is more accurate. + interpOpSetF2C = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'F2C'); + interpOpSetC2F = sbp.InterpAWW(ms(1),ms(2),orders(1),orders(2),'C2F'); + if grid_ratio < 1 + % Local is coarser than neighbour + I_neighbour2local_us = interpOpSetC2F.IF2C; + I_neighbour2local_ds = interpOpSetF2C.IF2C; + I_local2neighbour_us = interpOpSetC2F.IC2F; + I_local2neighbour_ds = interpOpSetF2C.IC2F; + elseif grid_ratio > 1 + % Local is finer than neighbour + I_neighbour2local_us = interpOpSetF2C.IC2F; + I_neighbour2local_ds = interpOpSetC2F.IC2F; + I_local2neighbour_us = interpOpSetF2C.IF2C; + I_local2neighbour_ds = interpOpSetC2F.IF2C; + end + otherwise + error(['Interpolation type ' obj.interpolation_type ... + ' is not available.' ]); + end + + else + % No interpolation required + I_neighbour2local_us = speye(m,m); + I_neighbour2local_ds = speye(m,m); + end + + int_damp_us = obj.interpolation_damping{1}; + int_damp_ds = obj.interpolation_damping{2}; + + I = speye(m,m); + I_back_forth_us = I_neighbour2local_us*I_local2neighbour_us; + I_back_forth_ds = I_neighbour2local_ds*I_local2neighbour_ds; + + + switch boundary + case {'w','W','west','West'} + tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y; + closure = obj.Hi*tau*obj.e_w'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; + + beta = int_damp_ds*obj.a{1}... + *obj.e_w*obj.H_y; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; + case {'e','E','east','East'} + tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; + closure = obj.Hi*tau*obj.e_e'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; + + beta = int_damp_us*obj.a{1}... + *obj.e_e*obj.H_y; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; + case {'s','S','south','South'} + tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; + closure = obj.Hi*tau*obj.e_s'; + penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; + + beta = int_damp_ds*obj.a{2}... + *obj.e_s*obj.H_x; + closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; + case {'n','N','north','North'} + tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x; + closure = obj.Hi*tau*obj.e_n'; + penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; + + beta = int_damp_us*obj.a{2}... + *obj.e_n*obj.H_x; + closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; + end + + + end + + function N = size(obj) + N = obj.m; + end + + end + + methods(Static) + % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u + % and bound_v of scheme schm_v. + % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') + function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) + [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); + [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); + end + end +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rkparameters/rk4.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,12 @@ +function [a,b,c,s] = rk4() + +% Butcher tableau for classical RK$ +s = 4; +a = sparse(s,s); +a(2,1) = 1/2; +a(3,2) = 1/2; +a(4,3) = 1; +b = 1/6*[1; 2; 2; 1]; +c = [0; 1/2; 1/2; 1]; + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/ExplicitRungeKuttaDiscreteData.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,116 @@ +classdef ExplicitRungeKuttaDiscreteData < time.Timestepper + properties + D + S % Function handle for time-dependent data + data % Matrix of data vectors, one column per stage + F + k + t + v + m + n + order + a, b, c, s % Butcher tableau + K % Stage rates + U % Stage approximations + T % Stage times + end + + + methods + function obj = ExplicitRungeKuttaDiscreteData(D, S, data, k, t0, v0, order) + default_arg('order', 4); + default_arg('S', []); + default_arg('data', []); + + obj.D = D; + obj.S = S; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.m = length(v0); + obj.n = 0; + obj.order = order; + obj.data = data; + + switch order + case 4 + [obj.a, obj.b, obj.c, obj.s] = time.rkparameters.rk4(); + otherwise + error('That RK method is not available'); + end + + obj.K = sparse(obj.m, obj.s); + obj.U = sparse(obj.m, obj.s); + + end + + function [v,t,U,T,K] = getV(obj) + v = obj.v; + t = obj.t; + U = obj.U; % Stage approximations in previous time step. + T = obj.T; % Stage times in previous time step. + K = obj.K; % Stage rates in previous time step. + end + + function [a,b,c,s] = getTableau(obj) + a = obj.a; + b = obj.b; + c = obj.c; + s = obj.s; + end + + function obj = step(obj) + v = obj.v; + a = obj.a; + b = obj.b; + c = obj.c; + s = obj.s; + S = obj.S; + dt = obj.k; + K = obj.K; + U = obj.U; + D = obj.D; + data = obj.data; + + for i = 1:s + U(:,i) = v; + for j = 1:i-1 + U(:,i) = U(:,i) + dt*a(i,j)*K(:,j); + end + + K(:,i) = D*U(:,i); + obj.T(i) = obj.t + c(i)*dt; + + % Data from continuos function and discrete time-points. + if ~isempty(S) + K(:,i) = K(:,i) + S(obj.T(i)); + end + if ~isempty(data) + K(:,i) = K(:,i) + data(:,obj.n*s + i); + end + + end + + obj.v = v + dt*K*b; + obj.t = obj.t + dt; + obj.n = obj.n + 1; + obj.U = U; + obj.K = K; + end + end + + + methods (Static) + function k = getTimeStep(lambda) + + switch obj.order + case 4 + k = rk4.get_rk4_time_step(lambda); + otherwise + error('Time-step function not available for this order'); + end + end + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscr.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,97 @@ +function ret = diracDiscr(x_0in , x , m_order, s_order, H, h) + +m = length(x); + +% Return zeros if x0 is outside grid +if(x_0in < x(1) || x_0in > x(end) ) + + ret = zeros(size(x)); + +else + + fnorm = diag(H); + eta = abs(x-x_0in); + tot = m_order+s_order; + S = []; + M = []; + + % Get interior grid spacing + middle = floor(m/2); + h = x(middle+1) - x(middle); + + poss = find(tot*h/2 >= eta); + + % Ensure that poss is not too long + if length(poss) == (tot + 2) + poss = poss(2:end-1); + elseif length(poss) == (tot + 1) + poss = poss(1:end-1); + end + + % Use first tot grid points + if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h; + index=1:tot; + pol=(x(1:tot)-x(1))/(x(tot)-x(1)); + x_0=(x_0in-x(1))/(x(tot)-x(1)); + norm=fnorm(1:tot)/h; + + % Use last tot grid points + elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h; + index = length(x)-tot+1:length(x); + pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1)); + norm = fnorm(end-tot+1:end)/h; + x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1)); + + % Interior, compensate for round-off errors. + elseif length(poss) < tot + if poss(end)<m + poss = [poss; poss(end)+1]; + else + poss = [poss(1)-1; poss]; + end + pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); + x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); + norm = fnorm(poss)/h; + index = poss; + + % Interior + else + pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); + x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); + norm = fnorm(poss)/h; + index = poss; + end + + h_pol = pol(2)-pol(1); + b = zeros(m_order+s_order,1); + + for i = 1:m_order + b(i,1) = x_0^(i-1); + end + + for i = 1:(m_order+s_order) + for j = 1:m_order + M(j,i) = pol(i)^(j-1)*h_pol*norm(i); + end + end + + for i = 1:(m_order+s_order) + for j = 1:s_order + S(j,i) = (-1)^(i-1)*pol(i)^(j-1); + end + end + + A = [M;S]; + + d = A\b; + ret = x*0; + ret(index) = d/h*h_pol; +end + +end + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscrTest.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,379 @@ +function tests = diracDiscrTest() + tests = functiontests(localfunctions); +end + +function testLeftGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test left boundary grid points + x0s = xl + [0, h, 2*h]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testLeftRandom(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test random points near left boundary + x0s = xl + 2*h*rand(1,10); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testRightGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test right boundary grid points + x0s = xr-[0, h, 2*h]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testRightRandom(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test random points near right boundary + x0s = xr - 2*h*rand(1,10); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testInteriorGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test interior grid points + m_half = round(m/2); + x0s = xl + (m_half-1:m_half+1)*h; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testInteriorRandom(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test random points in interior + x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +% x0 outside grid should yield 0 integral! +function testX0OutsideGrid(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test points outisde grid + x0s = [xl-1.1*h, xr+1.1*h]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - 0); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test all grid points + x0s = x; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testHalfGP(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + + % Test halfway between all grid points + x0s = 1/2*( x(2:end)+x(1:end-1) ); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllGPStaggered(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + + % Test all grid points + x0s = x; + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testHalfGPStaggered(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + + % Test halfway between all grid points + x0s = 1/2*( x(2:end)+x(1:end-1) ); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testRandomStaggered(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + + % Test random points within grid boundaries + x0s = xl + (xr-xl)*rand(1,300); + + for j = 1:length(fs) + f = fs{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral - f(x0)); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + + +% ============== Setup functions ======================= +function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond) + + % Grid + xl = -3; + xr = 900; + L = xr-xl; + m = 101; + h = (xr-xl)/(m-1); + g = grid.equidistant(m, {xl, xr}); + x = g.points(); + + % Quadrature + ops = sbp.D2Standard(m, {xl, xr}, order); + H = ops.H; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x) (x/L).^p; + end + +end + +function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond) + + % Grid + xl = -3; + xr = 900; + L = xr-xl; + m = 101; + [~, g_dual] = grid.primalDual1D(m, {xl, xr}); + x = g_dual.points(); + h = g_dual.h; + + % Quadrature + ops = sbp.D1Staggered(m, {xl, xr}, order); + H = ops.H_dual; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x) (x/L).^p; + end + +end + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/kroneckerDelta.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,6 @@ +function d = kroneckerDelta(i,j) + +d = 0; +if i==j + d = 1; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagVariable.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,17 @@ +function A = spdiagVariable(a,i) + default_arg('i',0); + + if isrow(a) + a = a'; + end + + n = length(a)+abs(i); + + if i > 0 + a = [sparse(i,1); a]; + elseif i < 0 + a = [a; sparse(abs(i),1)]; + end + + A = spdiags(a,i,n,n); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagsVariablePeriodic.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,42 @@ +function A = spdiagsVariablePeriodic(vals,diags) + % Creates an m x m periodic discretization matrix. + % vals - m x ndiags matrix of values + % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on + % vals that are not on main diagonal are going to spill over to + % off-diagonal corners. + + default_arg('diags',0); + + [m, ~] = size(vals); + + A = sparse(m,m); + + for i = 1:length(diags) + + d = diags(i); + a = vals(:,i); + + % Sub-diagonals + if d < 0 + a_bulk = a(1+abs(d):end); + a_corner = a(1:1+abs(d)-1); + corner_diag = m-abs(d); + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Super-diagonals + elseif d > 0 + a_bulk = a(1:end-d); + a_corner = a(end-d+1:end); + corner_diag = -m + d; + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Main diagonal + else + A = A + spdiagVariable(a, 0); + end + + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stripeMatrixPeriodic.m Sun Jun 17 15:29:46 2018 -0700 @@ -0,0 +1,8 @@ +% Creates a periodic discretization matrix of size n x n +% with the values of val on the diagonals diag. +% A = stripeMatrix(val,diags,n) +function A = stripeMatrixPeriodic(val,diags,n) + + D = ones(n,1)*val; + A = spdiagsVariablePeriodic(D,diags); +end \ No newline at end of file