Line data Source code
1 : module test_support
2 :
3 : use num_def
4 : use num_lib
5 : use math_lib
6 : use const_def, only: dp, arg_not_provided
7 : use utils_lib, only: mesa_error
8 :
9 : implicit none
10 :
11 : contains
12 :
13 31 : subroutine show_results(nv, y, expect, show_all)
14 : integer, intent(in) :: nv
15 : real(dp), dimension(nv), intent(in) :: y, expect
16 : logical, intent(in) :: show_all
17 : integer :: i
18 : include 'formats'
19 31 : if (show_all) then
20 0 : do i = 1, nv
21 0 : write (*, 2) 'calculated', i, y(i)
22 : end do
23 0 : do i = 1, nv
24 0 : write (*, 2) 'reference', i, expect(i)
25 0 : write (*, 2) 'lg(abs rel diff)', i, &
26 0 : log10(abs(y(i) - expect(i))/max(1d-299, abs(expect(i))))
27 : end do
28 : else
29 221 : do i = 1, nv
30 221 : write (*, 2) 'calculated', i, y(i)
31 : end do
32 221 : do i = 1, nv
33 221 : write (*, 2) 'reference', i, expect(i)
34 : end do
35 : end if
36 31 : write (*, *)
37 31 : end subroutine show_results
38 :
39 29 : subroutine show_statistics(nfcn, njac, nstep, show_all)
40 : integer, intent(in) :: nfcn, njac, nstep
41 : logical, intent(in) :: show_all
42 29 : if (.not. show_all) return
43 0 : write (*, *) 'number of steps ', nstep
44 0 : write (*, *) 'number of function evals', nfcn
45 0 : write (*, *) 'number of jacobians ', njac
46 0 : write (*, *) 'functions + jacobians ', nfcn + njac
47 0 : write (*, *)
48 : end subroutine show_statistics
49 :
50 8 : subroutine check_results(nv, y, expect, tol, ierr)
51 : integer, intent(in) :: nv
52 : real(dp), dimension(nv), intent(in) :: y, expect
53 : real(dp), intent(in) :: tol
54 : integer, intent(out) :: ierr
55 : integer :: i
56 : logical :: okay
57 : include 'formats'
58 8 : okay = .true.
59 8 : ierr = 0
60 3208 : do i = 1, nv
61 3200 : if (abs(expect(i)) < tol) cycle
62 1560 : if (abs(y(i) - expect(i)) > tol) then
63 0 : write (*, 2) 'check results result#', i
64 0 : write (*, 1) 'log10 abs diff', log10(abs(y(i) - expect(i)))
65 0 : write (*, 1) 'y(i)', y(i)
66 0 : write (*, 1) 'expect(i)', expect(i)
67 0 : write (*, 1) 'log10 tol', log10(tol)
68 0 : write (*, *)
69 0 : ierr = -1
70 0 : okay = .false.
71 : end if
72 : end do
73 8 : if (okay) return
74 0 : write (*, *)
75 0 : do i = 1, nv
76 0 : write (*, 2) 'y expected', i, y(i), expect(i)
77 : end do
78 0 : write (*, *)
79 0 : call mesa_error(__FILE__, __LINE__)
80 0 : return
81 : end subroutine check_results
82 :
83 8 : real(dp) function f(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
84 : integer, intent(in) :: lrpar, lipar
85 : real(dp), intent(in) :: x
86 : real(dp), intent(out) :: dfdx
87 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
88 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
89 : integer, intent(out) :: ierr
90 8 : ierr = 0
91 8 : f = x - 3*sin(1 - x)
92 8 : dfdx = 1 + 3*cos(1 - x)
93 8 : end function f
94 :
95 1 : subroutine test_root_with_brackets
96 : integer, parameter :: lrpar = 0, lipar = 0
97 : real(dp) :: x, dfdx
98 : real(dp) :: x1, x3 ! bounds for x
99 : ! values of f at x1 and x3 must have opposite sign
100 : ! return value for safe_root will be bracketed by x1 and x3
101 : real(dp) :: y1, y3 ! f(x1) and f(x3)
102 : integer :: imax ! max number of iterations for search
103 : real(dp) :: epsx, epsy
104 : ! stop searching when x is determined to within epsx
105 : ! or when abs(f(x)) is less than epsy
106 : integer :: ierr
107 : real(dp), parameter :: expected_root = 0.74800611d0
108 : real(dp), target :: rpar_ary(lrpar)
109 : integer, target :: ipar_ary(lipar)
110 : integer, pointer :: ipar(:) ! (lipar)
111 : real(dp), pointer :: rpar(:) ! (lrpar)
112 : include 'formats'
113 1 : ipar => ipar_ary
114 1 : rpar => rpar_ary
115 : ierr = 0
116 1 : imax = 100
117 1 : x1 = 0
118 1 : x3 = 10
119 1 : y1 = f(x1, dfdx, lrpar, rpar, lipar, ipar, ierr)
120 1 : y3 = f(x3, dfdx, lrpar, rpar, lipar, ipar, ierr)
121 1 : epsx = 1d-6
122 1 : epsy = 1d-6
123 : x = safe_root_with_brackets( &
124 1 : f, x1, x3, y1, y3, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
125 1 : if (abs(x - expected_root) > 1d-6) call mesa_error(__FILE__, __LINE__)
126 1 : write (*, 1) 'root', x
127 1 : end subroutine test_root_with_brackets
128 :
129 34 : real(dp) function test_f(x, dfdx, lrpar, rpar, lipar, ipar, ierr)
130 : ! returns with ierr = 0 if was able to evaluate f and df/dx at x
131 : real(dp), intent(in) :: x
132 : real(dp), intent(out) :: dfdx
133 : integer, intent(in) :: lipar, lrpar
134 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
135 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
136 : integer, intent(out) :: ierr
137 34 : test_f = tanh(x) - 0.4621171572600098d0
138 34 : dfdx = 1/cosh(x)**2
139 34 : ierr = 0
140 34 : end function test_f
141 :
142 1 : subroutine test_root2
143 : real(dp) :: x ! provide starting guess on input
144 : real(dp) :: x1, x3 ! bounds for x
145 : real(dp) :: y1, y3 ! f(x1) and f(x3)
146 : integer, parameter :: imax = 50, lipar = 0, lrpar = 0
147 : real(dp) :: dx
148 : real(dp), parameter :: epsx = 1d-10, epsy = 1d-10
149 : integer :: ierr
150 : real(dp), target :: rpar_ary(lrpar)
151 : integer, target :: ipar_ary(lipar)
152 : integer, pointer :: ipar(:) ! (lipar)
153 : real(dp), pointer :: rpar(:) ! (lrpar)
154 : include 'formats'
155 1 : ipar => ipar_ary
156 1 : rpar => rpar_ary
157 1 : x = -1d0
158 1 : dx = 0.1d0
159 : ierr = 0
160 1 : write (*, *) 'test_root2'
161 1 : call look_for_brackets(x, dx, x1, x3, test_f, y1, y3, imax, lrpar, rpar, lipar, ipar, ierr)
162 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
163 1 : write (*, 1) 'x1', x1
164 1 : write (*, 1) 'x3', x3
165 1 : write (*, 1) 'y1', y1
166 1 : write (*, 1) 'y3', y3
167 : x = safe_root_with_brackets( &
168 1 : test_f, x1, x3, y1, y3, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
169 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
170 1 : write (*, 1) 'safe_root', x
171 1 : write (*, *)
172 1 : end subroutine test_root2
173 :
174 1 : subroutine test_root3
175 : real(dp) :: x ! provide starting guess on input
176 : real(dp) :: x1, x3 ! bounds for x
177 : real(dp) :: y1, y3 ! f(x1) and f(x3)
178 : integer, parameter :: newt_imax = 10, imax = 50, lipar = 0, lrpar = 0
179 : real(dp) :: dx
180 : real(dp), parameter :: epsx = 1d-10, epsy = 1d-10
181 : integer :: ierr
182 : real(dp), target :: rpar_ary(lrpar)
183 : integer, target :: ipar_ary(lipar)
184 : integer, pointer :: ipar(:) ! (lipar)
185 : real(dp), pointer :: rpar(:) ! (lrpar)
186 : include 'formats'
187 1 : ipar => ipar_ary
188 1 : rpar => rpar_ary
189 1 : dx = 0.1d0
190 1 : x1 = arg_not_provided
191 1 : x3 = arg_not_provided
192 1 : y1 = arg_not_provided
193 1 : y3 = arg_not_provided
194 : ierr = 0
195 1 : write (*, *) 'test_root3'
196 1 : x = 0.1d0 ! not too bad initial guess. newton should find it okay.
197 : x = safe_root_with_guess( &
198 : test_f, x, dx, x1, x3, y1, y3, &
199 1 : newt_imax, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
200 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
201 1 : write (*, 1) 'first safe_root_with_guess', x
202 1 : x = -1d0 ! really bad guess will make it give up on newton
203 : x = safe_root_with_guess( &
204 : test_f, x, dx, x1, x3, y1, y3, &
205 1 : newt_imax, imax, epsx, epsy, lrpar, rpar, lipar, ipar, ierr)
206 1 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
207 1 : write (*, 1) 'second safe_root_with_guess', x
208 1 : write (*, *)
209 1 : end subroutine test_root3
210 :
211 17231 : subroutine van_der_Pol_derivs(n, x, h, y, f, lrpar, rpar, lipar, ipar, ierr)
212 : integer, intent(in) :: n, lrpar, lipar
213 : real(dp), intent(in) :: x, h
214 : real(dp), intent(inout) :: y(:) ! (n)
215 : real(dp), intent(inout) :: f(:) ! (n)
216 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
217 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
218 : integer, intent(out) :: ierr ! nonzero means retry with smaller timestep.
219 : include 'formats'
220 17231 : ierr = 0
221 17231 : f(1) = y(2)
222 17231 : f(2) = ((1 - y(1)*y(1))*y(2) - y(1))/rpar(1)
223 : ! the derivatives do not depend on x
224 17231 : end subroutine van_der_Pol_derivs
225 :
226 1842 : subroutine solout(nr, xold, x, n, y, rwork, iwork, interp_y, lrpar, rpar, lipar, ipar, irtrn)
227 : ! nr is the step number.
228 : ! x is the current x value; xold is the previous x value.
229 : ! y is the current y value.
230 : ! irtrn negative means terminate integration.
231 : ! rwork and iwork hold info for
232 : integer, intent(in) :: nr, n, lrpar, lipar
233 : real(dp), intent(in) :: xold, x
234 : real(dp), intent(inout) :: y(:) ! (n)
235 : real(dp), intent(inout), target :: rwork(*)
236 : integer, intent(inout), target :: iwork(*)
237 : integer, intent(inout), pointer :: ipar(:) ! (lipar)
238 : real(dp), intent(inout), pointer :: rpar(:) ! (lrpar)
239 : interface
240 : ! this subroutine can be called from your solout routine.
241 : ! it computes interpolated values for y components during the just completed step.
242 : real(dp) function interp_y(i, s, rwork, iwork, ierr)
243 : use const_def, only: dp
244 : implicit none
245 : integer, intent(in) :: i ! result is interpolated approximation of y(i) at x=s.
246 : real(dp), intent(in) :: s ! interpolation x value (between xold and x).
247 : real(dp), intent(inout), target :: rwork(*)
248 : integer, intent(inout), target :: iwork(*)
249 : integer, intent(out) :: ierr
250 : end function interp_y
251 : end interface
252 : integer, intent(out) :: irtrn
253 : ! --- prints solution at equidistant output-points
254 : ! --- by using "contd8", the continuous collocation solution
255 1842 : real(dp) :: xout, y1, y2
256 : integer :: ierr
257 1842 : xout = rpar(2)
258 1842 : irtrn = 1
259 1842 : if (ipar(1) /= 1) return ! no output
260 :
261 0 : if (nr == 1) then
262 0 : write (6, 99) x, y(1), y(2), nr - 1
263 0 : xout = x + 0.2d0
264 : else
265 0 : do
266 0 : if (x >= xout - 1d-10) then
267 : ierr = 0
268 0 : y1 = interp_y(1, xout, rwork, iwork, ierr)
269 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
270 0 : y2 = interp_y(2, xout, rwork, iwork, ierr)
271 0 : if (ierr /= 0) call mesa_error(__FILE__, __LINE__)
272 0 : write (6, 99) xout, y1, y2, nr - 1
273 0 : xout = xout + 0.2d0
274 : cycle
275 : end if
276 : exit
277 : end do
278 : end if
279 : 99 format(1x, 'x =', f5.2, ' y =', 2e18.10, ' nstep =', i6)
280 0 : rpar(2) = xout
281 : end subroutine solout
282 :
283 2 : subroutine test_dopri(do_853, show_all)
284 : logical, intent(in) :: do_853, show_all
285 : integer, parameter :: nv = 2 ! the number of variables in the van der Pol system of ODEs
286 : real(dp), parameter :: eps = 1d-3 ! stiffness parameter for van der Pol
287 4 : real(dp) :: rtol(1) ! relative error tolerance(s)
288 4 : real(dp) :: atol(1) ! absolute error tolerance(s)
289 2 : real(dp) :: x ! starting value for the interval of integration
290 2 : real(dp) :: xend ! ending value for the interval of integration
291 6 : real(dp) :: expect(nv)
292 : integer, parameter :: lrpar = 2, lipar = 1, nrdens = nv
293 : integer, parameter :: liwork = nrdens + 100, lwork = 11*nv + 8*nrdens + 100
294 2 : real(dp) :: h, max_step_size
295 : integer :: lout, iout, idid, itol, j
296 : integer :: check_liwork, check_lwork, max_steps, ierr
297 6 : real(dp), target :: y_ary(nv)
298 : real(dp), pointer :: y(:)
299 284 : real(dp), target :: rpar_ary(lrpar), work_ary(lwork)
300 : integer, target :: ipar_ary(lipar), iwork_ary(liwork)
301 : integer, pointer :: ipar(:) ! (lipar)
302 : real(dp), pointer :: rpar(:) ! (lrpar)
303 : integer, pointer :: iwork(:)
304 : real(dp), pointer :: work(:)
305 :
306 2 : ipar => ipar_ary
307 2 : rpar => rpar_ary
308 2 : work => work_ary
309 2 : iwork => iwork_ary
310 2 : y => y_ary
311 :
312 2 : write (*, *)
313 2 : write (*, *) 'vdpol'
314 2 : if (do_853) then
315 1 : write (*, *) 'dop853'
316 : else
317 1 : write (*, *) 'dopri5'
318 : end if
319 :
320 2 : x = 0
321 2 : xend = 2.0
322 2 : y(1) = 2
323 2 : y(2) = 0
324 :
325 2 : lout = 0
326 2 : max_steps = 0
327 2 : max_step_size = 9
328 :
329 2 : itol = 0 ! scalar tolerances
330 2 : iout = 2 ! want dense output
331 :
332 2 : rtol(1) = 1d-4
333 2 : atol(1) = 1d-4
334 2 : h = 1d-6
335 :
336 2 : rpar(1) = eps
337 2 : rpar(2) = 0
338 2 : if (show_all) then
339 0 : ipar(1) = 1
340 : else
341 2 : ipar(1) = 0
342 : end if
343 :
344 206 : iwork = 0
345 278 : work = 0
346 :
347 2 : iwork(5) = nrdens ! want dense output for all components
348 2 : iwork(4) = 1 ! test for stiffness at each step
349 :
350 2 : if (do_853) then
351 1 : call dopri5_work_sizes(nv, nrdens, check_liwork, check_lwork)
352 : else
353 1 : call dop853_work_sizes(nv, nrdens, check_liwork, check_lwork)
354 : end if
355 :
356 2 : if (check_liwork > liwork .or. check_lwork > lwork) then
357 0 : write (*, *) 'need to enlarge work arrays for dopri5'
358 0 : call mesa_error(__FILE__, __LINE__)
359 : end if
360 :
361 2 : ierr = 0
362 2 : if (do_853) then
363 : call dop853( &
364 : nv, van_der_Pol_derivs, x, y, xend, &
365 : h, max_step_size, max_steps, &
366 : rtol, atol, itol, &
367 : solout, iout, work, lwork, iwork, liwork, &
368 1 : lrpar, rpar, lipar, ipar, lout, idid)
369 : else
370 : call dopri5( &
371 : nv, van_der_Pol_derivs, x, y, xend, &
372 : h, max_step_size, max_steps, &
373 : rtol, atol, itol, &
374 : solout, iout, work, lwork, iwork, liwork, &
375 1 : lrpar, rpar, lipar, ipar, lout, idid)
376 : end if
377 :
378 2 : if (idid /= 1) then ! trouble
379 0 : write (*, *) 'idid', idid
380 0 : call mesa_error(__FILE__, __LINE__)
381 : end if
382 :
383 2 : expect(1:2) = [1.7632345401889102d+00, -8.3568868191466206d-01]
384 :
385 2 : call show_results(nv, y, expect, show_all)
386 2 : if (.not. show_all) return
387 :
388 : ! typical: fcn= 21530 step= 1468 accpt= 1345 rejct= 122
389 0 : write (6, 91) (iwork(j), j=17, 20)
390 : 91 format(' fcn=', i8, ' step=', i6, ' accpt=', i6, ' rejct=', i5)
391 :
392 0 : write (*, *)
393 :
394 : end subroutine test_dopri
395 :
396 1 : subroutine test_binary_search
397 : integer, parameter :: n = 100
398 : integer :: k, loc(3)
399 :
400 4 : real(dp) :: val(3)
401 101 : real(dp), target :: vec_ary(n)
402 : real(dp), pointer :: vec(:)
403 : include 'formats'
404 1 : vec => vec_ary
405 :
406 101 : do k = 1, n
407 101 : vec(k) = dble(k)*dble(k)
408 : end do
409 :
410 1 : write (*, *)
411 1 : write (*, *) 'binary_search, increasing values'
412 :
413 : loc = -1
414 4 : val = [0d0, FLOOR(n/3d0)**2 + 2d0, vec(n) + 1d0]
415 4 : do k = 1, 3
416 3 : loc(k) = binary_search(n, vec, 0, val(k))
417 3 : if (loc(k) == 0 .and. val(k) < vec(1)) then
418 1 : write (*, 2) 'val is less than vec(1):', k, val(k), vec(1)
419 2 : else if (loc(k) == n .and. val(k) > vec(n)) then
420 1 : write (*, 1) 'val is greater than vec(n):', val(k), vec(n)
421 1 : else if (vec(loc(k)) <= val(k) .and. val(k) < vec(loc(k) + 1)) then
422 1 : write (*, 1) 'vec(result)', vec(loc(k))
423 1 : write (*, 1) 'val', val(k)
424 1 : write (*, 1) 'vec(result+1)', vec(loc(k) + 1)
425 : else
426 0 : write (*, *) 'binary_search failed for increasing-value array'
427 0 : call mesa_error(__FILE__, __LINE__)
428 : end if
429 4 : write (*, *) 'okay'
430 : end do
431 :
432 : ! test decreasing values
433 : loc = -1
434 101 : where (vec /= 0d0) vec = -vec
435 4 : where (val /= 0d0) val = -val
436 :
437 1 : write (*, *)
438 1 : write (*, *) 'binary_search, decreasing values'
439 4 : do k = 1, 3
440 3 : loc(k) = binary_search(n, vec, 0, val(k))
441 3 : if (loc(k) == 0 .and. val(k) > vec(1)) then
442 1 : write (*, 2) 'val is greater than vec(n):', k, val(k), vec(1)
443 2 : else if (loc(k) == n .and. val(k) < vec(n)) then
444 1 : write (*, 1) 'val is less than vec(1):', val(k), vec(n)
445 1 : else if (vec(loc(k)) >= val(k) .and. val(k) > vec(loc(k) + 1)) then
446 1 : write (*, 1) 'vec(result+1)', vec(loc(k) + 1)
447 1 : write (*, 1) 'val', val(k)
448 1 : write (*, 1) 'vec(result)', vec(loc(k))
449 : else
450 0 : write (*, *) 'binary_search failed for decreasing-value array'
451 0 : call mesa_error(__FILE__, __LINE__)
452 : end if
453 4 : write (*, *) 'okay'
454 : end do
455 1 : write (*, *)
456 :
457 1 : end subroutine test_binary_search
458 :
459 1 : subroutine test_qsort
460 : use const_def
461 : integer, parameter :: n = 100
462 : integer :: ord(n), i
463 101 : real(dp) :: a(n)
464 : include 'formats'
465 1 : write (*, *)
466 1 : write (*, *) 'qsort into increasing order'
467 101 : do i = 1, n
468 101 : a(i) = sinpi(2.1*dble(i)/dble(n))
469 : end do
470 1 : call qsort(ord, n, a)
471 101 : do i = 1, n
472 101 : write (*, 3) 'ord, a(ord)', i, ord(i), a(ord(i))
473 : end do
474 1 : write (*, *)
475 1 : end subroutine test_qsort
476 :
477 2 : real(dp) function g(x) result(y)
478 : real(dp), intent(in) :: x
479 2 : y = (x - 3)*(x - 8)
480 0 : end function g
481 :
482 1 : subroutine test_find0_quadratic
483 : real(dp) :: xx1, yy1, xx2, yy2, xx3, yy3, x, y
484 : integer :: ierr
485 : include 'formats'
486 1 : write (*, *) 'test_find0_quadratic'
487 1 : xx1 = 1
488 1 : yy1 = g(xx1)
489 1 : xx2 = 2
490 1 : yy2 = g(xx2)
491 1 : xx3 = 4
492 1 : yy3 = g(xx3)
493 1 : x = find0_quadratic(xx1, yy1, xx2, yy2, xx3, yy3, ierr)
494 1 : if (ierr /= 0) then
495 0 : call mesa_error(__FILE__, __LINE__)
496 : end if
497 1 : y = g(x)
498 1 : write (*, 1) 'x', x
499 1 : write (*, 1) 'y', y
500 :
501 1 : xx1 = 6
502 1 : yy1 = g(xx1)
503 1 : xx2 = 7
504 1 : yy2 = g(xx2)
505 1 : xx3 = 8
506 1 : yy3 = g(xx3)
507 1 : x = find0_quadratic(xx1, yy1, xx2, yy2, xx3, yy3, ierr)
508 1 : if (ierr /= 0) then
509 0 : call mesa_error(__FILE__, __LINE__)
510 : end if
511 1 : y = g(x)
512 1 : write (*, 1) 'x', x
513 1 : write (*, 1) 'y', y
514 1 : write (*, *)
515 1 : end subroutine test_find0_quadratic
516 :
517 2 : subroutine test_find_max_quadratic
518 1 : real(dp) :: x1, y1, x2, y2, x3, y3, dx1, dx2, xmax, ymax
519 : integer :: ierr
520 : include 'formats'
521 1 : write (*, *) 'test_find_max_quadratic'
522 1 : dx1 = 1; dx2 = 2; y1 = 1; y2 = 10; y3 = 8
523 1 : x1 = 0; x2 = x1 + dx1; x3 = x2 + dx2
524 : ierr = 0
525 1 : call find_max_quadratic(x1, y1, x2, y2, x3, y3, xmax, ymax, ierr)
526 1 : if (ierr /= 0) then
527 0 : call mesa_error(__FILE__, __LINE__)
528 : end if
529 1 : write (*, 1) 'xmax', xmax, 1.85d0
530 1 : write (*, 1) 'ymax', ymax, 12.4083d0
531 1 : write (*, *)
532 1 : end subroutine test_find_max_quadratic
533 :
534 : end module test_support
|
