module helpers_m !! Module implementing useful routines for unit testing use precision_m, only : SP, FP, DP, FP_EPS, SP_EPS, DP_EPS use csrmat_m, only: csrmat_t use screen_io_m, only : get_stdout use constants_m, only: pi implicit none public :: digit_sum_csr public :: exclude_third_quadrant_circular public :: exclude_third_quadrant_cerfons public :: excl_3rd_quadrant_segment_circtor interface almost_equal !! Check equality between single or double precision numbers module procedure almost_equal_single module procedure almost_equal_double end interface contains function almost_equal_single(a, b, rtol, atol) result(res) !! Checks whether (a-b) / max(a,eps) within tolerance real(SP), intent(in):: a, b !! values to be compared real(SP) :: rtol !! relative tolerance real(SP) :: atol !! minumum value for normalisation (~absolute tolerance for rtol=0) logical:: res real(SP) :: maxdiff maxdiff = max(abs(rtol*a), atol) if( abs(a - b) < maxdiff) then res = .true. else res = .false. endif if (.not.res) then write(get_stdout(), *) 'almost_equal failed at' write(get_stdout(), *) 'rtol =', rtol,', atol=',atol write(get_stdout(), *) 'a = ',a write(get_stdout(), *) 'b = ',b write(get_stdout(), *) 'a-b = ',a-b endif end function function almost_equal_double(a, b, rtol, atol) result(res) !! Checks whether (a-b) / max(a,eps) within tolerance real(DP), intent(in):: a, b !! values to be compared real(DP) :: rtol !! relative tolerance real(DP) :: atol !! minumum value for normalisation (~absolute tolerance for rtol=0) logical:: res real(DP) :: maxdiff maxdiff = max(abs(rtol*a), atol) if( abs(a - b) < maxdiff) then res = .true. else res = .false. endif if (.not.res) then write(get_stdout(), *) 'almost_equal failed at' write(get_stdout(), *) 'rtol =', rtol,', atol=',atol write(get_stdout(), *) 'a = ',a write(get_stdout(), *) 'b = ',b write(get_stdout(), *) 'a-b = ',a-b endif end function subroutine digit_sum_csr(mat, isum, jsum, valsum) !! Calculates the digit sum of the csr matrix mat type(csrmat_t), intent(in) :: mat !! csr matrix integer, intent(out) :: isum !! sum over i index of matrix integer, intent(out) :: jsum !! sum over j index of matrix real(FP), intent(out) :: valsum !! sum over values of matrix integer :: k isum = 0 do k = 1, mat%ndim+1 isum = isum + mat%i(k) enddo jsum = 0 valsum = 0.0_DP do k = 1, mat%nnz jsum = jsum + mat%j(k) valsum = valsum + mat%val(k) enddo end subroutine logical function exclude_third_quadrant_circular(x, y, phi) !! Excludes third quadrant in field line tracing for circular equilibrium real(FP), intent(in) :: x, y, phi exclude_third_quadrant_circular = .false. if ((x <= 0.0_FP).and. (y <= 0.0_FP)) then exclude_third_quadrant_circular = .true. else exclude_third_quadrant_circular = .false. endif end function logical function exclude_third_quadrant_cerfons(x, y, phi) !! Excludes third quadrant in field line tracing for cerfons equilibrium real(FP), intent(in) :: x, y, phi exclude_third_quadrant_cerfons = .false. if ((x <= 1.0_FP).and. (y <= 0.0_FP)) then exclude_third_quadrant_cerfons = .true. else exclude_third_quadrant_cerfons = .false. endif end function logical function excl_3rd_quadrant_segment_circtor(x, y, phi) !! Excludes third quadrant in field line tracing for the circular !! toroidal equilibrium between phi = [pi/4,pi/2]. real(FP), intent(in) :: x, y, phi excl_3rd_quadrant_segment_circtor = .false. if ((x <= 1.0_FP).and. (y <= 0.0_FP) .and. (phi > 0.25_FP * pi) .and. & (phi < 0.5_FP * pi)) then excl_3rd_quadrant_segment_circtor = .true. else excl_3rd_quadrant_segment_circtor = .false. endif end function end module helpers_m