Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 1128:3a9262c045d0 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Copy diracDiscr.m from feature/poroelastic
author Martin Almquist <malmquist@stanford.edu>
date Tue, 21 May 2019 17:59:30 -0700
parents
children 52d774e69b1f
comparison
equal deleted inserted replaced
1127:0aed89043ad6 1128:3a9262c045d0
1 function tests = diracDiscrTest()
2 tests = functiontests(localfunctions);
3 end
4
5 function testLeftGP(testCase)
6
7 orders = [2, 4, 6];
8 mom_conds = orders;
9
10 for o = 1:length(orders)
11 order = orders(o);
12 mom_cond = mom_conds(o);
13 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
14
15 % Test left boundary grid points
16 x0s = xl + [0, h, 2*h];
17
18 for j = 1:length(fs)
19 f = fs{j};
20 fx = f(x);
21 for i = 1:length(x0s)
22 x0 = x0s(i);
23 delta = diracDiscr(x0, x, mom_cond, 0, H);
24 integral = delta'*H*fx;
25 err = abs(integral - f(x0));
26 testCase.verifyLessThan(err, 1e-12);
27 end
28 end
29 end
30 end
31
32 function testLeftRandom(testCase)
33
34 orders = [2, 4, 6];
35 mom_conds = orders;
36
37 for o = 1:length(orders)
38 order = orders(o);
39 mom_cond = mom_conds(o);
40 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
41
42 % Test random points near left boundary
43 x0s = xl + 2*h*rand(1,10);
44
45 for j = 1:length(fs)
46 f = fs{j};
47 fx = f(x);
48 for i = 1:length(x0s)
49 x0 = x0s(i);
50 delta = diracDiscr(x0, x, mom_cond, 0, H);
51 integral = delta'*H*fx;
52 err = abs(integral - f(x0));
53 testCase.verifyLessThan(err, 1e-12);
54 end
55 end
56 end
57 end
58
59 function testRightGP(testCase)
60
61 orders = [2, 4, 6];
62 mom_conds = orders;
63
64 for o = 1:length(orders)
65 order = orders(o);
66 mom_cond = mom_conds(o);
67 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
68
69 % Test right boundary grid points
70 x0s = xr-[0, h, 2*h];
71
72 for j = 1:length(fs)
73 f = fs{j};
74 fx = f(x);
75 for i = 1:length(x0s)
76 x0 = x0s(i);
77 delta = diracDiscr(x0, x, mom_cond, 0, H);
78 integral = delta'*H*fx;
79 err = abs(integral - f(x0));
80 testCase.verifyLessThan(err, 1e-12);
81 end
82 end
83 end
84 end
85
86 function testRightRandom(testCase)
87
88 orders = [2, 4, 6];
89 mom_conds = orders;
90
91 for o = 1:length(orders)
92 order = orders(o);
93 mom_cond = mom_conds(o);
94 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
95
96 % Test random points near right boundary
97 x0s = xr - 2*h*rand(1,10);
98
99 for j = 1:length(fs)
100 f = fs{j};
101 fx = f(x);
102 for i = 1:length(x0s)
103 x0 = x0s(i);
104 delta = diracDiscr(x0, x, mom_cond, 0, H);
105 integral = delta'*H*fx;
106 err = abs(integral - f(x0));
107 testCase.verifyLessThan(err, 1e-12);
108 end
109 end
110 end
111 end
112
113 function testInteriorGP(testCase)
114
115 orders = [2, 4, 6];
116 mom_conds = orders;
117
118 for o = 1:length(orders)
119 order = orders(o);
120 mom_cond = mom_conds(o);
121 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
122
123 % Test interior grid points
124 m_half = round(m/2);
125 x0s = xl + (m_half-1:m_half+1)*h;
126
127 for j = 1:length(fs)
128 f = fs{j};
129 fx = f(x);
130 for i = 1:length(x0s)
131 x0 = x0s(i);
132 delta = diracDiscr(x0, x, mom_cond, 0, H);
133 integral = delta'*H*fx;
134 err = abs(integral - f(x0));
135 testCase.verifyLessThan(err, 1e-12);
136 end
137 end
138 end
139 end
140
141 function testInteriorRandom(testCase)
142
143 orders = [2, 4, 6];
144 mom_conds = orders;
145
146 for o = 1:length(orders)
147 order = orders(o);
148 mom_cond = mom_conds(o);
149 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
150
151 % Test random points in interior
152 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
153
154 for j = 1:length(fs)
155 f = fs{j};
156 fx = f(x);
157 for i = 1:length(x0s)
158 x0 = x0s(i);
159 delta = diracDiscr(x0, x, mom_cond, 0, H);
160 integral = delta'*H*fx;
161 err = abs(integral - f(x0));
162 testCase.verifyLessThan(err, 1e-12);
163 end
164 end
165 end
166 end
167
168 % x0 outside grid should yield 0 integral!
169 function testX0OutsideGrid(testCase)
170
171 orders = [2, 4, 6];
172 mom_conds = orders;
173
174 for o = 1:length(orders)
175 order = orders(o);
176 mom_cond = mom_conds(o);
177 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
178
179 % Test points outisde grid
180 x0s = [xl-1.1*h, xr+1.1*h];
181
182 for j = 1:length(fs)
183 f = fs{j};
184 fx = f(x);
185 for i = 1:length(x0s)
186 x0 = x0s(i);
187 delta = diracDiscr(x0, x, mom_cond, 0, H);
188 integral = delta'*H*fx;
189 err = abs(integral - 0);
190 testCase.verifyLessThan(err, 1e-12);
191 end
192 end
193 end
194 end
195
196 function testAllGP(testCase)
197
198 orders = [2, 4, 6];
199 mom_conds = orders;
200
201 for o = 1:length(orders)
202 order = orders(o);
203 mom_cond = mom_conds(o);
204 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
205
206 % Test all grid points
207 x0s = x;
208
209 for j = 1:length(fs)
210 f = fs{j};
211 fx = f(x);
212 for i = 1:length(x0s)
213 x0 = x0s(i);
214 delta = diracDiscr(x0, x, mom_cond, 0, H);
215 integral = delta'*H*fx;
216 err = abs(integral - f(x0));
217 testCase.verifyLessThan(err, 1e-12);
218 end
219 end
220 end
221 end
222
223 function testHalfGP(testCase)
224
225 orders = [2, 4, 6];
226 mom_conds = orders;
227
228 for o = 1:length(orders)
229 order = orders(o);
230 mom_cond = mom_conds(o);
231 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond);
232
233 % Test halfway between all grid points
234 x0s = 1/2*( x(2:end)+x(1:end-1) );
235
236 for j = 1:length(fs)
237 f = fs{j};
238 fx = f(x);
239 for i = 1:length(x0s)
240 x0 = x0s(i);
241 delta = diracDiscr(x0, x, mom_cond, 0, H);
242 integral = delta'*H*fx;
243 err = abs(integral - f(x0));
244 testCase.verifyLessThan(err, 1e-12);
245 end
246 end
247 end
248 end
249
250 % function testAllGPStaggered(testCase)
251
252 % orders = [2, 4, 6];
253 % mom_conds = orders;
254
255 % for o = 1:length(orders)
256 % order = orders(o);
257 % mom_cond = mom_conds(o);
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});
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 % ======================================================
476 % ============== Setup functions =======================
477 % ======================================================
478 function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond)
479
480 % Grid
481 xl = -3;
482 xr = 900;
483 L = xr-xl;
484 m = 101;
485 h = (xr-xl)/(m-1);
486 g = grid.equidistant(m, {xl, xr});
487 x = g.points();
488
489 % Quadrature
490 ops = sbp.D2Standard(m, {xl, xr}, order);
491 H = ops.H;
492
493 % Moment conditions
494 fs = cell(mom_cond,1);
495 for p = 0:mom_cond-1
496 fs{p+1} = @(x) (x/L).^p;
497 end
498
499 end
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
574 function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond)
575
576 % Grid
577 xl = -3;
578 xr = 900;
579 L = xr-xl;
580 m = 101;
581 [~, g_dual] = grid.primalDual1D(m, {xl, xr});
582 x = g_dual.points();
583 h = g_dual.h;
584
585 % Quadrature
586 ops = sbp.D1Staggered(m, {xl, xr}, order);
587 H = ops.H_dual;
588
589 % Moment conditions
590 fs = cell(mom_cond,1);
591 for p = 0:mom_cond-1
592 fs{p+1} = @(x) (x/L).^p;
593 end
594
595 end