Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 1232:52d774e69b1f feature/dirac_discr

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up diracDiscr, remove obsolete tests.
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Tue, 19 Nov 2019 13:54:41 -0800
parents 86ee5648e384
children f1806475498b
comparison
equal deleted inserted replaced
1229:86ee5648e384 1232:52d774e69b1f
1 function tests = diracDiscrTest() 1 function tests = diracDiscrTest()
2 tests = functiontests(localfunctions); 2 tests = functiontests(localfunctions);
3 end 3 end
4 4
5 function testLeftGP(testCase) 5 %TODO: Test discretizing with smoothness conditions.
6 6 % Only discretization with moment conditions currently tested.
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 7
32 function testLeftRandom(testCase) 8 function testLeftRandom(testCase)
33 9
34 orders = [2, 4, 6]; 10 orders = [2, 4, 6];
35 mom_conds = orders; 11 mom_conds = orders;
36 12 rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
37 for o = 1:length(orders) 13
38 order = orders(o); 14 for o = 1:length(orders)
39 mom_cond = mom_conds(o); 15 order = orders(o);
40 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); 16 mom_cond = mom_conds(o);
17 [xl, ~, h, x, H, fs] = setup1D(order, mom_cond);
41 18
42 % Test random points near left boundary 19 % Test random points near left boundary
43 x0s = xl + 2*h*rand(1,10); 20 x0s = xl + 2*h*rand(1,10);
44 21
45 for j = 1:length(fs) 22 for j = 1:length(fs)
54 end 31 end
55 end 32 end
56 end 33 end
57 end 34 end
58 35
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) 36 function testRightRandom(testCase)
87 37
88 orders = [2, 4, 6]; 38 orders = [2, 4, 6];
89 mom_conds = orders; 39 mom_conds = orders;
90 40 rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
91 for o = 1:length(orders) 41
92 order = orders(o); 42 for o = 1:length(orders)
93 mom_cond = mom_conds(o); 43 order = orders(o);
94 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); 44 mom_cond = mom_conds(o);
45 [~, xr, h, x, H, fs] = setup1D(order, mom_cond);
95 46
96 % Test random points near right boundary 47 % Test random points near right boundary
97 x0s = xr - 2*h*rand(1,10); 48 x0s = xr - 2*h*rand(1,10);
98 49
99 for j = 1:length(fs) 50 for j = 1:length(fs)
108 end 59 end
109 end 60 end
110 end 61 end
111 end 62 end
112 63
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) 64 function testInteriorRandom(testCase)
142 65
143 orders = [2, 4, 6]; 66 orders = [2, 4, 6];
144 mom_conds = orders; 67 mom_conds = orders;
145 68 rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
146 for o = 1:length(orders) 69
147 order = orders(o); 70 for o = 1:length(orders)
148 mom_cond = mom_conds(o); 71 order = orders(o);
149 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); 72 mom_cond = mom_conds(o);
73 [xl, xr, h, x, H, fs] = setup1D(order, mom_cond);
150 74
151 % Test random points in interior 75 % Test random points in interior
152 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); 76 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
153 77
154 for j = 1:length(fs) 78 for j = 1:length(fs)
172 mom_conds = orders; 96 mom_conds = orders;
173 97
174 for o = 1:length(orders) 98 for o = 1:length(orders)
175 order = orders(o); 99 order = orders(o);
176 mom_cond = mom_conds(o); 100 mom_cond = mom_conds(o);
177 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); 101 [xl, xr, h, x, H, fs] = setup1D(order, mom_cond);
178 102
179 % Test points outisde grid 103 % Test points outisde grid
180 x0s = [xl-1.1*h, xr+1.1*h]; 104 x0s = [xl-1.1*h, xr+1.1*h];
181 105
182 for j = 1:length(fs) 106 for j = 1:length(fs)
199 mom_conds = orders; 123 mom_conds = orders;
200 124
201 for o = 1:length(orders) 125 for o = 1:length(orders)
202 order = orders(o); 126 order = orders(o);
203 mom_cond = mom_conds(o); 127 mom_cond = mom_conds(o);
204 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); 128 [~, ~, ~, x, H, fs] = setup1D(order, mom_cond);
205 129
206 % Test all grid points 130 % Test all grid points
207 x0s = x; 131 x0s = x;
208 132
209 for j = 1:length(fs) 133 for j = 1:length(fs)
226 mom_conds = orders; 150 mom_conds = orders;
227 151
228 for o = 1:length(orders) 152 for o = 1:length(orders)
229 order = orders(o); 153 order = orders(o);
230 mom_cond = mom_conds(o); 154 mom_cond = mom_conds(o);
231 [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); 155 [~, ~, ~, x, H, fs] = setup1D(order, mom_cond);
232 156
233 % Test halfway between all grid points 157 % Test halfway between all grid points
234 x0s = 1/2*( x(2:end)+x(1:end-1) ); 158 x0s = 1/2*( x(2:end)+x(1:end-1) );
235 159
236 for j = 1:length(fs) 160 for j = 1:length(fs)
245 end 169 end
246 end 170 end
247 end 171 end
248 end 172 end
249 173
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 ============================== 174 %=============== 2D tests ==============================
332 function testAllGP2D(testCase) 175 function testAllGP2D(testCase)
333 176
334 orders = [2, 4, 6]; 177 orders = [2, 4, 6];
335 mom_conds = orders; 178 mom_conds = orders;
336 179
337 for o = 1:length(orders) 180 for o = 1:length(orders)
338 order = orders(o); 181 order = orders(o);
339 mom_cond = mom_conds(o); 182 mom_cond = mom_conds(o);
340 [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond); 183 [~, ~, x, X, ~, H, fs] = setup2D(order, mom_cond);
341 H_global = kron(H{1}, H{2}); 184 H_global = kron(H{1}, H{2});
342 185
343 % Test all grid points 186 % Test all grid points
344 x0s = X; 187 x0s = X;
345 188
359 202
360 function testAllRandom2D(testCase) 203 function testAllRandom2D(testCase)
361 204
362 orders = [2, 4, 6]; 205 orders = [2, 4, 6];
363 mom_conds = orders; 206 mom_conds = orders;
364 207 rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
365 for o = 1:length(orders) 208
366 order = orders(o); 209 for o = 1:length(orders)
367 mom_cond = mom_conds(o); 210 order = orders(o);
368 [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond); 211 mom_cond = mom_conds(o);
212 [xlims, ylims, x, X, h, H, fs] = setup2D(order, mom_cond);
369 H_global = kron(H{1}, H{2}); 213 H_global = kron(H{1}, H{2});
370 214
371 xl = xlims{1}; 215 xl = xlims{1};
372 xr = xlims{2}; 216 xr = xlims{2};
373 yl = ylims{1}; 217 yl = ylims{1};
405 mom_conds = orders; 249 mom_conds = orders;
406 250
407 for o = 1:length(orders) 251 for o = 1:length(orders)
408 order = orders(o); 252 order = orders(o);
409 mom_cond = mom_conds(o); 253 mom_cond = mom_conds(o);
410 [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); 254 [~, ~, ~, x, X, ~, H, fs] = setup3D(order, mom_cond);
411 H_global = kron(kron(H{1}, H{2}), H{3}); 255 H_global = kron(kron(H{1}, H{2}), H{3});
412 256
413 % Test all grid points 257 % Test all grid points
414 x0s = X; 258 x0s = X;
415 259
429 273
430 function testAllRandom3D(testCase) 274 function testAllRandom3D(testCase)
431 275
432 orders = [2, 4, 6]; 276 orders = [2, 4, 6];
433 mom_conds = orders; 277 mom_conds = orders;
434 278 rng(0xDABBAD00) % Set seed. Jabba-dabba-doooo!
435 for o = 1:length(orders) 279
436 order = orders(o); 280 for o = 1:length(orders)
437 mom_cond = mom_conds(o); 281 order = orders(o);
438 [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); 282 mom_cond = mom_conds(o);
283 [xlims, ylims, zlims, x, X, h, H, fs] = setup3D(order, mom_cond);
439 H_global = kron(kron(H{1}, H{2}), H{3}); 284 H_global = kron(kron(H{1}, H{2}), H{3});
440 285
441 xl = xlims{1}; 286 xl = xlims{1};
442 xr = xlims{2}; 287 xr = xlims{2};
443 yl = ylims{1}; 288 yl = ylims{1};
473 318
474 319
475 % ====================================================== 320 % ======================================================
476 % ============== Setup functions ======================= 321 % ============== Setup functions =======================
477 % ====================================================== 322 % ======================================================
478 function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond) 323 function [xl, xr, h, x, H, fs] = setup1D(order, mom_cond)
479 324
480 % Grid 325 % Grid
481 xl = -3; 326 xl = -3;
482 xr = 900; 327 xr = 900;
483 L = xr-xl; 328 L = xr-xl;
496 fs{p+1} = @(x) (x/L).^p; 341 fs{p+1} = @(x) (x/L).^p;
497 end 342 end
498 343
499 end 344 end
500 345
501 function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond) 346 function [xlims, ylims, x, X, h, H, fs] = setup2D(order, mom_cond)
502 347
503 % Grid 348 % Grid
504 xlims = {-3, 20}; 349 xlims = {-3, 20};
505 ylims = {-11,5}; 350 ylims = {-11,5};
506 Lx = xlims{2} - xlims{1}; 351 Lx = xlims{2} - xlims{1};
530 hy = x{2}(mm(2)+1) - x{2}(mm(2)); 375 hy = x{2}(mm(2)+1) - x{2}(mm(2));
531 h = {hx, hy}; 376 h = {hx, hy};
532 377
533 end 378 end
534 379
535 function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond) 380 function [xlims, ylims, zlims, x, X, h, H, fs] = setup3D(order, mom_cond)
536 381
537 % Grid 382 % Grid
538 xlims = {-3, 20}; 383 xlims = {-3, 20};
539 ylims = {-11,5}; 384 ylims = {-11,5};
540 zlims = {2,4}; 385 zlims = {2,4};
568 hy = x{2}(mm(2)+1) - x{2}(mm(2)); 413 hy = x{2}(mm(2)+1) - x{2}(mm(2));
569 hz = x{3}(mm(3)+1) - x{3}(mm(3)); 414 hz = x{3}(mm(3)+1) - x{3}(mm(3));
570 h = {hx, hy, hz}; 415 h = {hx, hy, hz};
571 416
572 end 417 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