Mercurial > repos > public > sbplib
changeset 1131:ea225a4659fe feature/laplace_curvilinear_test
Add 2d tests for derivative of delta function, and confirm that they work.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 22 May 2019 15:31:26 -0700 |
| parents | 99fd66ffe714 |
| children | 375f73edbbd4 |
| files | diracPrimDiscrTest.m |
| diffstat | 1 files changed, 107 insertions(+), 84 deletions(-) [+] |
line wrap: on
line diff
--- a/diracPrimDiscrTest.m Tue May 21 18:44:01 2019 -0700 +++ b/diracPrimDiscrTest.m Wed May 22 15:31:26 2019 -0700 @@ -258,74 +258,88 @@ %=============== 2D tests ============================== -% function testAllGP2D(testCase) +function testAllGP2D(testCase) -% orders = [2, 4, 6]; -% mom_conds = orders; + 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}); + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, m, x, X, ~, H, fs, dfdxs, dfdys] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); -% % Test all grid points -% x0s = X; + % 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 + for j = 1:length(fs) + f = fs{j}; + dfdx = dfdxs{j}; + dfdy = dfdys{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + + delta_x = diracPrimDiscr(x0, x, mom_cond, 0, H, 1); + integral1 = delta_x'*H_global*fx; -% function testAllRandom2D(testCase) + delta_y = diracPrimDiscr(x0, x, mom_cond, 0, H, 2); + integral2 = delta_y'*H_global*fx; + + err = abs(integral1 + dfdx(x0(1), x0(2))) + abs(integral2 + dfdy(x0(1), x0(2))); -% orders = [2, 4, 6]; -% mom_conds = orders; + 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}); + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, m, x, X, h, H, fs, dfdxs, dfdys] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); -% xl = xlims{1}; -% xr = xlims{2}; -% yl = ylims{1}; -% yr = ylims{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) ]; + % 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; + for j = 1:length(fs) + f = fs{j}; + dfdx = dfdxs{j}; + dfdy = dfdys{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + + delta_x = diracPrimDiscr(x0, x, mom_cond, 0, H, 1); + integral1 = delta_x'*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 + delta_y = diracPrimDiscr(x0, x, mom_cond, 0, H, 2); + integral2 = delta_y'*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(integral1 - 0) + abs(integral2 - 0); + else + err = abs(integral1 + dfdx(x0(1), x0(2))) + abs(integral2 + dfdy(x0(1), x0(2))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end %=============== 3D tests ============================== % function testAllGP3D(testCase) @@ -434,39 +448,48 @@ end -% function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond) +function [xlims, ylims, m, x, X, h, H, fs, fxs, fys] = setup2D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; -% % 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; -% 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}; -% % 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,1); + fxs = cell(mom_cond+1,1); + fys = cell(mom_cond+1,1); + for p = 0:mom_cond + fs{p+1} = @(x,y) (x/Lx + y/Ly).^p; + end -% % Moment conditions -% fs = cell(mom_cond,1); -% for p = 0:mom_cond-1 -% fs{p+1} = @(x,y) (x/Lx + y/Ly).^p; -% end + fxs{1} = @(x,y) 0*x; + fys{1} = @(x,y) 0*y; + for p = 1:mom_cond + fxs{p+1} = @(x,y) p/Lx*(x/Lx + y/Ly).^(p-1); + fys{p+1} = @(x,y) p/Ly*(x/Lx + y/Ly).^(p-1); + 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}; + % 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 +end % function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond)