Mercurial > repos > public > sbplib

comparison diracDiscrTest.m @ 836:049e4c6fa630 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add dirac delta function, with corresponding test. Commented out tests for staggered grid because they do not yet exist on this branch.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 05 Sep 2018 14:46:39 -0700
parents
children 43a1c3ac07b1
comparison
equal deleted inserted replaced
811:9d4246ac94c0 836:049e4c6fa630
1 function tests = diracDiscrTest()
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] = setupStuff(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
32 function testLeftRandom(testCase)
33
34 orders = [2, 4, 6];
35 mom_conds = orders;
36
37 for o = 1:length(orders)
38 order = orders(o);
39 mom_cond = mom_conds(o);
40 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
41
42 % Test random points near left boundary
43 x0s = xl + 2*h*rand(1,10);
44
45 for j = 1:length(fs)
46 f = fs{j};
47 fx = f(x);
48 for i = 1:length(x0s)
49 x0 = x0s(i);
50 delta = diracDiscr(x0, x, mom_cond, 0, H);
51 integral = delta'*H*fx;
52 err = abs(integral - f(x0));
53 testCase.verifyLessThan(err, 1e-12);
54 end
55 end
56 end
57 end
58
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] = setupStuff(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)
87
88 orders = [2, 4, 6];
89 mom_conds = orders;
90
91 for o = 1:length(orders)
92 order = orders(o);
93 mom_cond = mom_conds(o);
94 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
95
96 % Test random points near right boundary
97 x0s = xr - 2*h*rand(1,10);
98
99 for j = 1:length(fs)
100 f = fs{j};
101 fx = f(x);
102 for i = 1:length(x0s)
103 x0 = x0s(i);
104 delta = diracDiscr(x0, x, mom_cond, 0, H);
105 integral = delta'*H*fx;
106 err = abs(integral - f(x0));
107 testCase.verifyLessThan(err, 1e-12);
108 end
109 end
110 end
111 end
112
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] = setupStuff(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)
142
143 orders = [2, 4, 6];
144 mom_conds = orders;
145
146 for o = 1:length(orders)
147 order = orders(o);
148 mom_cond = mom_conds(o);
149 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
150
151 % Test random points in interior
152 x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20);
153
154 for j = 1:length(fs)
155 f = fs{j};
156 fx = f(x);
157 for i = 1:length(x0s)
158 x0 = x0s(i);
159 delta = diracDiscr(x0, x, mom_cond, 0, H);
160 integral = delta'*H*fx;
161 err = abs(integral - f(x0));
162 testCase.verifyLessThan(err, 1e-12);
163 end
164 end
165 end
166 end
167
168 % x0 outside grid should yield 0 integral!
169 function testX0OutsideGrid(testCase)
170
171 orders = [2, 4, 6];
172 mom_conds = orders;
173
174 for o = 1:length(orders)
175 order = orders(o);
176 mom_cond = mom_conds(o);
177 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
178
179 % Test points outisde grid
180 x0s = [xl-1.1*h, xr+1.1*h];
181
182 for j = 1:length(fs)
183 f = fs{j};
184 fx = f(x);
185 for i = 1:length(x0s)
186 x0 = x0s(i);
187 delta = diracDiscr(x0, x, mom_cond, 0, H);
188 integral = delta'*H*fx;
189 err = abs(integral - 0);
190 testCase.verifyLessThan(err, 1e-12);
191 end
192 end
193 end
194 end
195
196 function testAllGP(testCase)
197
198 orders = [2, 4, 6];
199 mom_conds = orders;
200
201 for o = 1:length(orders)
202 order = orders(o);
203 mom_cond = mom_conds(o);
204 [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond);
205
206 % Test all grid points
207 x0s = x;
208
209 for j = 1:length(fs)
210 f = fs{j};
211 fx = f(x);
212 for i = 1:length(x0s)
213 x0 = x0s(i);
214 delta = diracDiscr(x0, x, mom_cond, 0, H);
215 integral = delta'*H*fx;
216 err = abs(integral - f(x0));
217 testCase.verifyLessThan(err, 1e-12);
218 end
219 end
220 end
221 end
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 =======================
333 function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond)
334
335 % Grid
336 xl = -3;
337 xr = 900;
338 L = xr-xl;
339 m = 101;
340 h = (xr-xl)/(m-1);
341 g = grid.equidistant(m, {xl, xr});
342 x = g.points();
343
344 % Quadrature
345 ops = sbp.D2Standard(m, {xl, xr}, order);
346 H = ops.H;
347
348 % Moment conditions
349 fs = cell(mom_cond,1);
350 for p = 0:mom_cond-1
351 fs{p+1} = @(x) (x/L).^p;
352 end
353
354 end
355
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