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