Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 754:5264ce57b573 feature/d1_staggered

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Bugfix diracDiscr to make it work for grids that are non-equidistant near boundaries.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 15 Jun 2018 14:01:13 -0700
parents 4ee7d15bd8e6
children c70131daaa6e
comparison
equal deleted inserted replaced
753:44c46bd6913a 754:5264ce57b573
218 end 218 end
219 end 219 end
220 end 220 end
221 end 221 end
222 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] = setupStuff(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
332 % ============== Setup functions =======================
223 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond) 333 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
224 334
225 % Grid 335 % Grid
226 xl = -3; 336 xl = -3;
227 xr = 2; 337 xr = 900;
338 L = xr-xl;
228 m = 101; 339 m = 101;
229 h = (xr-xl)/(m-1); 340 h = (xr-xl)/(m-1);
230 g = grid.equidistant(m, {xl, xr}); 341 g = grid.equidistant(m, {xl, xr});
231 x = g.points(); 342 x = g.points();
232 343
235 H = ops.H; 346 H = ops.H;
236 347
237 % Moment conditions 348 % Moment conditions
238 fs = cell(mom_cond,1); 349 fs = cell(mom_cond,1);
239 for p = 0:mom_cond-1 350 for p = 0:mom_cond-1
240 fs{p+1} = @(x) x.^p; 351 fs{p+1} = @(x) (x/L).^p;
241 end 352 end
242 353
243 end 354 end
244 355
245 356 function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond)
357
358 % Grid
359 xl = -3;
360 xr = 900;
361 L = xr-xl;
362 m = 101;
363 [~, g_dual] = grid.primalDual1D(m, {xl, xr});
364 x = g_dual.points();
365 h = g_dual.h;
366
367 % Quadrature
368 ops = sbp.D1Staggered(m, {xl, xr}, order);
369 H = ops.H_dual;
370
371 % Moment conditions
372 fs = cell(mom_cond,1);
373 for p = 0:mom_cond-1
374 fs{p+1} = @(x) (x/L).^p;
375 end
376
377 end
378
379