From a8347e673246a7f2ccbe70aa02c93aa38d9970c1 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Thu, 16 Jan 2025 11:16:25 -0600 Subject: [PATCH] feat(disu-idomain): write idomain to binary grid file for disu models (#2144) The binary grid file written by MODFLOW 6 for DISU models did not include the IDOMAIN array. The binary grid file now includes IDOMAIN for all discretization types, including DISU. --- autotest/test_gwf_disu.py | 7 ++++++- doc/ReleaseNotes/develop.tex | 2 +- doc/mf6io/framework/binaryoutput.tex | 26 ++++++++++++++------------ src/Model/Discretization/Disu.f90 | 6 +++++- 4 files changed, 26 insertions(+), 15 deletions(-) diff --git a/autotest/test_gwf_disu.py b/autotest/test_gwf_disu.py index 3bb34f7217d..242d8fd65cd 100644 --- a/autotest/test_gwf_disu.py +++ b/autotest/test_gwf_disu.py @@ -47,7 +47,7 @@ def build_models(idx, test): ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0) npf = flopy.mf6.ModflowGwfnpf(gwf) spd = {0: [[(0,), 1.0], [(nrow * ncol - 1), 0.0]]} - chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, stress_period_data=spd) + chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=spd) return sim, None @@ -59,6 +59,10 @@ def check_output(idx, test): nodes = grbobj._datadict["NODES"] ia = grbobj._datadict["IA"] ja = grbobj._datadict["JA"] + idomain = grbobj._datadict["IDOMAIN"] + + if idx == 0: + assert np.array_equal(idomain, np.array(27 * [1]), int) if idx == 1: assert np.array_equal(ia[0:4], np.array([1, 4, 4, 7])) @@ -66,6 +70,7 @@ def check_output(idx, test): assert ia[-1] == 127 assert ia.shape[0] == 28, "ia should have size of 28" assert ja.shape[0] == 126, "ja should have size of 126" + assert np.array_equal(idomain, np.array([1, 0] + 25 * [1]), int) @pytest.mark.parametrize("idx, name", enumerate(cases)) diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 99c7dec2323..18a428a612b 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -6,7 +6,7 @@ %\underline{NEW FUNCTIONALITY} %\begin{itemize} \item Support for adjusting time step lengths using the adaptive time stepping (ATS) capability was added to the GWT Advection (ADV) Package of the Groundwater Transport (GWT) Model in release 6.6.0. The same functionality that was added to GWT is now available with the Groundwater Energy Transport (GWE) Model. A description of how this functionality works and how to activate it can be found in the release notes for version 6.6.0 (Appendix A) and in the MODFLOW 6 input-output guide. -% \item xxx + \item The binary grid file written by MODFLOW 6 for DISU models did not include the IDOMAIN array. The binary grid file now includes IDOMAIN for all discretization types, including DISU. % \item xxx %\end{itemize} diff --git a/doc/mf6io/framework/binaryoutput.tex b/doc/mf6io/framework/binaryoutput.tex index f8c2282c672..c89b9dc0676 100644 --- a/doc/mf6io/framework/binaryoutput.tex +++ b/doc/mf6io/framework/binaryoutput.tex @@ -157,15 +157,16 @@ \subsubsection{DISU Grids} \noindent Definition 7: \texttt{`BOT DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ \noindent Definition 8: \texttt{`IA INTEGER NDIM 1 ncells+1'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ \noindent Definition 9: \texttt{`JA INTEGER NDIM 1 nja'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ -\noindent Definition 10: \texttt{`ICELLTYPE INTEGER NDIM 1 ncells'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 10: \texttt{`IDOMAIN INTEGER NDIM 1 ncells'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 11: \texttt{`ICELLTYPE INTEGER NDIM 1 ncells'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ \vspace{5mm} \noindent If vertices are provided in the DISU Package, then 5 additional definitions are included: \\ -\noindent Definition 11: \texttt{`VERTICES DOUBLE NDIM 2 2 nvert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ -\noindent Definition 12: \texttt{`CELLX DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ -\noindent Definition 13: \texttt{`CELLY DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ -\noindent Definition 14: \texttt{`IAVERT INTEGER NDIM 1 nodes+1'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ -\noindent Definition 15: \texttt{`JAVERT INTEGER NDIM 1 njavert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 12: \texttt{`VERTICES DOUBLE NDIM 2 2 nvert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 13: \texttt{`CELLX DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 14: \texttt{`CELLY DOUBLE NDIM 1 nodes'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 15: \texttt{`IAVERT INTEGER NDIM 1 nodes+1'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ +\noindent Definition 16: \texttt{`JAVERT INTEGER NDIM 1 njavert'} {\color{red} \footnotesize{CHARACTER(LEN=LENTXT)}} \\ \vspace{5mm} \noindent Read \texttt{NDAT} data variables using the definitions defined above. \\ @@ -178,15 +179,16 @@ \subsubsection{DISU Grids} \noindent Record 7: \texttt{((BOT(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}} \\ \noindent Record 8: \texttt{(IA(J),J=1,NODES+1)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NODES+1)}} \\ \noindent Record 9: \texttt{(JA(J),J=1,NJA)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NJA)}} \\ -\noindent Record 10: \texttt{(ICELLTYPE(J),J=1,NCELLS)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NCELLS)}} \\ +\noindent Record 10: \texttt{(IDOMAIN(J),J=1,NCELLS)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NCELLS)}} \\ +\noindent Record 11: \texttt{(ICELLTYPE(J),J=1,NCELLS)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NCELLS)}} \\ \vspace{5mm} \noindent If vertices are provided in the DISU Package, then 5 additional records are included: \\ -\noindent Record 11: \texttt{((VERT(J,K),J=1,2),K=1,NVERT)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(2,NVERT)}} \\ -\noindent Record 12: \texttt{(CELLX(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}}\\ -\noindent Record 13: \texttt{(CELLY(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}} \\ -\noindent Record 14: \texttt{(IAVERT(J),J=1,NODES+1)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NODES+1)}} \\ -\noindent Record 15: \texttt{(JAVERT(J),J=1,NJAVERT)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NJAVERT)}} \\ +\noindent Record 12: \texttt{((VERT(J,K),J=1,2),K=1,NVERT)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(2,NVERT)}} \\ +\noindent Record 13: \texttt{(CELLX(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}}\\ +\noindent Record 14: \texttt{(CELLY(J),J=1,NODES)} {\color{red} \footnotesize{DOUBLE PRECISION ARRAY SIZE(NODES)}} \\ +\noindent Record 15: \texttt{(IAVERT(J),J=1,NODES+1)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NODES+1)}} \\ +\noindent Record 16: \texttt{(JAVERT(J),J=1,NJAVERT)} {\color{red} \footnotesize{INTEGER ARRAY SIZE(NJAVERT)}} \\ \newpage diff --git a/src/Model/Discretization/Disu.f90 b/src/Model/Discretization/Disu.f90 index 7b07e24a3b3..15b2df64058 100644 --- a/src/Model/Discretization/Disu.f90 +++ b/src/Model/Discretization/Disu.f90 @@ -923,7 +923,7 @@ subroutine write_grb(this, icelltype) &/,6X,'UNIT NUMBER: ', I0,/,6X, 'FILE NAME: ', A)" ! ! -- Initialize - ntxt = 10 + ntxt = 11 if (this%nvert > 0) ntxt = ntxt + 5 ! ! -- Open the file @@ -975,6 +975,9 @@ subroutine write_grb(this, icelltype) write (txt, '(3a, i0)') 'JA ', 'INTEGER ', 'NDIM 1 ', this%con%nja txt(lentxt:lentxt) = new_line('a') write (iunit) txt + write (txt, '(3a, i0)') 'IDOMAIN ', 'INTEGER ', 'NDIM 1 ', this%nodesuser + txt(lentxt:lentxt) = new_line('a') + write (iunit) txt write (txt, '(3a, i0)') 'ICELLTYPE ', 'INTEGER ', 'NDIM 1 ', this%nodesuser txt(lentxt:lentxt) = new_line('a') write (iunit) txt @@ -1008,6 +1011,7 @@ subroutine write_grb(this, icelltype) write (iunit) this%bot1d ! bot write (iunit) this%con%iausr ! ia write (iunit) this%con%jausr ! ja + write (iunit) this%idomain ! idomain write (iunit) icelltype ! icelltype ! ! -- if vertices have been read then write additional data