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