Mercurial > repos > public > sbplib

comparison diracPrimDiscrTest.m @ 1130:99fd66ffe714 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add derivative of delta functions and corresponding tests, tested for 1D.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 21 May 2019 18:44:01 -0700
parents
children ea225a4659fe
comparison
equal deleted inserted replaced
1129:b29892853daf 1130:99fd66ffe714
1 function tests = diracPrimDiscrTest()
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, fps] = 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 fp = fps{j};
21 fx = f(x);
22 for i = 1:length(x0s)
23 x0 = x0s(i);
24 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
25 integral = delta'*H*fx;
26 err = abs(integral + fp(x0));
27 testCase.verifyLessThan(err, 1e-12);
28 end
29 end
30 end
31 end
32
33 function testLeftRandom(testCase)
34
35 orders = [2, 4, 6];
36 mom_conds = orders;
37
38 for o = 1:length(orders)
39 order = orders(o);
40 mom_cond = mom_conds(o);
41 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
42
43 % Test random points near left boundary
44 x0s = xl + 2*h*rand(1,10);
45
46 for j = 1:length(fs)
47 f = fs{j};
48 fp = fps{j};
49 fx = f(x);
50 for i = 1:length(x0s)
51 x0 = x0s(i);
52 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
53 integral = delta'*H*fx;
54 err = abs(integral + fp(x0));
55 testCase.verifyLessThan(err, 1e-12);
56 end
57 end
58 end
59 end
60
61 function testRightGP(testCase)
62
63 orders = [2, 4, 6];
64 mom_conds = orders;
65
66 for o = 1:length(orders)
67 order = orders(o);
68 mom_cond = mom_conds(o);
69 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
70
71 % Test right boundary grid points
72 x0s = xr-[0, h, 2*h];
73
74 for j = 1:length(fs)
75 f = fs{j};
76 fp = fps{j};
77 fx = f(x);
78 for i = 1:length(x0s)
79 x0 = x0s(i);
80 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
81 integral = delta'*H*fx;
82 err = abs(integral + fp(x0));
83 testCase.verifyLessThan(err, 1e-12);
84 end
85 end
86 end
87 end
88
89 function testRightRandom(testCase)
90
91 orders = [2, 4, 6];
92 mom_conds = orders;
93
94 for o = 1:length(orders)
95 order = orders(o);
96 mom_cond = mom_conds(o);
97 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
98
99 % Test random points near right boundary
100 x0s = xr - 2*h*rand(1,10);
101
102 for j = 1:length(fs)
103 f = fs{j};
104 fp = fps{j};
105 fx = f(x);
106 for i = 1:length(x0s)
107 x0 = x0s(i);
108 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
109 integral = delta'*H*fx;
110 err = abs(integral + fp(x0));
111 testCase.verifyLessThan(err, 1e-12);
112 end
113 end
114 end
115 end
116
117 function testInteriorGP(testCase)
118
119 orders = [2, 4, 6];
120 mom_conds = orders;
121
122 for o = 1:length(orders)
123 order = orders(o);
124 mom_cond = mom_conds(o);
125 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
126
127 % Test interior grid points
128 m_half = round(m/2);
129 x0s = xl + (m_half-1:m_half+1)*h;
130
131 for j = 1:length(fs)
132 f = fs{j};
133 fp = fps{j};
134 fx = f(x);
135 for i = 1:length(x0s)
136 x0 = x0s(i);
137 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
138 integral = delta'*H*fx;
139 err = abs(integral + fp(x0));
140 testCase.verifyLessThan(err, 1e-12);
141 end
142 end
143 end
144 end
145
146 function testInteriorRandom(testCase)
147
148 orders = [2, 4, 6];
149 mom_conds = orders;
150
151 for o = 1:length(orders)
152 order = orders(o);
153 mom_cond = mom_conds(o);
154 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
155
156 % Test random points in interior
157 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
158
159 for j = 1:length(fs)
160 f = fs{j};
161 fp = fps{j};
162 fx = f(x);
163 for i = 1:length(x0s)
164 x0 = x0s(i);
165 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
166 integral = delta'*H*fx;
167 err = abs(integral + fp(x0));
168 testCase.verifyLessThan(err, 1e-12);
169 end
170 end
171 end
172 end
173
174 % x0 outside grid should yield 0 integral!
175 function testX0OutsideGrid(testCase)
176
177 orders = [2, 4, 6];
178 mom_conds = orders;
179
180 for o = 1:length(orders)
181 order = orders(o);
182 mom_cond = mom_conds(o);
183 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
184
185 % Test points outisde grid
186 x0s = [xl-1.1*h, xr+1.1*h];
187
188 for j = 1:length(fs)
189 f = fs{j};
190 fp = fps{j};
191 fx = f(x);
192 for i = 1:length(x0s)
193 x0 = x0s(i);
194 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
195 integral = delta'*H*fx;
196 err = abs(integral - 0);
197 testCase.verifyLessThan(err, 1e-12);
198 end
199 end
200 end
201 end
202
203 function testAllGP(testCase)
204
205 orders = [2, 4, 6];
206 mom_conds = orders;
207
208 for o = 1:length(orders)
209 order = orders(o);
210 mom_cond = mom_conds(o);
211 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
212
213 % Test all grid points
214 x0s = x;
215
216 for j = 1:length(fs)
217 f = fs{j};
218 fp = fps{j};
219 fx = f(x);
220 for i = 1:length(x0s)
221 x0 = x0s(i);
222 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
223 integral = delta'*H*fx;
224 err = abs(integral + fp(x0));
225 testCase.verifyLessThan(err, 1e-12);
226 end
227 end
228 end
229 end
230
231 function testHalfGP(testCase)
232
233 orders = [2, 4, 6];
234 mom_conds = orders;
235
236 for o = 1:length(orders)
237 order = orders(o);
238 mom_cond = mom_conds(o);
239 [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond);
240
241 % Test halfway between all grid points
242 x0s = 1/2*( x(2:end)+x(1:end-1) );
243
244 for j = 1:length(fs)
245 f = fs{j};
246 fp = fps{j};
247 fx = f(x);
248 for i = 1:length(x0s)
249 x0 = x0s(i);
250 delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
251 integral = delta'*H*fx;
252 err = abs(integral + fp(x0));
253 testCase.verifyLessThan(err, 1e-12);
254 end
255 end
256 end
257 end
258
259
260 %=============== 2D tests ==============================
261 % function testAllGP2D(testCase)
262
263 % orders = [2, 4, 6];
264 % mom_conds = orders;
265
266 % for o = 1:length(orders)
267 % order = orders(o);
268 % mom_cond = mom_conds(o);
269 % [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond);
270 % H_global = kron(H{1}, H{2});
271
272 % % Test all grid points
273 % x0s = X;
274
275 % for j = 1:length(fs)
276 % f = fs{j};
277 % fx = f(X(:,1), X(:,2));
278 % for i = 1:length(x0s)
279 % x0 = x0s(i,:);
280 % delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
281 % integral = delta'*H_global*fx;
282 % err = abs(integral - f(x0(1), x0(2)));
283 % testCase.verifyLessThan(err, 1e-12);
284 % end
285 % end
286 % end
287 % end
288
289 % function testAllRandom2D(testCase)
290
291 % orders = [2, 4, 6];
292 % mom_conds = orders;
293
294 % for o = 1:length(orders)
295 % order = orders(o);
296 % mom_cond = mom_conds(o);
297 % [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond);
298 % H_global = kron(H{1}, H{2});
299
300 % xl = xlims{1};
301 % xr = xlims{2};
302 % yl = ylims{1};
303 % yr = ylims{2};
304
305 % % Test random points, even outside grid
306 % Npoints = 100;
307 % x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ...
308 % (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ];
309
310 % for j = 1:length(fs)
311 % f = fs{j};
312 % fx = f(X(:,1), X(:,2));
313 % for i = 1:length(x0s)
314 % x0 = x0s(i,:);
315 % delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
316 % integral = delta'*H_global*fx;
317
318 % % Integral should be 0 if point is outside grid
319 % if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr
320 % err = abs(integral - 0);
321 % else
322 % err = abs(integral - f(x0(1), x0(2)));
323 % end
324 % testCase.verifyLessThan(err, 1e-12);
325 % end
326 % end
327 % end
328 % end
329
330 %=============== 3D tests ==============================
331 % function testAllGP3D(testCase)
332
333 % orders = [2, 4, 6];
334 % mom_conds = orders;
335
336 % for o = 1:length(orders)
337 % order = orders(o);
338 % mom_cond = mom_conds(o);
339 % [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
340 % H_global = kron(kron(H{1}, H{2}), H{3});
341
342 % % Test all grid points
343 % x0s = X;
344
345 % for j = 1:length(fs)
346 % f = fs{j};
347 % fx = f(X(:,1), X(:,2), X(:,3));
348 % for i = 1:length(x0s)
349 % x0 = x0s(i,:);
350 % delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
351 % integral = delta'*H_global*fx;
352 % err = abs(integral - f(x0(1), x0(2), x0(3)));
353 % testCase.verifyLessThan(err, 1e-12);
354 % end
355 % end
356 % end
357 % end
358
359 % function testAllRandom3D(testCase)
360
361 % orders = [2, 4, 6];
362 % mom_conds = orders;
363
364 % for o = 1:length(orders)
365 % order = orders(o);
366 % mom_cond = mom_conds(o);
367 % [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond);
368 % H_global = kron(kron(H{1}, H{2}), H{3});
369
370 % xl = xlims{1};
371 % xr = xlims{2};
372 % yl = ylims{1};
373 % yr = ylims{2};
374 % zl = zlims{1};
375 % zr = zlims{2};
376
377 % % Test random points, even outside grid
378 % Npoints = 200;
379 % x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ...
380 % (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1), ...
381 % (zl-3*h{3}) + (zr-zl+6*h{3})*rand(Npoints,1) ];
382
383 % for j = 1:length(fs)
384 % f = fs{j};
385 % fx = f(X(:,1), X(:,2), X(:,3));
386 % for i = 1:length(x0s)
387 % x0 = x0s(i,:);
388 % delta = diracPrimDiscr(x0, x, mom_cond, 0, H);
389 % integral = delta'*H_global*fx;
390
391 % % Integral should be 0 if point is outside grid
392 % if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr
393 % err = abs(integral - 0);
394 % else
395 % err = abs(integral - f(x0(1), x0(2), x0(3)));
396 % end
397 % testCase.verifyLessThan(err, 1e-12);
398 % end
399 % end
400 % end
401 % end
402
403
404 % ======================================================
405 % ============== Setup functions =======================
406 % ======================================================
407 function [xl, xr, m, h, x, H, fs, fps] = setup1D(order, mom_cond)
408
409 % Grid
410 xl = -3;
411 xr = 900;
412 L = xr-xl;
413 m = 101;
414 h = (xr-xl)/(m-1);
415 g = grid.equidistant(m, {xl, xr});
416 x = g.points();
417
418 % Quadrature
419 ops = sbp.D2Standard(m, {xl, xr}, order);
420 H = ops.H;
421
422 % Moment conditions
423 fs = cell(mom_cond+1,1);
424 fps = cell(mom_cond+1,1);
425
426 for p = 0:mom_cond
427 fs{p+1} = @(x) (x/L).^p;
428 end
429
430 fps{1} = @(x) 0*x;
431 for p = 1:mom_cond
432 fps{p+1} = @(x) p/L*(x/L).^(p-1);
433 end
434
435 end
436
437 % function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond)
438
439 % % Grid
440 % xlims = {-3, 20};
441 % ylims = {-11,5};
442 % Lx = xlims{2} - xlims{1};
443 % Ly = ylims{2} - ylims{1};
444
445 % m = [15, 16];
446 % g = grid.equidistant(m, xlims, ylims);
447 % X = g.points();
448 % x = g.x;
449
450 % % Quadrature
451 % opsx = sbp.D2Standard(m(1), xlims, order);
452 % opsy = sbp.D2Standard(m(2), ylims, order);
453 % Hx = opsx.H;
454 % Hy = opsy.H;
455 % H = {Hx, Hy};
456
457 % % Moment conditions
458 % fs = cell(mom_cond,1);
459 % for p = 0:mom_cond-1
460 % fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
461 % end
462
463 % % Grid spacing in interior
464 % mm = round(m/2);
465 % hx = x{1}(mm(1)+1) - x{1}(mm(1));
466 % hy = x{2}(mm(2)+1) - x{2}(mm(2));
467 % h = {hx, hy};
468
469 % end
470
471 % function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond)
472
473 % % Grid
474 % xlims = {-3, 20};
475 % ylims = {-11,5};
476 % zlims = {2,4};
477 % Lx = xlims{2} - xlims{1};
478 % Ly = ylims{2} - ylims{1};
479 % Lz = zlims{2} - zlims{1};
480
481 % m = [13, 14, 15];
482 % g = grid.equidistant(m, xlims, ylims, zlims);
483 % X = g.points();
484 % x = g.x;
485
486 % % Quadrature
487 % opsx = sbp.D2Standard(m(1), xlims, order);
488 % opsy = sbp.D2Standard(m(2), ylims, order);
489 % opsz = sbp.D2Standard(m(3), zlims, order);
490 % Hx = opsx.H;
491 % Hy = opsy.H;
492 % Hz = opsz.H;
493 % H = {Hx, Hy, Hz};
494
495 % % Moment conditions
496 % fs = cell(mom_cond,1);
497 % for p = 0:mom_cond-1
498 % fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p;
499 % end
500
501 % % Grid spacing in interior
502 % mm = round(m/2);
503 % hx = x{1}(mm(1)+1) - x{1}(mm(1));
504 % hy = x{2}(mm(2)+1) - x{2}(mm(2));
505 % hz = x{3}(mm(3)+1) - x{3}(mm(3));
506 % h = {hx, hy, hz};
507
508 % end