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