Mercurial > repos > public > sbplib
changeset 1229:86ee5648e384 feature/dirac_discr
Add multi-d dirac discretization with tests
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 19 Nov 2019 10:56:57 -0800 |
| parents | 2a3bd78edb0e |
| children | 8a456f6e54cc 52d774e69b1f |
| files | diracDiscr.m diracDiscrTest.m |
| diffstat | 2 files changed, 725 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscr.m Tue Nov 19 10:56:57 2019 -0800 @@ -0,0 +1,130 @@ + +function d = diracDiscr(x_s, x, m_order, s_order, H) + % n-dimensional delta function + % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z]. + % x: cell array of grid point column vectors for each dimension. + % m_order: Number of moment conditions + % s_order: Number of smoothness conditions + % H: cell array of 1D norm matrices + + dim = length(x_s); + d_1D = cell(dim,1); + + % If 1D, non-cell input is accepted + if dim == 1 && ~iscell(x) + d = diracDiscr1D(x_s, x, m_order, s_order, H); + + else + for i = 1:dim + d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i}); + end + + d = d_1D{dim}; + for i = dim-1: -1: 1 + % Perform outer product, transpose, and then turn into column vector + d = (d_1D{i}*d')'; + d = d(:); + end + end + +end + + +% Helper function for 1D delta functions +function ret = diracDiscr1D(x_0in , x , m_order, s_order, 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 Tue Nov 19 10:56:57 2019 -0800 @@ -0,0 +1,595 @@ +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] = setup1D(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] = setup1D(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] = setup1D(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] = setup1D(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] = setup1D(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] = setup1D(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] = setup1D(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] = setup1D(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] = setup1D(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 + +%=============== 2D tests ============================== +function testAllGP2D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + + % Test all grid points + x0s = X; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2))); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllRandom2D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + + xl = xlims{1}; + xr = xlims{2}; + yl = ylims{1}; + yr = ylims{2}; + + % Test random points, even outside grid + Npoints = 100; + x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ... + (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H_global*fx; + + % Integral should be 0 if point is outside grid + if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +%=============== 3D tests ============================== +function testAllGP3D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + + % Test all grid points + x0s = X; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2), X(:,3)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2), x0(3))); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllRandom3D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + + xl = xlims{1}; + xr = xlims{2}; + yl = ylims{1}; + yr = ylims{2}; + zl = zlims{1}; + zr = zlims{2}; + + % Test random points, even outside grid + Npoints = 200; + x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ... + (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1), ... + (zl-3*h{3}) + (zr-zl+6*h{3})*rand(Npoints,1) ]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2), X(:,3)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H_global*fx; + + % Integral should be 0 if point is outside grid + if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2), x0(3))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + + +% ====================================================== +% ============== Setup functions ======================= +% ====================================================== +function [xl, xr, m, h, x, H, fs] = setup1D(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 [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; + + m = [15, 16]; + g = grid.equidistant(m, xlims, ylims); + X = g.points(); + x = g.x; + + % Quadrature + opsx = sbp.D2Standard(m(1), xlims, order); + opsy = sbp.D2Standard(m(2), ylims, order); + Hx = opsx.H; + Hy = opsy.H; + H = {Hx, Hy}; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x,y) (x/Lx + y/Ly).^p; + end + + % Grid spacing in interior + mm = round(m/2); + hx = x{1}(mm(1)+1) - x{1}(mm(1)); + hy = x{2}(mm(2)+1) - x{2}(mm(2)); + h = {hx, hy}; + +end + +function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + zlims = {2,4}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; + Lz = zlims{2} - zlims{1}; + + m = [13, 14, 15]; + g = grid.equidistant(m, xlims, ylims, zlims); + X = g.points(); + x = g.x; + + % Quadrature + opsx = sbp.D2Standard(m(1), xlims, order); + opsy = sbp.D2Standard(m(2), ylims, order); + opsz = sbp.D2Standard(m(3), zlims, order); + Hx = opsx.H; + Hy = opsy.H; + Hz = opsz.H; + H = {Hx, Hy, Hz}; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p; + end + + % Grid spacing in interior + mm = round(m/2); + hx = x{1}(mm(1)+1) - x{1}(mm(1)); + hy = x{2}(mm(2)+1) - x{2}(mm(2)); + hz = x{3}(mm(3)+1) - x{3}(mm(3)); + h = {hx, hy, hz}; + +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