From 0e5cbb9e03369cd8827d0a63fb804f4f86920b18 Mon Sep 17 00:00:00 2001 From: wpbonelli Date: Sun, 19 Jan 2025 18:28:37 -0500 Subject: [PATCH] feat(grb): add binary grid file reader, support grid file in fmi (WIP) --- autotest/test_prt_voronoi2.py | 2 + doc/mf6io/mf6ivar/dfn/prt-fmi.dfn | 2 +- make/makefile | 1 + msvs/mf6core.vfproj | 1 + .../ModelUtilities/FlowModelInterface.f90 | 74 ++--- src/Model/TransportModel/tsp-fmi.f90 | 1 + src/Utilities/GridFileReader.f90 | 272 ++++++++++++++++++ src/meson.build | 1 + 8 files changed, 309 insertions(+), 45 deletions(-) create mode 100644 src/Utilities/GridFileReader.f90 diff --git a/autotest/test_prt_voronoi2.py b/autotest/test_prt_voronoi2.py index 59ff282a2bb..d5004e8cfad 100644 --- a/autotest/test_prt_voronoi2.py +++ b/autotest/test_prt_voronoi2.py @@ -175,9 +175,11 @@ def build_prt_sim(name, gwf_ws, prt_ws, targets, cell_ids): ) gwf_budget_file = gwf_ws / f"{gwf_name}.bud" gwf_head_file = gwf_ws / f"{gwf_name}.hds" + grb_file = gwf_ws / f"{gwf_name}.disv.grb" flopy.mf6.ModflowPrtfmi( prt, packagedata=[ + ("GWFGRID", grb_file), ("GWFHEAD", gwf_head_file), ("GWFBUDGET", gwf_budget_file), ], diff --git a/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn b/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn index 3deb36646fd..61a99fe0746 100644 --- a/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn +++ b/doc/mf6io/mf6ivar/dfn/prt-fmi.dfn @@ -25,7 +25,7 @@ type string tagged false reader urword longname flow type -description is the word GWFBUDGET or GWFHEAD. If GWFBUDGET is specified, then the corresponding file must be a budget file from a previous GWF Model run. +description is the word GWFBUDGET, GWFHEAD, or GRB. If GWFBUDGET is specified, then the corresponding file must be a budget file from a previous GWF Model run. If GWFHead is specified, the file must be a head file. If GWFGRID is specified, the file must be a binary grid file. block packagedata name filein diff --git a/make/makefile b/make/makefile index 27d99da2582..267a1aedef3 100644 --- a/make/makefile +++ b/make/makefile @@ -227,6 +227,7 @@ $(OBJDIR)/SourceCommon.o \ $(OBJDIR)/MemoryManagerExt.o \ $(OBJDIR)/ats.o \ $(OBJDIR)/GeomUtil.o \ +$(OBJDIR)/GridFileReader.o \ $(OBJDIR)/TimeSeriesLink.o \ $(OBJDIR)/TimeSeriesFileList.o \ $(OBJDIR)/tdis.o \ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index a90fddd986e..ae7b472b7fc 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -624,6 +624,7 @@ + diff --git a/src/Model/ModelUtilities/FlowModelInterface.f90 b/src/Model/ModelUtilities/FlowModelInterface.f90 index 60a4729e257..18db7c67e49 100644 --- a/src/Model/ModelUtilities/FlowModelInterface.f90 +++ b/src/Model/ModelUtilities/FlowModelInterface.f90 @@ -10,6 +10,7 @@ module FlowModelInterfaceModule use ListModule, only: ListType use BudgetFileReaderModule, only: BudgetFileReaderType use HeadFileReaderModule, only: HeadFileReaderType + use GridFileReaderModule, only: load_grb use PackageBudgetModule, only: PackageBudgetType use BudgetObjectModule, only: BudgetObjectType, budgetobject_cr_bfr @@ -24,6 +25,7 @@ module FlowModelInterfaceModule type(ListType), pointer :: gwfbndlist => null() !< list of gwf stress packages integer(I4B), pointer :: iflowsupdated => null() !< flows were updated for this time step integer(I4B), dimension(:), pointer, contiguous :: ibound => null() !< pointer to Model ibound + class(DisBaseType), pointer :: gwfdis => null() !< flow model discretization object real(DP), dimension(:), pointer, contiguous :: gwfflowja => null() !< pointer to the GWF flowja array real(DP), dimension(:, :), pointer, contiguous :: gwfspdis => null() !< pointer to npf specific discharge array real(DP), dimension(:), pointer, contiguous :: gwfhead => null() !< pointer to the GWF head array @@ -37,12 +39,13 @@ module FlowModelInterfaceModule integer(I4B), pointer :: iubud => null() !< unit number GWF budget file integer(I4B), pointer :: iuhds => null() !< unit number GWF head file integer(I4B), pointer :: iumvr => null() !< unit number GWF mover budget file + integer(I4B), pointer :: iugrb => null() !< unit number binary grid file integer(I4B), pointer :: nflowpack => null() !< number of GWF flow packages integer(I4B), dimension(:), pointer, contiguous :: igwfmvrterm => null() !< flag to indicate that gwf package is a mover term type(BudgetFileReaderType) :: bfr !< budget file reader type(HeadFileReaderType) :: hfr !< head file reader type(PackageBudgetType), dimension(:), allocatable :: gwfpackages !< used to get flows between a package and gwf - type(BudgetObjectType), pointer :: mvrbudobj => null() !< pointer to the mover budget budget object + type(BudgetObjectType), pointer :: mvrbudobj => null() !< pointer to the mover budget object character(len=16), dimension(:), allocatable :: flowpacknamearray !< array of boundary package names (e.g. LAK-1, SFR-3, etc.) character(len=LENVARNAME) :: depvartype = '' @@ -66,7 +69,6 @@ module FlowModelInterfaceModule procedure :: initialize_hfr procedure :: read_options procedure :: read_packagedata - end type FlowModelInterfaceType contains @@ -110,6 +112,7 @@ subroutine fmi_df(this, dis, idryinactive) if (this%inunit /= 0) then call this%read_options() end if + ! ! -- Read packagedata options if (this%inunit /= 0 .and. this%flows_from_file) then @@ -151,6 +154,7 @@ subroutine fmi_da(this) use MemoryManagerModule, only: mem_deallocate ! -- dummy class(FlowModelInterfaceType) :: this + ! -- todo: finalize hfr and bfr either here or in a finalize routine ! ! -- deallocate any memory stored with gwfpackages @@ -181,6 +185,7 @@ subroutine fmi_da(this) call mem_deallocate(this%iubud) call mem_deallocate(this%iuhds) call mem_deallocate(this%iumvr) + call mem_deallocate(this%iugrb) call mem_deallocate(this%nflowpack) call mem_deallocate(this%idryinactive) ! @@ -208,6 +213,7 @@ subroutine allocate_scalars(this) call mem_allocate(this%iubud, 'IUBUD', this%memoryPath) call mem_allocate(this%iuhds, 'IUHDS', this%memoryPath) call mem_allocate(this%iumvr, 'IUMVR', this%memoryPath) + call mem_allocate(this%iugrb, 'IUGRB', this%memoryPath) call mem_allocate(this%nflowpack, 'NFLOWPACK', this%memoryPath) call mem_allocate(this%idryinactive, "IDRYINACTIVE", this%memoryPath) ! @@ -220,6 +226,7 @@ subroutine allocate_scalars(this) this%iubud = 0 this%iuhds = 0 this%iumvr = 0 + this%iugrb = 0 this%nflowpack = 0 this%idryinactive = 1 end subroutine allocate_scalars @@ -346,12 +353,12 @@ subroutine read_packagedata(this) iapt = 0 blockrequired = .true. ! - ! -- get options block + ! -- get packagedata block call this%parser%GetBlock('PACKAGEDATA', isfound, ierr, & blockRequired=blockRequired, & supportOpenClose=.true.) ! - ! -- parse options block if detected + ! -- parse packagedata block if detected if (isfound) then write (this%iout, '(1x,a)') 'PROCESSING FMI PACKAGEDATA' do @@ -410,6 +417,19 @@ subroutine read_packagedata(this) call budgetobject_cr_bfr(this%mvrbudobj, 'MVT', this%iumvr, & this%iout) call this%mvrbudobj%fill_from_bfr(this%dis, this%iout) + case ('GWFGRID') + call this%parser%GetStringCaps(keyword) + if (keyword /= 'FILEIN') then + call store_error('GWFGRID KEYWORD MUST BE FOLLOWED BY '// & + '"FILEIN" then by filename.') + call this%parser%StoreErrorUnit() + end if + call this%parser%GetString(fname) + inunit = getunit() + call openfile(inunit, this%iout, fname, 'DATA(BINARY)', FORM, & + ACCESS, 'UNKNOWN') + this%iugrb = inunit + this%gwfdis => load_grb(this%iugrb, this%iout) case default write (errmsg, '(a,3(1x,a))') & 'UNKNOWN', trim(adjustl(this%text)), 'PACKAGEDATA:', trim(keyword) @@ -421,25 +441,18 @@ subroutine read_packagedata(this) end subroutine read_packagedata !> @brief Initialize the budget file reader - !< subroutine initialize_bfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy integer(I4B) :: ncrbud - ! - ! -- Initialize the budget file reader call this%bfr%initialize(this%iubud, this%iout, ncrbud) - ! - ! -- todo: need to run through the budget terms - ! and do some checking + ! todo: need to run through the budget terms + ! and do some checking end subroutine initialize_bfr !> @brief Advance the budget file reader !! !! Advance the budget file reader by reading the next chunk !! of information for the current time step and stress period. - !! !< subroutine advance_bfr(this) ! -- modules @@ -586,35 +599,22 @@ subroutine advance_bfr(this) end subroutine advance_bfr !> @brief Finalize the budget file reader - !< subroutine finalize_bfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy - ! - ! -- Finalize the budget file reader call this%bfr%finalize() - ! end subroutine finalize_bfr !> @brief Initialize the head file reader - !< subroutine initialize_hfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy - ! - ! -- Initialize the budget file reader call this%hfr%initialize(this%iuhds, this%iout) - ! - ! -- todo: need to run through the head terms - ! and do some checking + ! todo: need to run through the head terms + ! and do some checking end subroutine initialize_hfr !> @brief Advance the head file reader - !< subroutine advance_hfr(this) - ! -- modules + ! modules use TdisModule, only: kstp, kper class(FlowModelInterfaceType) :: this integer(I4B) :: nu, nr, i, ilay @@ -705,22 +705,15 @@ subroutine advance_hfr(this) end subroutine advance_hfr !> @brief Finalize the head file reader - !< subroutine finalize_hfr(this) - ! -- modules class(FlowModelInterfaceType) :: this - ! -- dummy - ! - ! -- Finalize the head file reader close (this%iuhds) - ! end subroutine finalize_hfr !> @brief Initialize gwf terms from budget file !! !! initialize terms and figure out how many !! different terms and packages are contained within the file - !! !< subroutine initialize_gwfterms_from_bfr(this) ! -- modules @@ -818,7 +811,6 @@ subroutine initialize_gwfterms_from_bfr(this) end subroutine initialize_gwfterms_from_bfr !> @brief Initialize gwf terms from a GWF exchange - !< subroutine initialize_gwfterms_from_gwfbndlist(this) ! -- modules use BndModule, only: BndType, GetBndFromList @@ -918,22 +910,16 @@ subroutine allocate_gwfpackages(this, ngwfterms) end subroutine allocate_gwfpackages !> @brief Deallocate memory in the gwfpackages array - !< subroutine deallocate_gwfpackages(this) - ! -- modules - ! -- dummy class(FlowModelInterfaceType) :: this - ! -- local integer(I4B) :: n - ! - ! -- initialize + do n = 1, this%nflowpack call this%gwfpackages(n)%da() end do end subroutine deallocate_gwfpackages !> @brief Find the package index for the package with the given name - !< subroutine get_package_index(this, name, idx) use BndModule, only: BndType, GetBndFromList class(FlowModelInterfaceType) :: this diff --git a/src/Model/TransportModel/tsp-fmi.f90 b/src/Model/TransportModel/tsp-fmi.f90 index bd696f98fab..caed1abf753 100644 --- a/src/Model/TransportModel/tsp-fmi.f90 +++ b/src/Model/TransportModel/tsp-fmi.f90 @@ -351,6 +351,7 @@ subroutine gwtfmi_da(this) call mem_deallocate(this%iubud) call mem_deallocate(this%iuhds) call mem_deallocate(this%iumvr) + call mem_deallocate(this%iugrb) call mem_deallocate(this%nflowpack) call mem_deallocate(this%idryinactive) ! diff --git a/src/Utilities/GridFileReader.f90 b/src/Utilities/GridFileReader.f90 new file mode 100644 index 00000000000..d82f408b76b --- /dev/null +++ b/src/Utilities/GridFileReader.f90 @@ -0,0 +1,272 @@ +module GridFileReaderModule + + use KindModule + use SimModule, only: store_error, store_error_unit + use ConstantsModule, only: LINELENGTH + use BaseDisModule, only: DisBaseType + use DisModule, only: DisType + use DisvModule, only: DisvType + use DisuModule, only: DisuType + use InputOutputModule, only: urword, read_line, upcase + use LongLineReaderModule, only: LongLineReaderType + + implicit none + + private :: load_dis, load_disv, load_disu + public :: load_grb + +contains + + !> @brief Load a discretization from a binary grid file. + function load_grb(inunit, iout) result(dis) + ! dummy + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in), optional :: iout + class(DisBaseType), pointer :: dis + ! local + type(LongLineReaderType) :: line_reader + character(len=:), allocatable :: line + integer(I4B) :: liout, lloc, istart, istop, ival, ierr + real(DP) :: rval + character(len=10) :: grid + integer(I4B) :: version, ntxt, lentxt + + if (present(iout)) then + liout = iout + else + liout = 0 + end if + + ! grid type + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + grid = line(istart:istop) + call upcase(grid) + + ! version + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + version = ival + + ! ntxt + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + ntxt = ival + + ! lentxt + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + lentxt = ival + + ! load grid-specific data + select case (grid) + case ('DIS') + dis => load_dis(inunit, liout) + case ('DISV') + dis => load_disv(inunit, liout) + case ('DISU') + dis => load_disu(inunit, liout) + end select + + ! close file + close (inunit) + end function load_grb + + function load_dis(inunit, iout) result(dis) + ! dummy + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in), optional :: iout + type(DisType), pointer :: dis + ! local + type(LongLineReaderType) :: line_reader + integer(I4B) :: liout + + if (present(iout)) then + liout = iout + else + liout = 0 + end if + + ! 1. TODO read header text + + end function load_dis + + function load_disv(inunit, iout) result(disv) + ! dummy + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in), optional :: iout + type(DisvType), pointer :: disv + ! local + type(LongLineReaderType) :: line_reader + character(len=:), allocatable :: line + integer(I4B) :: liout, dim, lloc, istart, istop, ival, ierr + real(DP) :: rval, xorigin, yorigin, angrot + integer(I4B) :: ncells, nlay, ncpl, nvert, njavert, nja + integer(I4B), allocatable :: shp(:) + + if (present(iout)) then + liout = iout + else + liout = 0 + end if + + ! 1. read header text + + ! ncells integer ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + ncells = ival + + ! nlay integer ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + nlay = ival + + ! ncpl integer ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + ncpl = ival + + ! nvert integer ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + nvert = ival + + ! njavert integer ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + njavert = ival + + ! nja integer ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 2, ival, rval, liout, 0) + nja = ival + + ! xorigin double ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 3, ival, rval, liout, 0) + xorigin = rval + + ! yorigin double ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 3, ival, rval, liout, 0) + yorigin = rval + + ! angrot double ndim 0 # + call line_reader%rdcom(inunit, liout, line, ierr) + lloc = 1 + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 1, ival, rval, liout, 0) + call urword(line, lloc, istart, istop, 3, ival, rval, liout, 0) + angrot = rval + + ! top + + ! botm + + ! vertices + + ! cellx + + ! celly + + ! iavert + + ! javert + + ! ia + + ! ja + + ! idomain + + ! icelltype + + ! 2. allocate + allocate (disv) + ! TODO allocate under fmi mempath? "floating" dis + ! owned by fmi, rather than by a model as normal? + + ! 3. TODO read/populate arrays + + end function load_disv + + function load_disu(inunit, iout) result(disu) + ! dummy + integer(I4B), intent(in) :: inunit + integer(I4B), intent(in), optional :: iout + type(DisuType), pointer :: disu + ! local + type(LongLineReaderType) :: line_reader + integer(I4B) :: liout + + if (present(iout)) then + liout = iout + else + liout = 0 + end if + + ! 1. TODO read header text + + end function load_disu + +end module GridFileReaderModule diff --git a/src/meson.build b/src/meson.build index 0cd60f953c5..7a15500bd7d 100644 --- a/src/meson.build +++ b/src/meson.build @@ -383,6 +383,7 @@ modflow_sources = files( 'Utilities' / 'DevFeature.f90', 'Utilities' / 'ErrorUtil.f90', 'Utilities' / 'GeomUtil.f90', + 'Utilities' / 'GridFileReader.f90', 'Utilities' / 'HashTable.f90', 'Utilities' / 'HeadFileReader.f90', 'Utilities' / 'HGeoUtil.f90',