diff --git a/autotest/TestMathUtil.f90 b/autotest/TestMathUtil.f90 index 199b6d9654c..331aa0d5c61 100644 --- a/autotest/TestMathUtil.f90 +++ b/autotest/TestMathUtil.f90 @@ -5,7 +5,8 @@ module TestMathUtil to_string, unittest_type use MathUtilModule, only: f1d, is_close, mod_offset, & zero_ch, zero_br, & - get_perturbation, linspace + get_perturbation, & + arange, linspace implicit none private public :: collect_mathutil @@ -26,6 +27,7 @@ subroutine collect_mathutil(testsuite) test_zero_br), & new_unittest("get_perturbation", & test_get_perturbation), & + new_unittest("arange", test_arange), & new_unittest("linspace", test_linspace) & ] end subroutine collect_mathutil @@ -190,14 +192,17 @@ subroutine test_zero_ch(error) z = zero_ch(-1.0_DP, 1.0_DP, f, 0.001_DP) call check(error, is_close(z, 0.0_DP, atol=1d-6), & 'expected 0, got: '//to_string(z)) + if (allocated(error)) return z = zero_ch(-4.0_DP, -1.0_DP, f, 0.001_DP) call check(error, is_close(z, -pi, atol=1d-6), & 'expected -pi, got: '//to_string(z)) + if (allocated(error)) return z = zero_ch(1.0_DP, 4.0_DP, f, 0.001_DP) call check(error, is_close(z, pi, atol=1d-6), & 'expected pi, got: '//to_string(z)) + if (allocated(error)) return end subroutine test_zero_ch subroutine test_zero_br(error) @@ -211,14 +216,17 @@ subroutine test_zero_br(error) z = zero_br(-1.0_DP, 1.0_DP, f, 0.001_DP) call check(error, is_close(z, 0.0_DP, atol=1d-6), & 'expected 0, got: '//to_string(z)) + if (allocated(error)) return z = zero_br(-4.0_DP, -1.0_DP, f, 0.001_DP) call check(error, is_close(z, -pi, atol=1d-6), & 'expected -pi, got: '//to_string(z)) + if (allocated(error)) return z = zero_br(1.0_DP, 4.0_DP, f, 0.001_DP) call check(error, is_close(z, pi, atol=1d-6), & 'expected pi, got: '//to_string(z)) + if (allocated(error)) return end subroutine test_zero_br subroutine test_get_perturbation(error) @@ -232,6 +240,7 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-12), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return ! test derivative calculation for sin(x) where x=1 x = 1.d0 @@ -241,6 +250,7 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-5), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return ! test derivative calculation for sin(x) where x=0 x = 0.d0 @@ -250,6 +260,7 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-5), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return ! test derivative calculation for x ** 2 x = 1.d6 @@ -259,9 +270,35 @@ subroutine test_get_perturbation(error) call check(error, & is_close(v1, v2, atol=1d-1), & 'expected '//to_string(v1)//' got: '//to_string(v2)) + if (allocated(error)) return end subroutine test_get_perturbation + subroutine test_arange(error) + type(error_type), allocatable, intent(out) :: error + real(DP), allocatable :: a(:) + integer(I4B) :: i + + a = arange(DZERO, 10.0_DP, DONE) + call check(error, size(a) == 10, "wrong size: "//to_string(size(a))) + if (allocated(error)) return + do i = 1, 10 + call check(error, is_close(a(i), real(i - 1, DP))) + if (allocated(error)) return + end do + + deallocate (a) + + a = arange(DZERO, DONE, 0.1_DP) + call check(error, size(a) == 10, "wrong size: "//to_string(size(a))) + if (allocated(error)) return + do i = 1, 10 + call check(error, is_close(a(i), real(i - 1, DP) / 10.0_DP)) + if (allocated(error)) return + end do + + end subroutine test_arange + subroutine test_linspace(error) type(error_type), allocatable, intent(out) :: error real(DP), allocatable :: a(:) @@ -269,8 +306,10 @@ subroutine test_linspace(error) a = linspace(DONE, 10.0_DP, 10) call check(error, size(a) == 10) + if (allocated(error)) return do i = 1, 10 call check(error, is_close(a(i), real(i, DP))) + if (allocated(error)) return end do end subroutine test_linspace diff --git a/autotest/__snapshots__/test_prt_release_timing/test_mf6model[prtreltfreq].npy b/autotest/__snapshots__/test_prt_release_timing/test_mf6model[prtreltfreq].npy index 1bb7d8859c5..ffa70cc253c 100644 Binary files a/autotest/__snapshots__/test_prt_release_timing/test_mf6model[prtreltfreq].npy and b/autotest/__snapshots__/test_prt_release_timing/test_mf6model[prtreltfreq].npy differ diff --git a/src/Model/ParticleTracking/prt-prp.f90 b/src/Model/ParticleTracking/prt-prp.f90 index 58f02b252d2..ace431f491d 100644 --- a/src/Model/ParticleTracking/prt-prp.f90 +++ b/src/Model/ParticleTracking/prt-prp.f90 @@ -25,7 +25,7 @@ module PrtPrpModule use DisModule, only: DisType use DisvModule, only: DisvType use ErrorUtilModule, only: pstop - use MathUtilModule, only: is_close, linspace + use MathUtilModule, only: arange, is_close use ArrayHandlersModule, only: ExpandArray implicit none @@ -1046,10 +1046,10 @@ subroutine prp_load_releasetimefrequency(this) if (this%rtfreq <= DZERO) return ! create array of regularly-spaced release times - times = linspace( & + times = arange( & start=DZERO, & stop=totalsimtime, & - num=ceiling(totalsimtime / this%rtfreq)) + step=this%rtfreq) ! register times with release schedule call this%schedule%time_select%extend(times) diff --git a/src/Utilities/MathUtil.f90 b/src/Utilities/MathUtil.f90 index b155552cea6..c0719175f38 100644 --- a/src/Utilities/MathUtil.f90 +++ b/src/Utilities/MathUtil.f90 @@ -8,7 +8,7 @@ module MathUtilModule implicit none private public :: f1d, is_close, mod_offset, zero_ch, zero_br, & - get_perturbation, linspace + get_perturbation, arange, linspace interface mod_offset module procedure :: mod_offset_int, mod_offset_dbl @@ -375,7 +375,38 @@ function get_perturbation(x) result(res) res = DPRECSQRT * max(abs(x), DONE) * sign(DONE, x) end function get_perturbation - !> @brief Generate evenly spaced reals over the specified interval. + !> @brief Return reals separated by the given step over the given interval. + !! + !! This function is designed to behave like NumPy's arange. + !! Adapted from https://stackoverflow.com/a/65839162/6514033. + !< + pure function arange(start, stop, step) result(array) + ! dummy + real(DP), intent(in) :: start, stop + real(DP), intent(in), optional :: step + real(DP), allocatable :: array(:) + ! local + real(DP) :: lstep + integer(I4B) :: i, n + + if (present(step)) then + lstep = step + else + lstep = DONE + end if + + if (step <= 0) then + allocate (array(0)) + else + n = ceiling((stop - start) / step) + allocate (array(n)) + do i = 0, n - 1 + array(i + 1) = start + step * i + end do + end if + end function arange + + !> @brief Return evenly spaced reals over the given interval. !! !! This function is designed to behave like NumPy's linspace. !! Adapted from https://stackoverflow.com/a/57211848/6514033.