Mercurial > repos > public > sbplib
changeset 1130:99fd66ffe714 feature/laplace_curvilinear_test
Add derivative of delta functions and corresponding tests, tested for 1D.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 21 May 2019 18:44:01 -0700 |
| parents | b29892853daf |
| children | ea225a4659fe |
| files | diracPrimDiscr.m diracPrimDiscr1D.m diracPrimDiscrTest.m |
| diffstat | 3 files changed, 638 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracPrimDiscr.m Tue May 21 18:44:01 2019 -0700 @@ -0,0 +1,34 @@ + +function d = diracPrimDiscr(x_s, x, m_order, s_order, H, derivDir) + % n-dimensional delta function, with derivative in direction number derivDir + % 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 + default_arg('derivDir', 1); + + dim = length(x_s); + d_1D = cell(dim,1); + + % If 1D, non-cell input is accepted + if dim == 1 && ~iscell(x) + d = diracPrimDiscr1D(x_s, x, m_order, s_order, H); + else + for i = 1:dim + if i == derivDir + d_1D{i} = diracPrimDiscr1D(x_s(i), x{i}, m_order, s_order, H{i}); + else + d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i}); + end + 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 \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracPrimDiscr1D.m Tue May 21 18:44:01 2019 -0700 @@ -0,0 +1,96 @@ +% Generates discretized derivative of delta function in 1D +function ret = diracPrimDiscr1D(x_0in, x, m_order, s_order, H) + + % diracPrim satisfies one more moment condition than dirac + m_order = m_order + 1; + + 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); + + b(1) = 0; + for i = 2:m_order + b(i) = -(i-1)*x_0^(i-2); + 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_pol/h)^2; + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracPrimDiscrTest.m Tue May 21 18:44:01 2019 -0700 @@ -0,0 +1,508 @@ +function tests = diracPrimDiscrTest() + 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, fps] = setup1D(order, mom_cond); + + % Test left boundary grid points + x0s = xl + [0, h, 2*h]; + + for j = 1:length(fs) + f = fs{j}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = 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}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = setup1D(order, mom_cond); + + % Test right boundary grid points + x0s = xr-[0, h, 2*h]; + + for j = 1:length(fs) + f = fs{j}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = 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}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = 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}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = 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}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = 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}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(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, fps] = setup1D(order, mom_cond); + + % Test all grid points + x0s = x; + + for j = 1:length(fs) + f = fs{j}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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, fps] = 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}; + fp = fps{j}; + fx = f(x); + for i = 1:length(x0s) + x0 = x0s(i); + delta = diracPrimDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H*fx; + err = abs(integral + fp(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 = diracPrimDiscr(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 = diracPrimDiscr(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 = diracPrimDiscr(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 = diracPrimDiscr(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, fps] = 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,1); + fps = cell(mom_cond+1,1); + + for p = 0:mom_cond + fs{p+1} = @(x) (x/L).^p; + end + + fps{1} = @(x) 0*x; + for p = 1:mom_cond + fps{p+1} = @(x) p/L*(x/L).^(p-1); + 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