fix time frequency impl (add/use arange instead of linspace)

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,18 +270,46 @@ 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(:)
     integer(I4B) :: i
 
     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

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)

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.