Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 1251:6424745e1b58 feature/volcano

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge in latest changes from default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 20 Nov 2019 14:26:57 -0800
parents 25efceb0c392
children 60c875c18de3
comparison
equal deleted inserted replaced
1248:10881b234f77 1251:6424745e1b58
1 function tests = diracDiscrTest()
2 tests = functiontests(localfunctions);
3 end
4
5 %TODO:
6 % 1. Test discretizing with smoothness conditions.
7 % Only discretization with moment conditions currently tested.
8 % 2. Test using other types of grids. Only equidistant grids currently used.
9
10 function testLeftRandom(testCase)
11
12 orders = [2, 4, 6];
13 mom_conds = orders;
14 rng(1) % Set seed
15
16 for o = 1:length(orders)
17 order = orders(o);
18 mom_cond = mom_conds(o);
19 [g, H, fs] = setup1D(order, mom_cond);
20 xl = g.lim{1}{1};
21 h = g.scaling();
22 x = g.x{1};
23
24 % Test random points near left boundary
25 x0s = xl + 2*h*rand(1,10);
26
27 for j = 1:length(fs)
28 f = fs{j};
29 fx = f(x);
30 for i = 1:length(x0s)
31 x0 = x0s(i);
32 delta = diracDiscr(g, x0, mom_cond, 0, H);
33 integral = delta'*H*fx;
34 err = abs(integral - f(x0));
35 testCase.verifyLessThan(err, 1e-12);
36 end
37 end
38 end
39 end
40
41 function testRightRandom(testCase)
42
43 orders = [2, 4, 6];
44 mom_conds = orders;
45 rng(1) % Set seed
46
47 for o = 1:length(orders)
48 order = orders(o);
49 mom_cond = mom_conds(o);
50 [g, H, fs] = setup1D(order, mom_cond);
51 xr = g.lim{1}{2};
52 h = g.scaling();
53 x = g.x{1};
54
55 % Test random points near right boundary
56 x0s = xr - 2*h*rand(1,10);
57
58 for j = 1:length(fs)
59 f = fs{j};
60 fx = f(x);
61 for i = 1:length(x0s)
62 x0 = x0s(i);
63 delta = diracDiscr(g, x0, mom_cond, 0, H);
64 integral = delta'*H*fx;
65 err = abs(integral - f(x0));
66 testCase.verifyLessThan(err, 1e-12);
67 end
68 end
69 end
70 end
71
72 function testInteriorRandom(testCase)
73
74 orders = [2, 4, 6];
75 mom_conds = orders;
76 rng(1) % Set seed
77
78 for o = 1:length(orders)
79 order = orders(o);
80 mom_cond = mom_conds(o);
81 [g, H, fs] = setup1D(order, mom_cond);
82 xl = g.lim{1}{1};
83 xr = g.lim{1}{2};
84 h = g.scaling();
85 x = g.x{1};
86
87 % Test random points in interior
88 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
89
90 for j = 1:length(fs)
91 f = fs{j};
92 fx = f(x);
93 for i = 1:length(x0s)
94 x0 = x0s(i);
95 delta = diracDiscr(g, x0, mom_cond, 0, H);
96 integral = delta'*H*fx;
97 err = abs(integral - f(x0));
98 testCase.verifyLessThan(err, 1e-12);
99 end
100 end
101 end
102 end
103
104 % x0 outside grid should yield 0 integral!
105 function testX0OutsideGrid(testCase)
106
107 orders = [2, 4, 6];
108 mom_conds = orders;
109
110 for o = 1:length(orders)
111 order = orders(o);
112 mom_cond = mom_conds(o);
113 [g, H, fs] = setup1D(order, mom_cond);
114 xl = g.lim{1}{1};
115 xr = g.lim{1}{2};
116 h = g.scaling();
117 x = g.x{1};
118
119 % Test points outisde grid
120 x0s = [xl-1.1*h, xr+1.1*h];
121
122 for j = 1:length(fs)
123 f = fs{j};
124 fx = f(x);
125 for i = 1:length(x0s)
126 x0 = x0s(i);
127 delta = diracDiscr(g, x0, mom_cond, 0, H);
128 integral = delta'*H*fx;
129 err = abs(integral - 0);
130 testCase.verifyLessThan(err, 1e-12);
131 end
132 end
133 end
134 end
135
136 function testAllGP(testCase)
137
138 orders = [2, 4, 6];
139 mom_conds = orders;
140
141 for o = 1:length(orders)
142 order = orders(o);
143 mom_cond = mom_conds(o);
144 [g, H, fs] = setup1D(order, mom_cond);
145 x = g.x{1};
146
147 % Test all grid points
148 x0s = x;
149
150 for j = 1:length(fs)
151 f = fs{j};
152 fx = f(x);
153 for i = 1:length(x0s)
154 x0 = x0s(i);
155 delta = diracDiscr(g, x0, mom_cond, 0, H);
156 integral = delta'*H*fx;
157 err = abs(integral - f(x0));
158 testCase.verifyLessThan(err, 1e-12);
159 end
160 end
161 end
162 end
163
164 function testHalfGP(testCase)
165
166 orders = [2, 4, 6];
167 mom_conds = orders;
168
169 for o = 1:length(orders)
170 order = orders(o);
171 mom_cond = mom_conds(o);
172 [g, H, fs] = setup1D(order, mom_cond);
173 x = g.x{1};
174
175 % Test halfway between all grid points
176 x0s = 1/2*(x(2:end)+x(1:end-1));
177
178 for j = 1:length(fs)
179 f = fs{j};
180 fx = f(x);
181 for i = 1:length(x0s)
182 x0 = x0s(i);
183 delta = diracDiscr(g, x0, mom_cond, 0, H);
184 integral = delta'*H*fx;
185 err = abs(integral - f(x0));
186 testCase.verifyLessThan(err, 1e-12);
187 end
188 end
189 end
190 end
191
192 %=============== 2D tests ==============================
193 function testAllGP2D(testCase)
194
195 orders = [2, 4, 6];
196 mom_conds = orders;
197
198 for o = 1:length(orders)
199 order = orders(o);
200 mom_cond = mom_conds(o);
201 [g, H, fs] = setup2D(order, mom_cond);
202 H_global = kron(H{1}, H{2});
203 X = g.points();
204 % Test all grid points
205 x0s = X;
206
207 for j = 1:length(fs)
208 f = fs{j};
209 fx = f(X(:,1), X(:,2));
210 for i = 1:length(x0s)
211 x0 = x0s(i,:);
212 delta = diracDiscr(g, x0, mom_cond, 0, H);
213 integral = delta'*H_global*fx;
214 err = abs(integral - f(x0(1), x0(2)));
215 testCase.verifyLessThan(err, 1e-12);
216 end
217 end
218 end
219 end
220
221 function testAllRandom2D(testCase)
222
223 orders = [2, 4, 6];
224 mom_conds = orders;
225 rng(1) % Set seed
226
227 for o = 1:length(orders)
228 order = orders(o);
229 mom_cond = mom_conds(o);
230 [g, H, fs] = setup2D(order, mom_cond);
231 H_global = kron(H{1}, H{2});
232 X = g.points();
233 xl = g.lim{1}{1};
234 xr = g.lim{1}{2};
235 yl = g.lim{2}{1};
236 yr = g.lim{2}{2};
237 h = g.scaling();
238
239
240 % Test random points, even outside grid
241 Npoints = 100;
242 x0s = [(xl-3*h(1)) + (xr-xl+6*h(1))*rand(Npoints,1), ...
243 (yl-3*h(2)) + (yr-yl+6*h(2))*rand(Npoints,1) ];
244
245 for j = 1:length(fs)
246 f = fs{j};
247 fx = f(X(:,1), X(:,2));
248 for i = 1:length(x0s)
249 x0 = x0s(i,:);
250 delta = diracDiscr(g, x0, mom_cond, 0, H);
251 integral = delta'*H_global*fx;
252
253 % Integral should be 0 if point is outside grid
254 if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr
255 err = abs(integral - 0);
256 else
257 err = abs(integral - f(x0(1), x0(2)));
258 end
259 testCase.verifyLessThan(err, 1e-12);
260 end
261 end
262 end
263 end
264
265 %=============== 3D tests ==============================
266 function testAllGP3D(testCase)
267
268 orders = [2, 4, 6];
269 mom_conds = orders;
270
271 for o = 1:length(orders)
272 order = orders(o);
273 mom_cond = mom_conds(o);
274 [g, H, fs] = setup3D(order, mom_cond);
275 H_global = kron(kron(H{1}, H{2}), H{3});
276 X = g.points();
277 % Test all grid points
278 x0s = X;
279
280 for j = 1:length(fs)
281 f = fs{j};
282 fx = f(X(:,1), X(:,2), X(:,3));
283 for i = 1:length(x0s)
284 x0 = x0s(i,:);
285 delta = diracDiscr(g, x0, mom_cond, 0, H);
286 integral = delta'*H_global*fx;
287 err = abs(integral - f(x0(1), x0(2), x0(3)));
288 testCase.verifyLessThan(err, 1e-12);
289 end
290 end
291 end
292 end
293
294 function testAllRandom3D(testCase)
295
296 orders = [2, 4, 6];
297 mom_conds = orders;
298 rng(1) % Set seed
299
300 for o = 1:length(orders)
301 order = orders(o);
302 mom_cond = mom_conds(o);
303 [g, H, fs] = setup3D(order, mom_cond);
304 H_global = kron(kron(H{1}, H{2}), H{3});
305 X = g.points();
306 xl = g.lim{1}{1};
307 xr = g.lim{1}{2};
308 yl = g.lim{2}{1};
309 yr = g.lim{2}{2};
310 zl = g.lim{3}{1};
311 zr = g.lim{3}{2};
312 h = g.scaling();
313
314 % Test random points, even outside grid
315 Npoints = 200;
316 x0s = [(xl-3*h(1)) + (xr-xl+6*h(1))*rand(Npoints,1), ...
317 (yl-3*h(2)) + (yr-yl+6*h(2))*rand(Npoints,1), ...
318 (zl-3*h(3)) + (zr-zl+6*h(3))*rand(Npoints,1) ];
319
320 for j = 1:length(fs)
321 f = fs{j};
322 fx = f(X(:,1), X(:,2), X(:,3));
323 for i = 1:length(x0s)
324 x0 = x0s(i,:);
325 delta = diracDiscr(g, x0, mom_cond, 0, H);
326 integral = delta'*H_global*fx;
327
328 % Integral should be 0 if point is outside grid
329 if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr
330 err = abs(integral - 0);
331 else
332 err = abs(integral - f(x0(1), x0(2), x0(3)));
333 end
334 testCase.verifyLessThan(err, 1e-12);
335 end
336 end
337 end
338 end
339
340
341 % ======================================================
342 % ============== Setup functions =======================
343 % ======================================================
344 function [g, H, fs] = setup1D(order, mom_cond)
345
346 % Grid
347 xl = -3;
348 xr = 900;
349 L = xr-xl;
350 m = 101;
351 g = grid.equidistant(m, {xl, xr});
352
353 % Quadrature
354 ops = sbp.D2Standard(m, {xl, xr}, order);
355 H = ops.H;
356
357 % Moment conditions
358 fs = cell(mom_cond,1);
359 for p = 0:mom_cond-1
360 fs{p+1} = @(x) (x/L).^p;
361 end
362
363 end
364
365 function [g, H, fs] = setup2D(order, mom_cond)
366
367 % Grid
368 xlims = {-3, 20};
369 ylims = {-11,5};
370 Lx = xlims{2} - xlims{1};
371 Ly = ylims{2} - ylims{1};
372 m = [15, 16];
373 g = grid.equidistant(m, xlims, ylims);
374
375 % Quadrature
376 opsx = sbp.D2Standard(m(1), xlims, order);
377 opsy = sbp.D2Standard(m(2), ylims, order);
378 Hx = opsx.H;
379 Hy = opsy.H;
380 H = {Hx, Hy};
381
382 % Moment conditions
383 fs = cell(mom_cond,1);
384 for p = 0:mom_cond-1
385 fs{p+1} = @(x,y) (x/Lx + y/Ly).^p;
386 end
387 end
388
389 function [g, H, fs] = setup3D(order, mom_cond)
390
391 % Grid
392 xlims = {-3, 20};
393 ylims = {-11,5};
394 zlims = {2,4};
395 Lx = xlims{2} - xlims{1};
396 Ly = ylims{2} - ylims{1};
397 Lz = zlims{2} - zlims{1};
398 m = [13, 14, 15];
399 g = grid.equidistant(m, xlims, ylims, zlims);
400
401 % Quadrature
402 opsx = sbp.D2Standard(m(1), xlims, order);
403 opsy = sbp.D2Standard(m(2), ylims, order);
404 opsz = sbp.D2Standard(m(3), zlims, order);
405 Hx = opsx.H;
406 Hy = opsy.H;
407 Hz = opsz.H;
408 H = {Hx, Hy, Hz};
409
410 % Moment conditions
411 fs = cell(mom_cond,1);
412 for p = 0:mom_cond-1
413 fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p;
414 end
415 end