Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 1234:f1806475498b feature/dirac_discr

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