Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 837:43a1c3ac07b1 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make diracDiscr work for n dimensions. Add tests up to 3D.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Sep 2018 18:09:18 -0700
parents 049e4c6fa630
children c70131daaa6e 60c875c18de3
comparison
equal deleted inserted replaced
836:049e4c6fa630 837:43a1c3ac07b1
8 mom_conds = orders; 8 mom_conds = orders;
9 9
10 for o = 1:length(orders) 10 for o = 1:length(orders)
11 order = orders(o); 11 order = orders(o);
12 mom_cond = mom_conds(o); 12 mom_cond = mom_conds(o);
13 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 13 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
14 14
15 % Test left boundary grid points 15 % Test left boundary grid points
16 x0s = xl + [0, h, 2*h]; 16 x0s = xl + [0, h, 2*h];
17 17
18 for j = 1:length(fs) 18 for j = 1:length(fs)
35 mom_conds = orders; 35 mom_conds = orders;
36 36
37 for o = 1:length(orders) 37 for o = 1:length(orders)
38 order = orders(o); 38 order = orders(o);
39 mom_cond = mom_conds(o); 39 mom_cond = mom_conds(o);
40 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 40 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
41 41
42 % Test random points near left boundary 42 % Test random points near left boundary
43 x0s = xl + 2*h*rand(1,10); 43 x0s = xl + 2*h*rand(1,10);
44 44
45 for j = 1:length(fs) 45 for j = 1:length(fs)
62 mom_conds = orders; 62 mom_conds = orders;
63 63
64 for o = 1:length(orders) 64 for o = 1:length(orders)
65 order = orders(o); 65 order = orders(o);
66 mom_cond = mom_conds(o); 66 mom_cond = mom_conds(o);
67 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 67 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
68 68
69 % Test right boundary grid points 69 % Test right boundary grid points
70 x0s = xr-[0, h, 2*h]; 70 x0s = xr-[0, h, 2*h];
71 71
72 for j = 1:length(fs) 72 for j = 1:length(fs)
89 mom_conds = orders; 89 mom_conds = orders;
90 90
91 for o = 1:length(orders) 91 for o = 1:length(orders)
92 order = orders(o); 92 order = orders(o);
93 mom_cond = mom_conds(o); 93 mom_cond = mom_conds(o);
94 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 94 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
95 95
96 % Test random points near right boundary 96 % Test random points near right boundary
97 x0s = xr - 2*h*rand(1,10); 97 x0s = xr - 2*h*rand(1,10);
98 98
99 for j = 1:length(fs) 99 for j = 1:length(fs)
116 mom_conds = orders; 116 mom_conds = orders;
117 117
118 for o = 1:length(orders) 118 for o = 1:length(orders)
119 order = orders(o); 119 order = orders(o);
120 mom_cond = mom_conds(o); 120 mom_cond = mom_conds(o);
121 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 121 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
122 122
123 % Test interior grid points 123 % Test interior grid points
124 m_half = round(m/2); 124 m_half = round(m/2);
125 x0s = xl + (m_half-1:m_half+1)*h; 125 x0s = xl + (m_half-1:m_half+1)*h;
126 126
144 mom_conds = orders; 144 mom_conds = orders;
145 145
146 for o = 1:length(orders) 146 for o = 1:length(orders)
147 order = orders(o); 147 order = orders(o);
148 mom_cond = mom_conds(o); 148 mom_cond = mom_conds(o);
149 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 149 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
150 150
151 % Test random points in interior 151 % Test random points in interior
152 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); 152 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
153 153
154 for j = 1:length(fs) 154 for j = 1:length(fs)
172 mom_conds = orders; 172 mom_conds = orders;
173 173
174 for o = 1:length(orders) 174 for o = 1:length(orders)
175 order = orders(o); 175 order = orders(o);
176 mom_cond = mom_conds(o); 176 mom_cond = mom_conds(o);
177 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 177 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
178 178
179 % Test points outisde grid 179 % Test points outisde grid
180 x0s = [xl-1.1*h, xr+1.1*h]; 180 x0s = [xl-1.1*h, xr+1.1*h];
181 181
182 for j = 1:length(fs) 182 for j = 1:length(fs)
199 mom_conds = orders; 199 mom_conds = orders;
200 200
201 for o = 1:length(orders) 201 for o = 1:length(orders)
202 order = orders(o); 202 order = orders(o);
203 mom_cond = mom_conds(o); 203 mom_cond = mom_conds(o);
204 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 204 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
205 205
206 % Test all grid points 206 % Test all grid points
207 x0s = x; 207 x0s = x;
208 208
209 for j = 1:length(fs) 209 for j = 1:length(fs)
226 mom_conds = orders; 226 mom_conds = orders;
227 227
228 for o = 1:length(orders) 228 for o = 1:length(orders)
229 order = orders(o); 229 order = orders(o);
230 mom_cond = mom_conds(o); 230 mom_cond = mom_conds(o);
231 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); 231 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
232 232
233 % Test halfway between all grid points 233 % Test halfway between all grid points
234 x0s = 1/2*( x(2:end)+x(1:end-1) ); 234 x0s = 1/2*( x(2:end)+x(1:end-1) );
235 235
236 for j = 1:length(fs) 236 for j = 1:length(fs)
326 % end 326 % end
327 % end 327 % end
328 % end 328 % end
329 % end 329 % end
330 330
331 331 %=============== 2D tests ==============================
332 function testAllGP2D(testCase)
333
334 orders = [2, 4, 6];
335 mom_conds = orders;
336
337 for o = 1:length(orders)
338 order = orders(o);
339 mom_cond = mom_conds(o);
340 [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond);
341 H_global = kron(H{1}, H{2});
342
343 % Test all grid points
344 x0s = X;
345
346 for j = 1:length(fs)
347 f = fs{j};
348 fx = f(X(:,1), X(:,2));
349 for i = 1:length(x0s)
350 x0 = x0s(i,:);
351 delta = diracDiscr(x0, x, mom_cond, 0, H);
352 integral = delta'*H_global*fx;
353 err = abs(integral - f(x0(1), x0(2)));
354 testCase.verifyLessThan(err, 1e-12);
355 end
356 end
357 end
358 end
359
360 function testAllRandom2D(testCase)
361
362 orders = [2, 4, 6];
363 mom_conds = orders;
364
365 for o = 1:length(orders)
366 order = orders(o);
367 mom_cond = mom_conds(o);
368 [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond);
369 H_global = kron(H{1}, H{2});
370
371 xl = xlims{1};
372 xr = xlims{2};
373 yl = ylims{1};
374 yr = ylims{2};
375
376 % Test random points, even outside grid
377 Npoints = 100;
378 x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ...
379 (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ];
380
381 for j = 1:length(fs)
382 f = fs{j};
383 fx = f(X(:,1), X(:,2));
384 for i = 1:length(x0s)
385 x0 = x0s(i,:);
386 delta = diracDiscr(x0, x, mom_cond, 0, H);
387 integral = delta'*H_global*fx;
388
389 % Integral should be 0 if point is outside grid
390 if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr
391 err = abs(integral - 0);
392 else
393 err = abs(integral - f(x0(1), x0(2)));
394 end
395 testCase.verifyLessThan(err, 1e-12);
396 end
397 end
398 end
399 end
400
401 %=============== 3D tests ==============================
402 function testAllGP3D(testCase)
403
404 orders = [2, 4, 6];
405 mom_conds = orders;
406
407 for o = 1:length(orders)
408 order = orders(o);
409 mom_cond = mom_conds(o);
410 [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
411 H_global = kron(kron(H{1}, H{2}), H{3});
412
413 % Test all grid points
414 x0s = X;
415
416 for j = 1:length(fs)
417 f = fs{j};
418 fx = f(X(:,1), X(:,2), X(:,3));
419 for i = 1:length(x0s)
420 x0 = x0s(i,:);
421 delta = diracDiscr(x0, x, mom_cond, 0, H);
422 integral = delta'*H_global*fx;
423 err = abs(integral - f(x0(1), x0(2), x0(3)));
424 testCase.verifyLessThan(err, 1e-12);
425 end
426 end
427 end
428 end
429
430 function testAllRandom3D(testCase)
431
432 orders = [2, 4, 6];
433 mom_conds = orders;
434
435 for o = 1:length(orders)
436 order = orders(o);
437 mom_cond = mom_conds(o);
438 [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
439 H_global = kron(kron(H{1}, H{2}), H{3});
440
441 xl = xlims{1};
442 xr = xlims{2};
443 yl = ylims{1};
444 yr = ylims{2};
445 zl = zlims{1};
446 zr = zlims{2};
447
448 % Test random points, even outside grid
449 Npoints = 200;
450 x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ...
451 (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1), ...
452 (zl-3*h{3}) + (zr-zl+6*h{3})*rand(Npoints,1) ];
453
454 for j = 1:length(fs)
455 f = fs{j};
456 fx = f(X(:,1), X(:,2), X(:,3));
457 for i = 1:length(x0s)
458 x0 = x0s(i,:);
459 delta = diracDiscr(x0, x, mom_cond, 0, H);
460 integral = delta'*H_global*fx;
461
462 % Integral should be 0 if point is outside grid
463 if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr
464 err = abs(integral - 0);
465 else
466 err = abs(integral - f(x0(1), x0(2), x0(3)));
467 end
468 testCase.verifyLessThan(err, 1e-12);
469 end
470 end
471 end
472 end
473
474
475 % ======================================================
332 % ============== Setup functions ======================= 476 % ============== Setup functions =======================
333 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond) 477 % ======================================================
478 function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond)
334 479
335 % Grid 480 % Grid
336 xl = -3; 481 xl = -3;
337 xr = 900; 482 xr = 900;
338 L = xr-xl; 483 L = xr-xl;
351 fs{p+1} = @(x) (x/L).^p; 496 fs{p+1} = @(x) (x/L).^p;
352 end 497 end
353 498
354 end 499 end
355 500
501 function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond)
502
503 % Grid
504 xlims = {-3, 20};
505 ylims = {-11,5};
506 Lx = xlims{2} - xlims{1};
507 Ly = ylims{2} - ylims{1};
508
509 m = [15, 16];
510 g = grid.equidistant(m, xlims, ylims);
511 X = g.points();
512 x = g.x;
513
514 % Quadrature
515 opsx = sbp.D2Standard(m(1), xlims, order);
516 opsy = sbp.D2Standard(m(2), ylims, order);
517 Hx = opsx.H;
518 Hy = opsy.H;
519 H = {Hx, Hy};
520
521 % Moment conditions
522 fs = cell(mom_cond,1);
523 for p = 0:mom_cond-1
524 fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
525 end
526
527 % Grid spacing in interior
528 mm = round(m/2);
529 hx = x{1}(mm(1)+1) - x{1}(mm(1));
530 hy = x{2}(mm(2)+1) - x{2}(mm(2));
531 h = {hx, hy};
532
533 end
534
535 function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond)
536
537 % Grid
538 xlims = {-3, 20};
539 ylims = {-11,5};
540 zlims = {2,4};
541 Lx = xlims{2} - xlims{1};
542 Ly = ylims{2} - ylims{1};
543 Lz = zlims{2} - zlims{1};
544
545 m = [13, 14, 15];
546 g = grid.equidistant(m, xlims, ylims, zlims);
547 X = g.points();
548 x = g.x;
549
550 % Quadrature
551 opsx = sbp.D2Standard(m(1), xlims, order);
552 opsy = sbp.D2Standard(m(2), ylims, order);
553 opsz = sbp.D2Standard(m(3), zlims, order);
554 Hx = opsx.H;
555 Hy = opsy.H;
556 Hz = opsz.H;
557 H = {Hx, Hy, Hz};
558
559 % Moment conditions
560 fs = cell(mom_cond,1);
561 for p = 0:mom_cond-1
562 fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p;
563 end
564
565 % Grid spacing in interior
566 mm = round(m/2);
567 hx = x{1}(mm(1)+1) - x{1}(mm(1));
568 hy = x{2}(mm(2)+1) - x{2}(mm(2));
569 hz = x{3}(mm(3)+1) - x{3}(mm(3));
570 h = {hx, hy, hz};
571
572 end
573
356 function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond) 574 function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond)
357 575
358 % Grid 576 % Grid
359 xl = -3; 577 xl = -3;
360 xr = 900; 578 xr = 900;
373 for p = 0:mom_cond-1 591 for p = 0:mom_cond-1
374 fs{p+1} = @(x) (x/L).^p; 592 fs{p+1} = @(x) (x/L).^p;
375 end 593 end
376 594
377 end 595 end
378
379