fix(prt): fix duplicate particles, tidier release mechanism (MODFLOW-…

…USGS#1928)

Introduce a module/type to manage the particle release schedule. This uses the time and time step selection (MODFLOW-USGS#1923) types recently added, allowing cleaner release logic in the PRP package. Make some quality of life improvements to the time selection type for easier usage. Add utilities arange and linspace to the math utility module, modeled after the eponymous NumPy functions. Use these in the PRP package to support regularly spaced release times via a new RELEASE_TIME_FREQUENCY option.

This includes a few breaking changes, which address correctness and/or do not involve a loss of functionality.

1. Remove the FRACTION period-block release setting option. The RELEASETIMES block (MODFLOW-USGS#1989) can achieve the same explicitly. This brings period-block settings to equivalence with (and shares the implementation of) OC settings:

    - ALL
    - FIRST
    - LAST
    - STEPS
    - FREQUENCY

If FRACTION is provided, the simulation will now fail. Add a removal notice to the DFN which will get picked up at release time by the automation and an entry will be inserted into the deprecations/removals table in the release notes.

2. Clarify the conceptual model for release points and fix some potential issues with the release mechanism. Previously a release was made for the union of explicit release times and period-block settings, even if an explicitly provided release time happened to coincide with a period block release setting's time. This could introduce duplicate particles, which is a bug given that we describe particles as uniquely identifiable via composite key consisting of model ID, PRP ID, release point ID, and release time. Only allow one release from each point at a time — by default within machine tolerance * 10^9, configurable via new RELEASE_TIME_TOLERANCE option. The release schedule is now the union of times explicitly provided in the RELEASETIMES block, regularly spaced times configured via RELEASE_TIME_FREQUENCY, and period-block time settings, where times within the configured tolerance of one another are merged into a single release time.

autotest/TestMathUtil.f90

@@ -1,11 +1,12 @@
 module TestMathUtil
   use KindModule, only: I4B, DP
-  use ConstantsModule, only: DNODATA, DZERO
+  use ConstantsModule, only: DNODATA, DZERO, DONE
   use testdrive, only: check, error_type, new_unittest, test_failed, &
                        to_string, unittest_type
   use MathUtilModule, only: f1d, is_close, mod_offset, &
                             zero_ch, zero_br, &
-                            get_perturbation
+                            get_perturbation, &
+                            arange, linspace
   implicit none
   private
   public :: collect_mathutil
@@ -25,7 +26,9 @@ subroutine collect_mathutil(testsuite)
                 new_unittest("zero_br", &
                              test_zero_br), &
                 new_unittest("get_perturbation", &
-                             test_get_perturbation) &
+                             test_get_perturbation), &
+                new_unittest("arange", test_arange), &
+                new_unittest("linspace", test_linspace) &
                 ]
   end subroutine collect_mathutil
 
@@ -189,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)
@@ -210,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)
@@ -231,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
@@ -240,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
@@ -249,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
@@ -258,7 +270,48 @@ 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
+
 end module TestMathUtil

autotest/TestTimeSelect.f90

@@ -13,21 +13,24 @@ subroutine collect_timeselect(testsuite)
     type(unittest_type), allocatable, intent(out) :: testsuite(:)
     testsuite = [ &
                 new_unittest("is_increasing", test_is_increasing), &
-                new_unittest("slice", test_slice) &
+                new_unittest("select", test_select), &
+                new_unittest("extend_and_sort", &
+                             test_extend_and_sort) &
                 ]
   end subroutine collect_timeselect
 
   subroutine test_is_increasing(error)
     type(error_type), allocatable, intent(out) :: error
     type(TimeSelectType) :: ts
 
+    call ts%init()
     call ts%expand(3)
 
     ! increasing
     ts%times = (/0.0_DP, 1.0_DP, 2.0_DP/)
     call check(error, ts%increasing())
 
-    ! not decreasing
+    ! not decreasing (duplicates)
     ts%times = (/0.0_DP, 0.0_DP, 2.0_DP/)
     call check(error,.not. ts%increasing())
 
@@ -36,11 +39,12 @@ subroutine test_is_increasing(error)
     call check(error,.not. ts%increasing())
   end subroutine
 
-  subroutine test_slice(error)
+  subroutine test_select(error)
     type(error_type), allocatable, intent(out) :: error
     type(TimeSelectType) :: ts
     logical(LGP) :: changed
 
+    call ts%init()
     call ts%expand(3)
     ts%times = (/0.0_DP, 1.0_DP, 2.0_DP/)
     call check( &
@@ -107,5 +111,19 @@ subroutine test_slice(error)
       "lb ub eq slice failed, got [" &
       //to_string(ts%selection(1))//","//to_string(ts%selection(2))//"]")
 
-  end subroutine test_slice
+  end subroutine test_select
+
+  subroutine test_extend_and_sort(error)
+    type(error_type), allocatable, intent(out) :: error
+    type(TimeSelectType) :: ts
+    real(DP) :: a(3)
+
+    a = (/0.0_DP, 2.0_DP, 1.0_DP/)
+
+    call ts%init()
+    call ts%extend(a)
+    call check(error, ts%increasing())
+
+  end subroutine test_extend_and_sort
+
 end module TestTimeSelect

autotest/test_prt_disv1.py

@@ -456,10 +456,9 @@ def check_output(idx, test):
     assert list_file.is_file()
     lines = open(list_file).readlines()
     lines = [l.strip() for l in lines]
-    assert (
-        "PARTICLE RELEASE:      TIME STEP(S) 1  AT OFFSET           0.000"
-        in lines
-    )
+    for iprp in range(1, 3):
+        i = lines.index(f"PARTICLE RELEASE FOR PRP {iprp}")
+        assert "RELEASE SCHEDULE:" in lines[i + 1]
 
     # load mp7 pathline results
     plf = PathlineFile(mp7_ws / mp7_pathline_file)