Fix collision precompute initialization and CQUAD IEEE handling

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

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

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)
 

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

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

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