Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 890:c70131daaa6e feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 21 Nov 2018 18:29:29 -0800
parents 5264ce57b573 43a1c3ac07b1
children
comparison
equal deleted inserted replaced
885:18e10217dca9 890:c70131daaa6e
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)
31 31
32 function testLeftRandom(testCase) 32 function testLeftRandom(testCase)
33 33
34 orders = [2, 4, 6]; 34 orders = [2, 4, 6];
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)
58 58
59 function testRightGP(testCase) 59 function testRightGP(testCase)
60 60
61 orders = [2, 4, 6]; 61 orders = [2, 4, 6];
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)
85 85
86 function testRightRandom(testCase) 86 function testRightRandom(testCase)
87 87
88 orders = [2, 4, 6]; 88 orders = [2, 4, 6];
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)
112 112
113 function testInteriorGP(testCase) 113 function testInteriorGP(testCase)
114 114
115 orders = [2, 4, 6]; 115 orders = [2, 4, 6];
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
140 140
141 function testInteriorRandom(testCase) 141 function testInteriorRandom(testCase)
142 142
143 orders = [2, 4, 6]; 143 orders = [2, 4, 6];
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)
168 % x0 outside grid should yield 0 integral! 168 % x0 outside grid should yield 0 integral!
169 function testX0OutsideGrid(testCase) 169 function testX0OutsideGrid(testCase)
170 170
171 orders = [2, 4, 6]; 171 orders = [2, 4, 6];
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)
195 195
196 function testAllGP(testCase) 196 function testAllGP(testCase)
197 197
198 orders = [2, 4, 6]; 198 orders = [2, 4, 6];
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)
222 222
223 function testHalfGP(testCase) 223 function testHalfGP(testCase)
224 224
225 orders = [2, 4, 6]; 225 orders = [2, 4, 6];
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)
245 end 245 end
246 end 246 end
247 end 247 end
248 end 248 end
249 249
250 function testAllGPStaggered(testCase) 250 % function testAllGPStaggered(testCase)
251 251
252 orders = [2, 4, 6]; 252 % orders = [2, 4, 6];
253 mom_conds = orders; 253 % mom_conds = orders;
254 254
255 for o = 1:length(orders) 255 % for o = 1:length(orders)
256 order = orders(o); 256 % order = orders(o);
257 mom_cond = mom_conds(o); 257 % mom_cond = mom_conds(o);
258 [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); 258 % [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
259
260 % % Test all grid points
261 % x0s = x;
262
263 % for j = 1:length(fs)
264 % f = fs{j};
265 % fx = f(x);
266 % for i = 1:length(x0s)
267 % x0 = x0s(i);
268 % delta = diracDiscr(x0, x, mom_cond, 0, H);
269 % integral = delta'*H*fx;
270 % err = abs(integral - f(x0));
271 % testCase.verifyLessThan(err, 1e-12);
272 % end
273 % end
274 % end
275 % end
276
277 % function testHalfGPStaggered(testCase)
278
279 % orders = [2, 4, 6];
280 % mom_conds = orders;
281
282 % for o = 1:length(orders)
283 % order = orders(o);
284 % mom_cond = mom_conds(o);
285 % [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
286
287 % % Test halfway between all grid points
288 % x0s = 1/2*( x(2:end)+x(1:end-1) );
289
290 % for j = 1:length(fs)
291 % f = fs{j};
292 % fx = f(x);
293 % for i = 1:length(x0s)
294 % x0 = x0s(i);
295 % delta = diracDiscr(x0, x, mom_cond, 0, H);
296 % integral = delta'*H*fx;
297 % err = abs(integral - f(x0));
298 % testCase.verifyLessThan(err, 1e-12);
299 % end
300 % end
301 % end
302 % end
303
304 % function testRandomStaggered(testCase)
305
306 % orders = [2, 4, 6];
307 % mom_conds = orders;
308
309 % for o = 1:length(orders)
310 % order = orders(o);
311 % mom_cond = mom_conds(o);
312 % [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond);
313
314 % % Test random points within grid boundaries
315 % x0s = xl + (xr-xl)*rand(1,300);
316
317 % for j = 1:length(fs)
318 % f = fs{j};
319 % fx = f(x);
320 % for i = 1:length(x0s)
321 % x0 = x0s(i);
322 % delta = diracDiscr(x0, x, mom_cond, 0, H);
323 % integral = delta'*H*fx;
324 % err = abs(integral - f(x0));
325 % testCase.verifyLessThan(err, 1e-12);
326 % end
327 % end
328 % end
329 % end
330
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});
259 342
260 % Test all grid points 343 % Test all grid points
261 x0s = x; 344 x0s = X;
262 345
263 for j = 1:length(fs) 346 for j = 1:length(fs)
264 f = fs{j}; 347 f = fs{j};
265 fx = f(x); 348 fx = f(X(:,1), X(:,2));
266 for i = 1:length(x0s) 349 for i = 1:length(x0s)
267 x0 = x0s(i); 350 x0 = x0s(i,:);
268 delta = diracDiscr(x0, x, mom_cond, 0, H); 351 delta = diracDiscr(x0, x, mom_cond, 0, H);
269 integral = delta'*H*fx; 352 integral = delta'*H_global*fx;
270 err = abs(integral - f(x0)); 353 err = abs(integral - f(x0(1), x0(2)));
271 testCase.verifyLessThan(err, 1e-12); 354 testCase.verifyLessThan(err, 1e-12);
272 end 355 end
273 end 356 end
274 end 357 end
275 end 358 end
276 359
277 function testHalfGPStaggered(testCase) 360 function testAllRandom2D(testCase)
278 361
279 orders = [2, 4, 6]; 362 orders = [2, 4, 6];
280 mom_conds = orders; 363 mom_conds = orders;
281 364
282 for o = 1:length(orders) 365 for o = 1:length(orders)
283 order = orders(o); 366 order = orders(o);
284 mom_cond = mom_conds(o); 367 mom_cond = mom_conds(o);
285 [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); 368 [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond);
286 369 H_global = kron(H{1}, H{2});
287 % Test halfway between all grid points 370
288 x0s = 1/2*( x(2:end)+x(1:end-1) ); 371 xl = xlims{1};
289 372 xr = xlims{2};
290 for j = 1:length(fs) 373 yl = ylims{1};
291 f = fs{j}; 374 yr = ylims{2};
292 fx = f(x); 375
293 for i = 1:length(x0s) 376 % Test random points, even outside grid
294 x0 = x0s(i); 377 Npoints = 100;
295 delta = diracDiscr(x0, x, mom_cond, 0, H); 378 x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ...
296 integral = delta'*H*fx; 379 (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ];
297 err = abs(integral - f(x0)); 380
298 testCase.verifyLessThan(err, 1e-12); 381 for j = 1:length(fs)
299 end 382 f = fs{j};
300 end 383 fx = f(X(:,1), X(:,2));
301 end 384 for i = 1:length(x0s)
302 end 385 x0 = x0s(i,:);
303 386 delta = diracDiscr(x0, x, mom_cond, 0, H);
304 function testRandomStaggered(testCase) 387 integral = delta'*H_global*fx;
305 388
306 orders = [2, 4, 6]; 389 % Integral should be 0 if point is outside grid
307 mom_conds = orders; 390 if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr
308 391 err = abs(integral - 0);
309 for o = 1:length(orders) 392 else
310 order = orders(o); 393 err = abs(integral - f(x0(1), x0(2)));
311 mom_cond = mom_conds(o); 394 end
312 [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); 395 testCase.verifyLessThan(err, 1e-12);
313 396 end
314 % Test random points within grid boundaries 397 end
315 x0s = xl + (xr-xl)*rand(1,300); 398 end
316 399 end
317 for j = 1:length(fs) 400
318 f = fs{j}; 401 %=============== 3D tests ==============================
319 fx = f(x); 402 function testAllGP3D(testCase)
320 for i = 1:length(x0s) 403
321 x0 = x0s(i); 404 orders = [2, 4, 6];
322 delta = diracDiscr(x0, x, mom_cond, 0, H); 405 mom_conds = orders;
323 integral = delta'*H*fx; 406
324 err = abs(integral - f(x0)); 407 for o = 1:length(orders)
325 testCase.verifyLessThan(err, 1e-12); 408 order = orders(o);
326 end 409 mom_cond = mom_conds(o);
327 end 410 [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
328 end 411 H_global = kron(kron(H{1}, H{2}), H{3});
329 end 412
330 413 % Test all grid points
331 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