helpers_m.f90 Source File

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
    implicit none

    public :: digit_sum_csr
    public :: exclude_third_quadrant_circular
    public :: exclude_third_quadrant_cerfons

    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

end module helpers_m