Line data Source code
1 : module interp_1d_support
2 :
3 : use const_def, only: dp, pi
4 : use math_lib
5 : use interp_1d_def
6 : use interp_1d_lib
7 : use utils_lib, only: mesa_error
8 :
9 : implicit none
10 :
11 : integer, parameter :: num_xpts = 31
12 : integer, parameter :: num_ypts = 35
13 : integer, parameter :: sz_per_pt = 4
14 :
15 : contains
16 :
17 1 : subroutine do_test
18 :
19 1 : write (*, *)
20 1 : call test_min(.true.)
21 :
22 1 : write (*, *)
23 1 : call test_min(.false.)
24 :
25 1 : write (*, *)
26 1 : call test1(.true.)
27 :
28 1 : write (*, *)
29 1 : call test1(.false.)
30 :
31 1 : write (*, *)
32 1 : call test2(.true.)
33 :
34 1 : write (*, *)
35 1 : call test2(.false.)
36 :
37 1 : write (*, *)
38 1 : call test3(.true.)
39 :
40 1 : write (*, *)
41 1 : call test3(.false.)
42 :
43 1 : write (*, *)
44 1 : call test_4pt
45 :
46 1 : write (*, *)
47 1 : call test_interp_3_to_2
48 :
49 1 : end subroutine do_test
50 :
51 6 : subroutine test_data(n, U)
52 : integer, intent(in) :: n
53 : real(dp), intent(out) :: U(n)
54 :
55 6 : U(1) = 0.3691382918814334D0
56 6 : U(2) = 0.3695899829574669D0
57 6 : U(3) = 0.3699633032242006D0
58 6 : U(4) = 0.3709860746618878D0
59 6 : U(5) = 0.3724175838188959D0
60 6 : U(6) = 0.3742482756438662D0
61 6 : U(7) = 0.3764768199547272D0
62 6 : U(8) = 0.3790910321656262D0
63 6 : U(9) = 0.3820871425096383D0
64 6 : U(10) = 0.3854500131142883D0
65 6 : U(11) = 0.3891727128876379D0
66 6 : U(12) = 0.3932377816410259D0
67 6 : U(13) = 0.3976355585696444D0
68 6 : U(14) = 0.4023461172079230D0
69 6 : U(15) = 0.4073570798427897D0
70 6 : U(16) = 0.4126461990878889D0
71 6 : U(17) = 0.4181985492965749D0
72 6 : U(18) = 0.4239897218739848D0
73 6 : U(19) = 0.4300024071530404D0
74 6 : U(20) = 0.4362102270543894D0
75 6 : U(21) = 0.4425936978527515D0
76 6 : U(22) = 0.4491247203356548D0
77 6 : U(23) = 0.4557819260544076D0
78 6 : U(24) = 0.4625358769752868D0
79 6 : U(25) = 0.4693638888826394D0
80 6 : U(26) = 0.4762362700345437D0
81 6 : U(27) = 0.4831315999154387D0
82 6 : U(28) = 0.4900124769496845D0
83 6 : U(29) = 0.4968509249135122D0
84 6 : U(30) = 0.5036231656397859D0
85 6 : U(31) = 0.5102978108992456D0
86 6 : U(32) = 0.5168503329341820D0
87 6 : U(33) = 0.5232493058916751D0
88 6 : U(34) = 0.5294708374154181D0
89 6 : U(35) = 0.5354843215442489D0
90 6 : U(36) = 0.5412671682413143D0
91 6 : U(37) = 0.5467901937111218D0
92 6 : U(38) = 0.5520326792853251D0
93 6 : U(39) = 0.5569674077513700D0
94 6 : U(40) = 0.5615760727117703D0
95 6 : U(41) = 0.5658339305000317D0
96 6 : U(42) = 0.5697255981615711D0
97 6 : U(43) = 0.5732292848999179D0
98 6 : U(44) = 0.5763330504836860D0
99 6 : U(45) = 0.5790185013283228D0
100 6 : U(46) = 0.5812774415045446D0
101 6 : U(47) = 0.5830949206727917D0
102 6 : U(48) = 0.5844671751637984D0
103 6 : U(49) = 0.5853828024764719D0
104 6 : U(50) = 0.5858432670435806D0
105 6 : U(51) = 0.5858432670435451D0
106 6 : U(52) = 0.5853828024764413D0
107 6 : U(53) = 0.5844671751637710D0
108 6 : U(54) = 0.5830949206727681D0
109 6 : U(55) = 0.5812774415045240D0
110 6 : U(56) = 0.5790185013283052D0
111 6 : U(57) = 0.5763330504836706D0
112 6 : U(58) = 0.5732292848999045D0
113 6 : U(59) = 0.5697255981615594D0
114 6 : U(60) = 0.5658339305000216D0
115 6 : U(61) = 0.5615760727117616D0
116 6 : U(62) = 0.5569674077513626D0
117 6 : U(63) = 0.5520326792853185D0
118 6 : U(64) = 0.5467901937111161D0
119 6 : U(65) = 0.5412671682413093D0
120 6 : U(66) = 0.5354843215442447D0
121 6 : U(67) = 0.5294708374154145D0
122 6 : U(68) = 0.5232493058916718D0
123 6 : U(69) = 0.5168503329341793D0
124 6 : U(70) = 0.5102978108992431D0
125 6 : U(71) = 0.5036231656397836D0
126 6 : U(72) = 0.4968509249135103D0
127 6 : U(73) = 0.4900124769496828D0
128 6 : U(74) = 0.4831315999154373D0
129 6 : U(75) = 0.4762362700345427D0
130 6 : U(76) = 0.4693638888826384D0
131 6 : U(77) = 0.4625358769752859D0
132 6 : U(78) = 0.4557819260544067D0
133 6 : U(79) = 0.4491247203356541D0
134 6 : U(80) = 0.4425936978527509D0
135 6 : U(81) = 0.4362102270543888D0
136 6 : U(82) = 0.4300024071530398D0
137 6 : U(83) = 0.4239897218739841D0
138 6 : U(84) = 0.4181985492965744D0
139 6 : U(85) = 0.4126461990878886D0
140 6 : U(86) = 0.4073570798427895D0
141 6 : U(87) = 0.4023461172079228D0
142 6 : U(88) = 0.3976355585696443D0
143 6 : U(89) = 0.3932377816410258D0
144 6 : U(90) = 0.3891727128876379D0
145 6 : U(91) = 0.3854500131142883D0
146 6 : U(92) = 0.3820871425096385D0
147 6 : U(93) = 0.3790910321656264D0
148 6 : U(94) = 0.3764768199547278D0
149 6 : U(95) = 0.3742482756438669D0
150 6 : U(96) = 0.3724175838188968D0
151 6 : U(97) = 0.3709860746618888D0
152 6 : U(98) = 0.3699633032242018D0
153 6 : U(99) = 0.3695899829574566D0
154 6 : U(100) = 0.3691382918814347D0
155 :
156 6 : end subroutine test_data
157 :
158 2 : subroutine test1(increasing)
159 : logical, intent(in) :: increasing
160 : integer, parameter :: n = 100
161 202 : real(dp) :: U(n)
162 :
163 : integer, parameter :: nvals = 4
164 : integer, parameter :: nwork = max(pm_work_size, mp_work_size)
165 2802 : real(dp), target :: work_ary(n*nwork)
166 2 : real(dp), pointer :: work(:)
167 222 : real(dp) :: init_xs(n), xs(nvals), vals(nvals), dx
168 : integer :: i, ierr
169 802 : real(dp), target :: f_ary(4*n)
170 2 : real(dp), pointer :: f1(:), f(:, :)
171 2 : f1 => f_ary
172 2 : f(1:4, 1:n) => f1(1:4*n)
173 :
174 : include 'formats'
175 :
176 2 : write (*, *) 'test1', increasing
177 :
178 2 : work => work_ary
179 2 : dx = 0.594059405940594d0
180 :
181 202 : do i = 1, n
182 202 : init_xs(i) = dx*(i - 1)
183 : end do
184 :
185 102 : if (.not. increasing) init_xs(:) = -init_xs(:)
186 :
187 2 : call test_data(n, U)
188 :
189 202 : f(1, 1:n) = U(1:n)
190 2 : call interp_m3q(init_xs, n, f1, nwork, work, 'test1', ierr)
191 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
192 2 : write (*, 1) 'f(2, 1)', f(2, 1)
193 2 : write (*, 1) 'f(2, 2)', f(2, 2)
194 2 : write (*, 1) 'f(2, 3)', f(2, 3)
195 2 : write (*, *)
196 :
197 2 : xs = [10d0, 10.4d0, 10.8d0, 11d0]
198 6 : if (.not. increasing) xs(:) = -xs(:)
199 2 : call interp_values(init_xs, n, f1, nvals, xs, vals, ierr)
200 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
201 2 : write (*, 1) 'z(10.0)', vals(1)
202 2 : write (*, 1) 'z(10.4)', vals(2)
203 2 : write (*, 1) 'z(10.8)', vals(3)
204 2 : write (*, 1) 'z(11.0)', vals(4)
205 2 : write (*, *)
206 :
207 2 : end subroutine test1
208 :
209 2 : subroutine test2(increasing)
210 : logical, intent(in) :: increasing
211 : integer, parameter :: n = 43
212 88 : real(dp) :: U(n)
213 :
214 : logical, parameter :: show_errors = .false.
215 : integer, parameter :: nvals = 7
216 120 : real(dp) :: init_xs(n), xs(nvals), vals(nvals), dx
217 : integer, parameter :: nwork = max(pm_work_size, mp_work_size)
218 1206 : real(dp), target :: work_ary(n*nwork)
219 2 : real(dp), pointer :: work(:)
220 : integer :: i, ierr
221 346 : real(dp), target :: f_ary(4*n)
222 2 : real(dp), pointer :: f1(:), f(:, :)
223 2 : f1 => f_ary
224 2 : f(1:4, 1:n) => f1(1:4*n)
225 :
226 : include 'formats'
227 :
228 2 : write (*, *) 'test2', increasing
229 :
230 2 : work => work_ary
231 2 : dx = 1d0/(n - 1)
232 88 : do i = 1, n
233 86 : init_xs(i) = 4*dx*(i - 1)
234 86 : if (.not. increasing) init_xs(i) = -init_xs(i)
235 88 : U(i) = sin(init_xs(i))
236 : end do
237 :
238 88 : f(1, 1:n) = U(1:n)
239 2 : call interp_m3q(init_xs, n, f1, nwork, work, 'test2', ierr)
240 2 : if (ierr /= 0) then
241 0 : call mesa_error(__FILE__, __LINE__)
242 : end if
243 2 : write (*, *)
244 :
245 2 : xs = [0d0, pi/4, pi/3, pi/2, 2*pi/3, 3*pi/4, 19*pi/20]
246 9 : if (.not. increasing) xs(:) = -xs(:)
247 2 : call interp_values(init_xs, n, f1, nvals, xs, vals, ierr)
248 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
249 :
250 : if (show_errors) then
251 : do i = 1, nvals
252 : write (*, 2) 'x val exact err', i, xs(i), vals(i), sin(xs(i)), vals(i) - sin(xs(i))
253 : end do
254 : else
255 16 : do i = 1, nvals
256 16 : write (*, 2) 'x val exact', i, xs(i), vals(i), sin(xs(i))
257 : end do
258 : end if
259 2 : write (*, *)
260 :
261 2 : call integrate_values(init_xs, n, f1, nvals, xs, vals, ierr)
262 2 : if (ierr /= 0) then
263 0 : call mesa_error(__FILE__, __LINE__)
264 : end if
265 :
266 : if (show_errors) then
267 : do i = 2, nvals
268 : write (*, 2) 'x val exact err', i, xs(i), vals(i), cos(xs(i - 1)) - cos(xs(i)), vals(i) - (cos(xs(i - 1)) - cos(xs(i)))
269 : end do
270 : else
271 14 : do i = 2, nvals
272 14 : write (*, 2) 'x val exact', i, xs(i), vals(i), cos(xs(i - 1)) - cos(xs(i))
273 : end do
274 : end if
275 2 : write (*, *)
276 :
277 2 : end subroutine test2
278 :
279 2 : subroutine test3(increasing)
280 : logical, intent(in) :: increasing
281 : integer, parameter :: n = 100
282 202 : real(dp) :: U(n)
283 :
284 : integer, parameter :: nvals = 3
285 : integer, parameter :: nwork = max(pm_work_size, mp_work_size)
286 2802 : real(dp), target :: work_ary(n*nwork)
287 2 : real(dp), pointer :: work(:)
288 218 : real(dp) :: init_xs(n), xs(nvals), vals(nvals), dx
289 : integer :: i, ierr
290 :
291 802 : real(dp), target :: f_ary(4*n)
292 2 : real(dp), pointer :: f1(:), f(:, :)
293 2 : f1 => f_ary
294 2 : f(1:4, 1:n) => f1(1:4*n)
295 :
296 : include 'formats'
297 :
298 2 : write (*, *) 'test3', increasing
299 :
300 2 : work => work_ary
301 2 : dx = 0.594059405940594d0
302 :
303 202 : do i = 1, n
304 202 : init_xs(i) = dx*(i - 1)
305 : end do
306 :
307 102 : if (.not. increasing) init_xs(:) = -init_xs(:)
308 :
309 2 : call test_data(n, U)
310 :
311 202 : f(1, 1:n) = U(1:n)
312 2 : call interp_pm(init_xs, n, f1, nwork, work, 'test3', ierr)
313 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
314 2 : write (*, 1) 'f(2, 1)', f(2, 1)
315 2 : write (*, 1) 'f(2, 2)', f(2, 2)
316 2 : write (*, 1) 'f(2, 3)', f(2, 3)
317 2 : write (*, *)
318 :
319 2 : xs = [10d0, 10.4d0, 10.8d0]
320 5 : if (.not. increasing) xs(:) = -xs(:)
321 2 : call interp_values(init_xs, n, f1, nvals, xs, vals, ierr)
322 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
323 2 : write (*, 1) 'z(10.0)', vals(1)
324 2 : write (*, 1) 'z(10.4)', vals(2)
325 2 : write (*, 1) 'z(10.8)', vals(3)
326 2 : write (*, *)
327 :
328 2 : end subroutine test3
329 :
330 2 : subroutine test_min(increasing)
331 : logical, intent(in) :: increasing
332 : integer, parameter :: n = 100
333 202 : real(dp) :: U(n)
334 :
335 : integer, parameter :: nvals = 2
336 : integer, parameter :: nwork = max(pm_work_size, mp_work_size)
337 2802 : real(dp), target :: work_ary(n*nwork)
338 2 : real(dp), pointer :: work(:)
339 214 : real(dp) :: init_xs(n), xs(nvals), vals(nvals), dx
340 : integer :: i, ierr
341 :
342 802 : real(dp), target :: f_ary(4*n)
343 2 : real(dp), pointer :: f1(:), f(:, :)
344 2 : f1 => f_ary
345 2 : f(1:4, 1:n) => f1(1:4*n)
346 :
347 : include 'formats'
348 :
349 2 : write (*, *) 'test3', increasing
350 :
351 2 : dx = 0.594059405940594d0
352 2 : work => work_ary
353 :
354 202 : do i = 1, n
355 202 : init_xs(i) = dx*(i - 1)
356 : end do
357 :
358 102 : if (.not. increasing) init_xs(:) = -init_xs(:)
359 :
360 2 : call test_data(n, U)
361 :
362 202 : f(1, 1:n) = U(1:n)
363 2 : call interp_pm(init_xs, n, f1, nwork, work, 'test_min', ierr)
364 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
365 2 : write (*, 1) 'f(2, 1)', f(2, 1)
366 2 : write (*, 1) 'f(2, 2)', f(2, 2)
367 2 : write (*, *)
368 :
369 2 : xs = [10d0, 10.8d0]
370 4 : if (.not. increasing) xs(:) = -xs(:)
371 2 : call interp_values(init_xs, n, f1, nvals, xs, vals, ierr)
372 2 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
373 2 : write (*, 1) 'z(10.0)', vals(1)
374 2 : write (*, 1) 'z(10.8)', vals(2)
375 2 : write (*, *)
376 :
377 2 : end subroutine test_min
378 :
379 2 : subroutine test_interp_3_to_2
380 1 : real(dp) :: pdqm1, pdq00, ndqm1, ndq00, pfm1, pf00, pfp1, nf00, nfp1
381 : integer :: ierr
382 :
383 : include 'formats'
384 1 : write (*, *) 'test_interp_3_to_2'
385 1 : pdqm1 = 2.0114947208182310D-03
386 1 : pdq00 = 2.0097307373083619D-03
387 1 : ndqm1 = 1.5784933404892154D-03
388 1 : ndq00 = 1.3270582904786616D-03
389 1 : pfm1 = 1.8984458478374221D+30
390 1 : pf00 = 1.4579738233866260D+30
391 1 : pfp1 = 1.1136297079131214D+30
392 1 : call interp_3_to_2(pdqm1, pdq00, ndqm1, ndq00, pfm1, pf00, pfp1, nf00, nfp1, 'test_interp_3_to_2', ierr)
393 1 : write (*, *) 'ierr', ierr
394 1 : write (*, 1) 'nf00', nf00
395 1 : write (*, 1) 'nfp1', nfp1
396 1 : end subroutine test_interp_3_to_2
397 :
398 1 : subroutine test_4pt
399 : integer, parameter :: n = 6
400 17 : real(dp) :: x(n), y(n), a(3), dx, result, exact
401 : integer :: i
402 : integer, parameter :: nwork = max(pm_work_size, mp_work_size)
403 85 : real(dp), target :: work_ary(n*nwork)
404 1 : real(dp), pointer :: work(:)
405 25 : real(dp), target :: f_ary(4*n)
406 1 : real(dp), pointer :: f1(:)
407 :
408 : include 'formats'
409 1 : f1 => f_ary
410 1 : work => work_ary
411 1 : x = [pi/4, pi/3, pi/2, 2*pi/3, 3*pi/4, 19*pi/20]
412 7 : do i = 1, n
413 7 : y(i) = sin(x(i))
414 : end do
415 1 : call interp_4pt_pm(x(2:5), y(2:5), a)
416 1 : dx = (x(4) - x(3))/2
417 1 : result = y(3) + dx*(a(1) + dx*(a(2) + dx*a(3)))
418 1 : exact = sin(x(3) + dx)
419 1 : write (*, *) 'test_4pt'
420 1 : write (*, 1) 'res exact err', result, exact, abs((result - exact)/exact)
421 :
422 1 : end subroutine test_4pt
423 :
424 : end module interp_1d_support
|
