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