Mercurial > repos > public > sbplib

comparison diracPrimDiscrTest.m @ 1131:ea225a4659fe feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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
comparison
equal deleted inserted replaced
1130:99fd66ffe714 1131:ea225a4659fe
256 end 256 end
257 end 257 end
258 258
259 259
260 %=============== 2D tests ============================== 260 %=============== 2D tests ==============================
261 % function testAllGP2D(testCase) 261 function testAllGP2D(testCase)
262 262
263 % orders = [2, 4, 6]; 263 orders = [2, 4, 6];
264 % mom_conds = orders; 264 mom_conds = orders;
265 265
266 % for o = 1:length(orders) 266 for o = 1:length(orders)
267 % order = orders(o); 267 order = orders(o);
268 % mom_cond = mom_conds(o); 268 mom_cond = mom_conds(o);
269 % [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond); 269 [xlims, ylims, m, x, X, ~, H, fs, dfdxs, dfdys] = setup2D(order, mom_cond);
270 % H_global = kron(H{1}, H{2}); 270 H_global = kron(H{1}, H{2});
271 271
272 % % Test all grid points 272 % Test all grid points
273 % x0s = X; 273 x0s = X;
274 274
275 % for j = 1:length(fs) 275 for j = 1:length(fs)
276 % f = fs{j}; 276 f = fs{j};
277 % fx = f(X(:,1), X(:,2)); 277 dfdx = dfdxs{j};
278 % for i = 1:length(x0s) 278 dfdy = dfdys{j};
279 % x0 = x0s(i,:); 279 fx = f(X(:,1), X(:,2));
280 % delta = diracPrimDiscr(x0, x, mom_cond, 0, H); 280 for i = 1:length(x0s)
281 % integral = delta'*H_global*fx; 281 x0 = x0s(i,:);
282 % err = abs(integral - f(x0(1), x0(2))); 282
283 % testCase.verifyLessThan(err, 1e-12); 283 delta_x = diracPrimDiscr(x0, x, mom_cond, 0, H, 1);
284 % end 284 integral1 = delta_x'*H_global*fx;
285 % end 285
286 % end 286 delta_y = diracPrimDiscr(x0, x, mom_cond, 0, H, 2);
287 % end 287 integral2 = delta_y'*H_global*fx;
288 288
289 % function testAllRandom2D(testCase) 289 err = abs(integral1 + dfdx(x0(1), x0(2))) + abs(integral2 + dfdy(x0(1), x0(2)));
290 290
291 % orders = [2, 4, 6]; 291 testCase.verifyLessThan(err, 1e-12);
292 % mom_conds = orders; 292 end
293 293 end
294 % for o = 1:length(orders) 294 end
295 % order = orders(o); 295 end
296 % mom_cond = mom_conds(o); 296
297 % [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond); 297 function testAllRandom2D(testCase)
298 % H_global = kron(H{1}, H{2}); 298
299 299 orders = [2, 4, 6];
300 % xl = xlims{1}; 300 mom_conds = orders;
301 % xr = xlims{2}; 301
302 % yl = ylims{1}; 302 for o = 1:length(orders)
303 % yr = ylims{2}; 303 order = orders(o);
304 304 mom_cond = mom_conds(o);
305 % % Test random points, even outside grid 305 [xlims, ylims, m, x, X, h, H, fs, dfdxs, dfdys] = setup2D(order, mom_cond);
306 % Npoints = 100; 306 H_global = kron(H{1}, H{2});
307 % x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ... 307
308 % (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ]; 308 xl = xlims{1};
309 309 xr = xlims{2};
310 % for j = 1:length(fs) 310 yl = ylims{1};
311 % f = fs{j}; 311 yr = ylims{2};
312 % fx = f(X(:,1), X(:,2)); 312
313 % for i = 1:length(x0s) 313 % Test random points, even outside grid
314 % x0 = x0s(i,:); 314 Npoints = 100;
315 % delta = diracPrimDiscr(x0, x, mom_cond, 0, H); 315 x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ...
316 % integral = delta'*H_global*fx; 316 (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ];
317 317
318 % % Integral should be 0 if point is outside grid 318 for j = 1:length(fs)
319 % if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr 319 f = fs{j};
320 % err = abs(integral - 0); 320 dfdx = dfdxs{j};
321 % else 321 dfdy = dfdys{j};
322 % err = abs(integral - f(x0(1), x0(2))); 322 fx = f(X(:,1), X(:,2));
323 % end 323 for i = 1:length(x0s)
324 % testCase.verifyLessThan(err, 1e-12); 324 x0 = x0s(i,:);
325 % end 325
326 % end 326 delta_x = diracPrimDiscr(x0, x, mom_cond, 0, H, 1);
327 % end 327 integral1 = delta_x'*H_global*fx;
328 % end 328
329 delta_y = diracPrimDiscr(x0, x, mom_cond, 0, H, 2);
330 integral2 = delta_y'*H_global*fx;
331
332 % Integral should be 0 if point is outside grid
333 if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr
334 err = abs(integral1 - 0) + abs(integral2 - 0);
335 else
336 err = abs(integral1 + dfdx(x0(1), x0(2))) + abs(integral2 + dfdy(x0(1), x0(2)));
337 end
338 testCase.verifyLessThan(err, 1e-12);
339 end
340 end
341 end
342 end
329 343
330 %=============== 3D tests ============================== 344 %=============== 3D tests ==============================
331 % function testAllGP3D(testCase) 345 % function testAllGP3D(testCase)
332 346
333 % orders = [2, 4, 6]; 347 % orders = [2, 4, 6];
432 fps{p+1} = @(x) p/L*(x/L).^(p-1); 446 fps{p+1} = @(x) p/L*(x/L).^(p-1);
433 end 447 end
434 448
435 end 449 end
436 450
437 % function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond) 451 function [xlims, ylims, m, x, X, h, H, fs, fxs, fys] = setup2D(order, mom_cond)
438 452
439 % % Grid 453 % Grid
440 % xlims = {-3, 20}; 454 xlims = {-3, 20};
441 % ylims = {-11,5}; 455 ylims = {-11,5};
442 % Lx = xlims{2} - xlims{1}; 456 Lx = xlims{2} - xlims{1};
443 % Ly = ylims{2} - ylims{1}; 457 Ly = ylims{2} - ylims{1};
444 458
445 % m = [15, 16]; 459 m = [15, 16];
446 % g = grid.equidistant(m, xlims, ylims); 460 g = grid.equidistant(m, xlims, ylims);
447 % X = g.points(); 461 X = g.points();
448 % x = g.x; 462 x = g.x;
449 463
450 % % Quadrature 464 % Quadrature
451 % opsx = sbp.D2Standard(m(1), xlims, order); 465 opsx = sbp.D2Standard(m(1), xlims, order);
452 % opsy = sbp.D2Standard(m(2), ylims, order); 466 opsy = sbp.D2Standard(m(2), ylims, order);
453 % Hx = opsx.H; 467 Hx = opsx.H;
454 % Hy = opsy.H; 468 Hy = opsy.H;
455 % H = {Hx, Hy}; 469 H = {Hx, Hy};
456 470
457 % % Moment conditions 471 % Moment conditions
458 % fs = cell(mom_cond,1); 472 fs = cell(mom_cond+1,1);
459 % for p = 0:mom_cond-1 473 fxs = cell(mom_cond+1,1);
460 % fs{p+1} = @(x,y) (x/Lx + y/Ly).^p; 474 fys = cell(mom_cond+1,1);
461 % end 475 for p = 0:mom_cond
462 476 fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
463 % % Grid spacing in interior 477 end
464 % mm = round(m/2); 478
465 % hx = x{1}(mm(1)+1) - x{1}(mm(1)); 479 fxs{1} = @(x,y) 0*x;
466 % hy = x{2}(mm(2)+1) - x{2}(mm(2)); 480 fys{1} = @(x,y) 0*y;
467 % h = {hx, hy}; 481 for p = 1:mom_cond
468 482 fxs{p+1} = @(x,y) p/Lx*(x/Lx + y/Ly).^(p-1);
469 % end 483 fys{p+1} = @(x,y) p/Ly*(x/Lx + y/Ly).^(p-1);
484 end
485
486 % Grid spacing in interior
487 mm = round(m/2);
488 hx = x{1}(mm(1)+1) - x{1}(mm(1));
489 hy = x{2}(mm(2)+1) - x{2}(mm(2));
490 h = {hx, hy};
491
492 end
470 493
471 % function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond) 494 % function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond)
472 495
473 % % Grid 496 % % Grid
474 % xlims = {-3, 20}; 497 % xlims = {-3, 20};