Skip to content

Commit

Permalink
refactor prp:
Browse files Browse the repository at this point in the history 
* simplify implementation
* only support scalar time offset (fraction)
* rename referencetime to releasetime
* check particle release location elevation
* modify array expansion to respect lower bound
  • Loading branch information
@wpbonelli
wpbonelli committed Dec 15, 2023
1 parent 6897472 commit c1067bf
Show file tree
Hide file tree
Showing 12 changed files with 555 additions and 709 deletions.
144 changes: 98 additions & 46 deletions autotest/TestArrayHandlers.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,11 @@ subroutine collect_arrayhandlers(testsuite)
new_unittest("ExpandArray_dbl", test_ExpandArray_dbl), &
new_unittest("ExpandArray_log", test_ExpandArray_log), &
new_unittest("ExpandArray2D_int", test_ExpandArray2D_int), &
new_unittest("ExpandArray2D_dbl", test_ExpandArray2D_dbl) &
new_unittest("ExpandArray2D_dbl", test_ExpandArray2D_dbl), &
new_unittest("ExpandArray2D_int_offset", &
test_ExpandArray2D_int_offset), &
new_unittest("ExpandArray2D_dbl_offset", &
test_ExpandArray2D_dbl_offset) &
]
end subroutine collect_arrayhandlers

Expand All @@ -29,19 +33,15 @@ subroutine test_ExpandArray_int(error)
array(1) = 0
array(2) = 1

! resize array
! resize array and check new size
call ExpandArray(array, 3)

! check that array has been resized
call check(error, size(array, 1) == 5, "1d int array resize failed")
if (allocated(error)) return

! set new array elements
! set new array elements and check new/old contents
array(3) = 2
array(4) = 3
array(5) = 4

! check array contents
call check(error, &
array(1) == 0 .and. &
array(2) == 1 .and. &
Expand All @@ -50,9 +50,7 @@ subroutine test_ExpandArray_int(error)
array(5) == 4, &
"1d int array repopulation failed")
if (allocated(error)) return

deallocate (array)

end subroutine test_ExpandArray_int

subroutine test_ExpandArray_dbl(error)
Expand All @@ -64,26 +62,20 @@ subroutine test_ExpandArray_dbl(error)
array(1) = 0.5_DP
array(2) = 0.7_DP

! resize array
! resize array and check new size
call ExpandArray(array, 1)

! check that array has been resized
call check(error, size(array, 1) == 3, "1d dbl array resize failed")
if (allocated(error)) return

! set new array element
! set new array element and check new/old contents
array(3) = 0.1_DP

! check array contents
call check(error, &
array(1) == 0.5_DP .and. &
array(2) == 0.7_DP .and. &
array(3) == 0.1_DP, &
"1d dbl array repopulation failed")
if (allocated(error)) return

deallocate (array)

end subroutine test_ExpandArray_dbl

subroutine test_ExpandArray_log(error)
Expand All @@ -95,54 +87,43 @@ subroutine test_ExpandArray_log(error)
array(1) = .true.
array(2) = .false.

! resize array
! resize array and check new size
call ExpandArray(array, 1)

! check that array has been resized
call check(error, size(array, 1) == 3, "1d logical array resize failed")
if (allocated(error)) return

! set an element in the array
! set an element in the array and check new/old contents
array(3) = .true.

! check array contents
call check(error, &
array(1) .and. &
.not. array(2) .and. &
array(3), &
"1d logical array repopulation failed")
if (allocated(error)) return

deallocate (array)

end subroutine test_ExpandArray_log

!> @brief Test 2D int array expansion with default lower bound of 1
subroutine test_ExpandArray2D_int(error)
type(error_type), allocatable, intent(out) :: error
integer(I4B), allocatable :: array(:, :)

! allocate array
! allocate array and check initial size
allocate (array(2, 2))
array(1, :) = (/1, 2/)
array(2, :) = (/2, 3/)

! check initial array size
call check(error, size(array, 1) == 2 .and. size(array, 2) == 2)
if (allocated(error)) return

! resize array
! resize array and check new size
call ExpandArray2D(array, 1, 1)

! check that array has been resized
call check(error, &
size(array, 1) == 3 .and. size(array, 2) == 3, &
"2d int array resize failed")
if (allocated(error)) return

! add new array elements
! add new array elements and check new/old contents
array(3, :) = (/3, 4, 5/)

! check array contents
call check(error, &
array(1, 1) == 1 .and. &
array(1, 2) == 2 .and. &
Expand All @@ -156,37 +137,30 @@ subroutine test_ExpandArray2D_int(error)
array(3, 2) == 4 .and. &
array(3, 3) == 5, &
"2d int array repopulation failed")

deallocate (array)

end subroutine test_ExpandArray2D_int

!> @brief Test 2D dbl array expansion with default lower bound of 1
subroutine test_ExpandArray2D_dbl(error)
type(error_type), allocatable, intent(out) :: error
real(DP), allocatable :: array(:, :)

! allocate array
! allocate array and check initial size
allocate (array(2, 2))
array(1, :) = (/1.0_DP, 2.0_DP/)
array(2, :) = (/2.0_DP, 3.0_DP/)

! check initial array size
call check(error, size(array, 1) == 2 .and. size(array, 2) == 2)
if (allocated(error)) return

! resize array
! resize array and check new size
call ExpandArray2D(array, 1, 1)

! check that array has been resized
call check(error, &
size(array, 1) == 3 .and. size(array, 2) == 3, &
"2d dbl array resize failed")
if (allocated(error)) return

! set new array elements
! set new array elements and check new/old contents
array(3, :) = (/3.0_DP, 4.0_DP, 5.0_DP/)

! check array contents
call check(error, &
array(1, 1) == 1.0_DP .and. &
array(1, 2) == 2.0_DP .and. &
Expand All @@ -200,8 +174,86 @@ subroutine test_ExpandArray2D_dbl(error)
array(3, 2) == 4.0_DP .and. &
array(3, 3) == 5.0_DP, &
"2d dbl array repopulation failed")
deallocate (array)
end subroutine test_ExpandArray2D_dbl

!> @brief Test 2D int array expansion when lower bound /= 1
subroutine test_ExpandArray2D_int_offset(error)
type(error_type), allocatable, intent(out) :: error
integer(I4B), allocatable :: array(:, :)

! allocate array and check initial size
allocate (array(0:1, 0:1))
array(0, :) = (/1, 2/)
array(1, :) = (/2, 3/)
call check(error, size(array, 1) == 2 .and. size(array, 2) == 2)
if (allocated(error)) return

! resize array and check new size
call ExpandArray2D(array, 1, 1)
call check(error, &
size(array, 1) == 3 .and. size(array, 2) == 3, &
"2d int array resize failed")
call check(error, &
lbound(array, 1) == 0 .and. lbound(array, 2) == 0, &
"2d int array unexpected lower bound")
if (allocated(error)) return

! add new array elements and check new/old contents
array(2, :) = (/3, 4, 5/)
call check(error, &
array(0, 0) == 1 .and. &
array(0, 1) == 2 .and. &
! can't guarantee unassigned item value
! array(1, 3) == 0 .and. &
array(1, 0) == 2 .and. &
array(1, 1) == 3 .and. &
! can't guarantee unassigned item value
! array(2, 3) == 0 .and. &
array(2, 0) == 3 .and. &
array(2, 1) == 4 .and. &
array(2, 2) == 5, &
"2d int array repopulation failed")
deallocate (array)
end subroutine test_ExpandArray2D_int_offset

end subroutine test_ExpandArray2D_dbl
!> @brief Test 2D dbl array expansion when lower bound /= 1
subroutine test_ExpandArray2D_dbl_offset(error)
type(error_type), allocatable, intent(out) :: error
real(DP), allocatable :: array(:, :)

! allocate array and check initial size
allocate (array(0:1, 0:1))
array(0, :) = (/1.0_DP, 2.0_DP/)
array(1, :) = (/2.0_DP, 3.0_DP/)
call check(error, size(array, 1) == 2 .and. size(array, 2) == 2)
if (allocated(error)) return

! resize array and check new size
call ExpandArray2D(array, 1, 1)
call check(error, &
size(array, 1) == 3 .and. size(array, 2) == 3, &
"2d dbl array resize failed")
call check(error, &
lbound(array, 1) == 0 .and. lbound(array, 2) == 0, &
"2d dbl array unexpected lower bound")
if (allocated(error)) return

! set new array elements and check new/old contents
array(2, :) = (/3.0_DP, 4.0_DP, 5.0_DP/)
call check(error, &
array(0, 0) == 1.0_DP .and. &
array(0, 1) == 2.0_DP .and. &
! can't guarantee unassigned item value
! array(1, 3) == 0.0_DP .and. &
array(1, 0) == 2.0_DP .and. &
array(1, 1) == 3.0_DP .and. &
! can't guarantee unassigned item value
! array(2, 3) == 0.0_DP .and. &
array(2, 0) == 3.0_DP .and. &
array(2, 1) == 4.0_DP .and. &
array(2, 2) == 5.0_DP, &
"2d dbl array repopulation failed")
deallocate (array)
end subroutine test_ExpandArray2D_dbl_offset
end module TestArrayHandlers
30 changes: 18 additions & 12 deletions autotest/prt/test_prt_fmi05.py
View file
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
"""
Test cases exercising release timing, 1st via
package-level REFERENCETIME option, then with
package-level RELEASETIME option, & then with
period-block config STEPS 1 and FRACTION 0.5.
The model is setup to release halfway through
the first and only time step of the first and
only stress period, with duration 1 time unit,
so the same value of 0.5 can be used for both
REFERENCETIME and FRACTION.
RELEASETIME and FRACTION.
Period-block FRACTION should work with FIRST
and ALL, but flopy hangs with either option.
Expand All @@ -18,7 +18,7 @@
Particles are released from the top left cell.
Results are compared against a MODPATH 7 model.
Referencetime 0.5 could be configured, but mp7
Telease time 0.5 could be configured, but mp7
reports relative times, so there is no reason
& mp7 results are converted before comparison.
"""
Expand All @@ -37,13 +37,19 @@
from flopy.plot.plotutil import to_mp7_pathlines
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from prt_test_utils import (all_equal, check_budget_data, check_track_data,
get_gwf_sim, get_model_name, get_partdata)
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
get_gwf_sim,
get_model_name,
get_partdata,
)

simname = "prtfmi05"
ex = [
cases = [
# options block options
f"{simname}reft", # REFERENCETIME 0.5
f"{simname}relt", # RELEASETIME 0.5
# period block options
# f"{simname}all", # ALL FRACTION 0.5 # todo debug flopy hanging
# f"{simname}frst", # FIRST FRACTION 0.5 # todo debug flopy hanging
Expand All @@ -52,7 +58,7 @@


def get_perioddata(name, periods=1, fraction=None) -> Optional[dict]:
if "reft" in name:
if "relt" in name:
return None
opt = [
"FIRST"
Expand Down Expand Up @@ -115,8 +121,8 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
prp_track_file = f"{prtname}.prp.trk"
prp_track_csv_file = f"{prtname}.prp.trk.csv"
pdat = get_perioddata(prtname, fraction=fraction)
# fraction 0.5 equiv. to reftime 0.5 since 1 period 1 step with length 1
reft = fraction if "reft" in prtname else None
# fraction 0.5 equiv. to release time 0.5 since 1 period 1 step with length 1
trelease = fraction if "relt" in prtname else None
flopy.mf6.ModflowPrtprp(
prt,
pname="prp1",
Expand All @@ -126,7 +132,7 @@ def build_prt_sim(ctx, name, ws, mf6, fraction=None):
perioddata=pdat,
track_filerecord=[prp_track_file],
trackcsv_filerecord=[prp_track_csv_file],
referencetime=reft,
releasetime=trelease,
)

# create output control package
Expand Down Expand Up @@ -195,7 +201,7 @@ def build_mp7_sim(ctx, name, ws, mp7, gwf):
return mp


@pytest.mark.parametrize("name", ex)
@pytest.mark.parametrize("name", cases)
@pytest.mark.parametrize("fraction", [0.5])
def test_mf6model(name, function_tmpdir, targets, fraction):
# workspace
Expand Down
13 changes: 10 additions & 3 deletions autotest/prt/test_prt_fmi06.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,9 +19,14 @@
from flopy.utils import PathlineFile
from flopy.utils.binaryfile import HeadFile
from flopy.utils.gridutil import get_disv_kwargs
from prt_test_utils import (all_equal, check_budget_data, check_track_data,
get_partdata, has_default_boundnames,
plot_nodes_and_vertices)
from prt_test_utils import (
all_equal,
check_budget_data,
check_track_data,
get_partdata,
has_default_boundnames,
plot_nodes_and_vertices,
)

simname = "prtfmi06"
ex = [f"{simname}", f"{simname}bprp"]
Expand Down Expand Up @@ -442,5 +447,7 @@ def test_mf6model(idx, name, function_tmpdir, targets):
del mp7_pls["node"]

# compare mf6 / mp7 pathline data
# import pdb
# pdb.set_trace()
assert mf6_pls.shape == mp7_pls.shape
assert np.allclose(mf6_pls, mp7_pls, atol=1e-3)
Loading

0 comments on commit c1067bf

Please sign in to comment.