diff --git a/COMMON/collop_compute.f90 b/COMMON/collop_compute.f90 index 7a3531d4..8f611279 100644 --- a/COMMON/collop_compute.f90 +++ b/COMMON/collop_compute.f90 @@ -56,8 +56,8 @@ module collop_compute !********************************************************** ! Species !********************************************************** - real(kind=dp) :: m_a - real(kind=dp) :: m_b + real(kind=dp) :: m_a = 1.0d0 + real(kind=dp) :: m_b = 1.0d0 character(len=3) :: tag_a character(len=3) :: tag_b @@ -955,7 +955,7 @@ end function func1d if (lsw_split_interval) then do k = 1, kmax-1 - n_sub = num_sub_intervals*max(int(log(m_b/m_a)), 1) + n_sub = scaled_num_sub_intervals() if ((.not. binknots) .and. (k .eq. kmax-1)) n_sub = num_sub_intervals_cutoff x_sub_del = (x_inter(k+1) - x_inter(k))/dble(n_sub) do k_sub = 1,n_sub @@ -1033,7 +1033,7 @@ end function func1d if (lsw_split_interval) then do k = 1, kmax-1 - n_sub = num_sub_intervals*max(int(log(m_b/m_a)), 1) + n_sub = scaled_num_sub_intervals() if ((.not. binknots) .and. (k .eq. kmax-1)) n_sub = num_sub_intervals_cutoff x_sub_del = (x_inter(k+1) - x_inter(k))/dble(n_sub) do k_sub = 1,n_sub @@ -1057,6 +1057,10 @@ end function func1d end function integrate_a_to_b_param + integer function scaled_num_sub_intervals() result(n_sub) + n_sub = num_sub_intervals*max(int(log(m_b/m_a)), 1) + end function scaled_num_sub_intervals + recursive function integrate_a_to_infinity(func1d, a) result(y) real(kind=dp) :: a real(kind=dp) :: y diff --git a/COMMON/gsl_integration_routines_mod.f90 b/COMMON/gsl_integration_routines_mod.f90 index 4fc2be95..0d3cda41 100644 --- a/COMMON/gsl_integration_routines_mod.f90 +++ b/COMMON/gsl_integration_routines_mod.f90 @@ -26,6 +26,8 @@ MODULE gsl_integration_routines_mod USE fgsl ! Fortran interface of the GSL Library USE, INTRINSIC :: iso_c_binding + USE, INTRINSIC :: ieee_arithmetic, ONLY: ieee_get_flag, ieee_get_halting_mode, & + ieee_invalid, ieee_is_finite, ieee_set_flag, ieee_set_halting_mode IMPLICIT NONE @@ -99,6 +101,7 @@ END FUNCTION func1d_param1 ! CQUAD doubly-adaptive integration PUBLIC fint1d_cquad PRIVATE fint1d_param0_cquad, fint1d_param1_cquad + PRIVATE begin_cquad_ieee, end_cquad_ieee, record_cquad_integrand INTERFACE fint1d_cquad MODULE PROCEDURE fint1d_param0_cquad, fint1d_param1_cquad END INTERFACE fint1d_cquad @@ -113,6 +116,38 @@ END FUNCTION func1d_param1 CONTAINS + SUBROUTINE begin_cquad_ieee(invalid_before, invalid_halting) + LOGICAL, INTENT(OUT) :: invalid_before, invalid_halting + + ! CQUAD raises invalid internally for exact-zero subintervals. Callbacks + ! record integrand exceptions separately, so invalid physics remains fatal. + CALL ieee_get_flag(ieee_invalid, invalid_before) + CALL ieee_get_halting_mode(ieee_invalid, invalid_halting) + CALL ieee_set_halting_mode(ieee_invalid, .FALSE.) + CALL ieee_set_flag(ieee_invalid, .FALSE.) + END SUBROUTINE begin_cquad_ieee + + SUBROUTINE end_cquad_ieee(invalid_before, invalid_halting, user_invalid, ra, rda) + LOGICAL, INTENT(IN) :: invalid_before, invalid_halting, user_invalid + REAL(fgsl_double), INTENT(IN) :: ra, rda + + CALL ieee_set_flag(ieee_invalid, invalid_before) + CALL ieee_set_halting_mode(ieee_invalid, invalid_halting) + IF(user_invalid) ERROR STOP 'CQUAD integrand raised IEEE invalid' + IF((.NOT. ieee_is_finite(ra)) .OR. (.NOT. ieee_is_finite(rda))) & + ERROR STOP 'CQUAD returned a nonfinite result' + END SUBROUTINE end_cquad_ieee + + SUBROUTINE record_cquad_integrand(value, user_invalid) + REAL(c_double), INTENT(IN) :: value + LOGICAL, INTENT(INOUT) :: user_invalid + LOGICAL :: invalid_raised + + CALL ieee_get_flag(ieee_invalid, invalid_raised) + user_invalid = user_invalid .OR. invalid_raised .OR. (.NOT. ieee_is_finite(value)) + CALL ieee_set_flag(ieee_invalid, .FALSE.) + END SUBROUTINE record_cquad_integrand + !--------------------------------------------------------------------------------------! ! Q(uadrature) A(daptive) G(eneral integrand) integration procedure RECURSIVE FUNCTION fint1d_param0_qag(func1d_param0_user,x_low,x_up,epsabs,epsrel,sw_qag_rule) result(res) @@ -624,6 +659,7 @@ END FUNCTION func1d_param0_user TYPE(fgsl_integration_cquad_workspace) :: integ_cq TYPE(c_ptr) :: param0_ptr ! This pointer holds the C-location of user-specified ! parameter (here is no parameter specified -> c_null_ptr) + LOGICAL :: invalid_before, invalid_halting, user_invalid ! Input/output parameter of the fgsl_integration_cquad function, ! thus better initialize. @@ -643,8 +679,11 @@ END FUNCTION func1d_param0_user ! Initialize solver 'fgsl_integration_cquad' to use the function 'stdfunc' and ! the user-specified parameters + CALL begin_cquad_ieee(invalid_before, invalid_halting) + user_invalid = .FALSE. statusval = fgsl_integration_cquad(stdfunc, x_low, x_up, & epsabs, epsrel, integ_cq, ra, rda, neval) + CALL end_cquad_ieee(invalid_before, invalid_halting, user_invalid, ra, rda) CALL check_error(statusval) ! Return the results (ra,rda,neval) @@ -672,7 +711,9 @@ RECURSIVE FUNCTION f2c_wrapper_func1d_param0(x, params) BIND(c) IF(C_ASSOCIATED(params)) STOP '***Error*** in f2c_wrapper_func1d_param0' ! Wrap user-specified function to a C-interoperable function + CALL ieee_set_flag(ieee_invalid, .FALSE.) f2c_wrapper_func1d_param0 = func1d_param0_user(x) + CALL record_cquad_integrand(f2c_wrapper_func1d_param0, user_invalid) END FUNCTION f2c_wrapper_func1d_param0 END FUNCTION fint1d_param0_cquad @@ -705,6 +746,7 @@ END FUNCTION func1d_param1_user TYPE(fgsl_integration_cquad_workspace) :: integ_cq TYPE(c_ptr) :: param1_ptr ! This pointer holds the C-location of user-specified ! ! parameter + LOGICAL :: invalid_before, invalid_halting, user_invalid ! Turn off error handler std = fgsl_set_error_handler_off() @@ -720,8 +762,11 @@ END FUNCTION func1d_param1_user ! Initialize solver 'fgsl_integration_cquad' to use the function 'stdfunc' and ! the user-specified parameters + CALL begin_cquad_ieee(invalid_before, invalid_halting) + user_invalid = .FALSE. status = fgsl_integration_cquad(stdfunc, x_low, x_up, & epsabs, epsrel, integ_cq, ra, rda, neval) + CALL end_cquad_ieee(invalid_before, invalid_halting, user_invalid, ra, rda) CALL check_error(status) ! Return the results (ra,rda,neval) @@ -753,7 +798,9 @@ RECURSIVE FUNCTION f2c_wrapper_func1d_param1(x, params) BIND(c) ! Cast C-pointer to the above-defined Fortran pointer CALL C_F_POINTER(params, p) ! Wrap user-specified function to a C-interoperable function + CALL ieee_set_flag(ieee_invalid, .FALSE.) f2c_wrapper_func1d_param1 = func1d_param1_user(x,p) + CALL record_cquad_integrand(f2c_wrapper_func1d_param1, user_invalid) END FUNCTION f2c_wrapper_func1d_param1 END FUNCTION fint1d_param1_cquad diff --git a/TEST/CMakeLists.txt b/TEST/CMakeLists.txt index 41e6b6ba..48ed2cfb 100644 --- a/TEST/CMakeLists.txt +++ b/TEST/CMakeLists.txt @@ -230,6 +230,45 @@ set_tests_properties(collop_zero_endpoint_test PROPERTIES TIMEOUT 30 ) +add_executable(test_collop_precompute_mass test_collop_precompute_mass.f90) +target_link_libraries(test_collop_precompute_mass common) + +add_test(NAME collop_precompute_mass_test + COMMAND test_collop_precompute_mass + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(collop_precompute_mass_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + +add_executable(test_collop_bspline_precompute test_collop_bspline_precompute.f90) +target_link_libraries(test_collop_bspline_precompute common) + +add_test(NAME collop_bspline_precompute_test + COMMAND test_collop_bspline_precompute + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(collop_bspline_precompute_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL;Floating point exception" + TIMEOUT 30 +) + +add_executable(test_cquad_ieee_guard test_cquad_ieee_guard.f90) +target_link_libraries(test_cquad_ieee_guard common) + +add_test(NAME cquad_ieee_guard_test + COMMAND test_cquad_ieee_guard + WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}) + +set_tests_properties(cquad_ieee_guard_test PROPERTIES + PASS_REGULAR_EXPRESSION "All tests passed!" + FAIL_REGULAR_EXPRESSION "FAIL" + TIMEOUT 30 +) + add_executable(test_spline_sparse_assembly test_spline_sparse_assembly.f90) target_link_libraries(test_spline_sparse_assembly common) diff --git a/TEST/test_collop_bspline_precompute.f90 b/TEST/test_collop_bspline_precompute.f90 new file mode 100644 index 00000000..62c49429 --- /dev/null +++ b/TEST/test_collop_bspline_precompute.f90 @@ -0,0 +1,25 @@ +program test_collop_bspline_precompute + use, intrinsic :: ieee_arithmetic, only: ieee_get_flag, ieee_invalid, & + ieee_set_flag, ieee_set_halting_mode + use collisionality_mod, only: collop_bspline_dist, collop_bspline_order, & + collop_bspline_taylor, phi_x_max + use collop_compute, only: init_collop + use rkstep_mod, only: lag, leg + implicit none + + logical :: invalid_raised + + lag = 3 + leg = 3 + phi_x_max = 4.0d0 + collop_bspline_order = 2 + collop_bspline_dist = 1.0d0 + collop_bspline_taylor = .true. + call ieee_set_flag(ieee_invalid, .false.) + call ieee_set_halting_mode(ieee_invalid, .true.) + call init_collop(11, 11, 0.0d0, 0.0d0) + call ieee_get_flag(ieee_invalid, invalid_raised) + + if (invalid_raised) error stop 'FAIL: B-spline precompute raised IEEE invalid' + print *, 'All tests passed!' +end program test_collop_bspline_precompute diff --git a/TEST/test_collop_precompute_mass.f90 b/TEST/test_collop_precompute_mass.f90 new file mode 100644 index 00000000..b288b7b0 --- /dev/null +++ b/TEST/test_collop_precompute_mass.f90 @@ -0,0 +1,18 @@ +program test_collop_precompute_mass + use, intrinsic :: ieee_arithmetic, only: ieee_get_flag, ieee_invalid, & + ieee_set_flag, ieee_set_halting_mode + use collop_compute, only: scaled_num_sub_intervals + implicit none + + integer :: n_sub + logical :: invalid_raised + + call ieee_set_flag(ieee_invalid, .false.) + call ieee_set_halting_mode(ieee_invalid, .true.) + n_sub = scaled_num_sub_intervals() + call ieee_get_flag(ieee_invalid, invalid_raised) + + if (invalid_raised) error stop 'FAIL: default subdivision used an undefined mass ratio' + if (n_sub /= 5) error stop 'FAIL: default subdivision is not the equal-mass value' + print *, 'All tests passed!' +end program test_collop_precompute_mass diff --git a/TEST/test_cquad_ieee_guard.f90 b/TEST/test_cquad_ieee_guard.f90 new file mode 100644 index 00000000..1fb305f9 --- /dev/null +++ b/TEST/test_cquad_ieee_guard.f90 @@ -0,0 +1,50 @@ +program test_cquad_ieee_guard + use, intrinsic :: ieee_arithmetic, only: ieee_get_flag, ieee_invalid, & + ieee_set_flag, ieee_set_halting_mode + use gsl_integration_routines_mod, only: fint1d_cquad + implicit none + + character(len=1024) :: argument, command, executable + integer :: command_status, exit_status + logical :: invalid_raised + real(kind(1.0d0)) :: result(2) + + call get_command_argument(1, argument) + if (trim(argument) == 'invalid-child') then + result = fint1d_cquad(invalid_integrand, 0.0d0, 1.0d0, 1.0d-13, 1.0d-13) + error stop 'FAIL: CQUAD accepted an invalid integrand' + end if + + call ieee_set_flag(ieee_invalid, .false.) + call ieee_set_halting_mode(ieee_invalid, .true.) + result = fint1d_cquad(zero_integrand, 0.0d0, 1.0d0, 1.0d-13, 1.0d-13) + call ieee_get_flag(ieee_invalid, invalid_raised) + if (invalid_raised) error stop 'FAIL: CQUAD leaked an internal invalid flag' + if (any(result /= 0.0d0)) error stop 'FAIL: zero integrand gave a nonzero result' + + call get_command_argument(0, executable) + command = '"' // trim(executable) // '" invalid-child' + call execute_command_line(trim(command), wait=.true., exitstat=exit_status, & + cmdstat=command_status) + if (command_status /= 0) error stop 'FAIL: invalid-integrand child did not launch' + if (exit_status == 0) error stop 'FAIL: invalid integrand was not rejected' + print *, 'All tests passed!' + +contains + + function zero_integrand(x) result(value) + real(kind(1.0d0)) :: x + real(kind(1.0d0)) :: value + + value = 0.0d0*x + end function zero_integrand + + function invalid_integrand(x) result(value) + real(kind(1.0d0)) :: x + real(kind(1.0d0)), volatile :: denominator + real(kind(1.0d0)) :: value + + denominator = x - x + value = denominator/denominator + end function invalid_integrand +end program test_cquad_ieee_guard