diff --git a/CITATION.cff b/CITATION.cff index 1595372ce9e..74d93133475 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -2,8 +2,8 @@ cff-version: 1.2.0 message: If you use this software, please cite the software itself. type: software title: MODFLOW 6 Modular Hydrologic Model -version: 6.4.3 -date-released: '2024-02-07' +version: 6.4.4 +date-released: '2024-02-13' doi: 10.5066/F76Q1VQV abstract: MODFLOW 6 is an object-oriented program and framework developed to provide a platform for supporting multiple models and multiple types of models within the diff --git a/README.md b/README.md index 5e561c379b6..2c8c5782f82 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,7 @@ This is the development repository for the USGS MODFLOW 6 Hydrologic Model. The official USGS distribution is available at [USGS Release Page](https://water.usgs.gov/ogw/modflow/MODFLOW.html). -### Version 6.4.3 +### Version 6.4.4 [![GitHub release](https://img.shields.io/github/release/MODFLOW-USGS/modflow6.svg)](https://github.com/MODFLOW-USGS/modflow6/releases/latest) [![MODFLOW 6 continuous integration](https://github.com/MODFLOW-USGS/modflow6/actions/workflows/ci.yml/badge.svg)](https://github.com/MODFLOW-USGS/modflow6/actions/workflows/ci.yml) diff --git a/autotest/get_exes.py b/autotest/get_exes.py index d5c789ceb11..2345ea4f59c 100644 --- a/autotest/get_exes.py +++ b/autotest/get_exes.py @@ -97,7 +97,9 @@ def rebuild(): def test_get_executables(downloaded_bin_path: Path): print(f"Installing MODFLOW-related executables to: {downloaded_bin_path}") downloaded_bin_path.mkdir(exist_ok=True, parents=True) - flopy.utils.get_modflow(str(downloaded_bin_path)) + # todo: remove release_id workaround when double-precision comparison issues fixed + # https://github.com/MODFLOW-USGS/modflow6/pull/1612 + flopy.utils.get_modflow(str(downloaded_bin_path), release_id="14.0") if __name__ == "__main__": diff --git a/code.json b/code.json index 0538ab48eb0..171ccdeb1a4 100644 --- a/code.json +++ b/code.json @@ -18,9 +18,9 @@ "email": "langevin@usgs.gov" }, "laborHours": -1, - "version": "6.4.3", + "version": "6.4.4", "date": { - "metadataLastUpdated": "2024-02-07" + "metadataLastUpdated": "2024-02-13" }, "organization": "U.S. Geological Survey", "permissions": { diff --git a/doc/ReleaseNotes/ReleaseNotes.tex b/doc/ReleaseNotes/ReleaseNotes.tex index 0ea54536e6c..39cef7608cb 100644 --- a/doc/ReleaseNotes/ReleaseNotes.tex +++ b/doc/ReleaseNotes/ReleaseNotes.tex @@ -177,6 +177,7 @@ \section{Release History} 6.4.1 & December 9, 2022 & \url{https://doi.org/10.5066/P9FL1JCC} \\ 6.4.2 & June 28, 2023 & \url{https://doi.org/10.5066/P9FL1JCC} \\ 6.4.3 & February 7, 2024 & \url{https://doi.org/10.5066/P9FL1JCC} \\ +6.4.4 & February 13, 2024 & \url{https://doi.org/10.5066/P9FL1JCC} \\ \hline \label{tab:releases} \end{tabular*} diff --git a/doc/ReleaseNotes/appendixA.tex b/doc/ReleaseNotes/appendixA.tex index 199d7151043..bcc041cb842 100644 --- a/doc/ReleaseNotes/appendixA.tex +++ b/doc/ReleaseNotes/appendixA.tex @@ -1,5 +1,6 @@ This appendix describes changes introduced into MODFLOW~6 in previous releases. These changes may substantially affect users. + \input{./previous/v6.4.3.tex} \input{./previous/v6.4.2.tex} \input{./previous/v6.4.1.tex} \input{./previous/v6.4.0.tex} diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 05253619fd4..c6e7876bdd9 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -2,50 +2,47 @@ % after a release has just been made. \item \currentmodflowversion - - \underline{NEW FUNCTIONALITY} - \begin{itemize} - \item The Input Data Processor (IDP), first released in version 6.4.2, is a general utility for reading user-provided input files. Package-specific routines for reading input files continue to be replaced by the IDP approach. For packages that use IDP for input, logging information is written to the simulation list file (mfsim.lst). Additional information on the IDP and the list of supported packages is contained in the MODFLOW 6 Description of Input and Output (mf6io.pdf) under a section titled ``Processing of Program Input.'' - \item The source code was refactored to support compilation of a parallel version of MODFLOW 6 based on the Message Passing Interface (MPI) and the Portable, Extensible Toolkit for Scientific Computation (PETSc) libraries. The parallel version of MODFLOW is considered preliminary. Limited testing of the parallel version has been performed on laptops, desktops, and supercomputers, but significant changes are expected in future releases. User support for the parallel version of MODFLOW 6 may be provided in the future. + + %\underline{NEW FUNCTIONALITY} + %\begin{itemize} % \item xxx - \end{itemize} + % \item xxx + % \item xxx + %\end{itemize} - \underline{EXAMPLES} - \begin{itemize} - \item A new exampled called ex-gwf-rad-disu was added. This new example uses a DISU grid to represent radial groundwater flow to a pumping well. - \item A new exampled called ex-gwf-curv-90 was added. This new example demonstrates use of a DISV grid to represent a curvilinear spatial discretization. For this example, the curvilinear grid is applied to one quarter of a radial groundwater flow system. - \item A new exampled called ex-gwf-curvilin was added. This new example uses a curvilinear grid, represented with the DISV Package, to simulate groundwater flow through a multi-region aquifer with bends in the domain boundaries. - \end{itemize} + %\underline{EXAMPLES} + %\begin{itemize} + % \item xxx + % \item xxx + % \item xxx + %\end{itemize} - \textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\ - \underline{BASIC FUNCTIONALITY} - \begin{itemize} - \item Improve error message if the size of data read from a binary array file is inconsistent with READARRAY control line and variable description keywords. - \item The area calculation for cells in the DISV package was inaccurate for some cases with very large cell vertex coordinates. The area calculation was improved by using transformed cell vertex coordinates prior to making the area calculation. - \item Auxiliary variables in RCH and EVT Array-Based input packages are now reset to zero when otherwise not specified in period input data and the auxiliary parameter is not controlled by a time-series. + %\textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\ + %\underline{BASIC FUNCTIONALITY} + %\begin{itemize} % \item xxx % \item xxx - \end{itemize} + % \item xxx + %\end{itemize} - \underline{INTERNAL FLOW PACKAGES} - \begin{itemize} - \item The data header in the binary output file written by the viscosity (VSC) package was printing `` VISCOSI'' instead of `` VISCOSITY''. The viscosity package now prints the full `` VISCOSITY'' header in the binary output file. - \item The CSUB Package did not support output of z-displacement arrays for models using the DISU package. The CSUB package was updated to support calculation of z-displacement arrays for DISU model grids. + %\underline{INTERNAL FLOW PACKAGES} + %\begin{itemize} % \item xxx - \end{itemize} + % \item xxx + % \item xxx + %\end{itemize} - \underline{STRESS PACKAGES} - \begin{itemize} - \item This release contains a fix for a longstanding issue related to the use of AUXMULTNAME and time series. Previous release notes included the following description of a known issue: \textit{``The AUXMULTNAME option can be used to scale input values, such as riverbed conductance, using values in an auxiliary column. When this AUXMULTNAME option is used, the multiplier value in the AUXMULTNAME column should not be represented with a time series unless the value to scale is also represented with a time series.''} With this release, the Input Data Processor (IDP) is now used to read stress package input files, and the limitation with AUXMULTNAME and time series no longer applies. + %\underline{STRESS PACKAGES} + %\begin{itemize} % \item xxx % \item xxx - \end{itemize} + % \item xxx + %\end{itemize} \underline{ADVANCED STRESS PACKAGES} \begin{itemize} - \item Added functionality to support zero values for each grid dimension when specifying the CELLID for SFR reaches that are not connected to an underlying groundwater grid cell. For example, for a DIS grid a CELLID of 0 0 0 should be specified for reaches with no connection to a groundwater cell. Warning messages will be issued if NONE is specified for the CELLID of an unconnected reach. Specifying a CELLID of NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error. - \item Added functionality to support specification of a DNODATA (3.0E+30) BEDLEAK value for LAK package connections. This DNODATA value is used to identify lake-GWF connections where conductance is solely a function of aquifer properties in the connected GWF cell. In this case, the lakebed sediments are assumed to be absent and all resistance to flow is assumed to be within the GWF cell. Warning messages are now issued if NONE is specified for LAK package connections. Specifying a BEDLEAK value equal to NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error. - \item SFR diversion would not be updated if the outflow of its upstream reach is zero. If diversion was not zero in the previous stress period, it would report mass balance error in the SFR budget. This bug was fixed by always updating the diversion. + \item Refactoring of the Water Mover package in version 6.4.3 introduced a reduction in performance for GWF models with a large number of movers. The program was corrected so that performance is similar to previous versions. + \item Deprecation warnings in the Lake and Streamflow Routing Packages, which can be issued by MODFLOW 6 at run time, incorrectly listed the version number as 6.5.0. The warning was modified to show 6.4.3 as the version number for the deprecation. % \item xxx \end{itemize} @@ -56,9 +53,17 @@ % \item xxx %\end{itemize} - \underline{EXCHANGES} + %\underline{EXCHANGES} + %\begin{itemize} + % \item xxx + % \item xxx + % \item xxx + %\end{itemize} + + \underline{PARALLEL} \begin{itemize} - \item A model budget error would occur when a constant-head (CHD) cell in one model had a direct connection to an active cell in another model. For the model budget to be calculated correctly a new term called ``FLOW-JA-FACE-CHD'' was added to the GWF model budget. This term is only included in the budget table when the GWF Model is connected to another GWF Model using a GWF-GWF Exchange. Additionally, the CHD budget calculation for a very specific (and rare) configuration was also incorrect. The incorrect budget calculation occurred when the following conditions were met: (1) a GWF model was connected to another GWF model with a GWF-GWF Exchange; (2) the model as well as the Exchange had the XT3D option enabled, and (3) the model was configured with a CHD cell that is either an Exchange cell, i.e. a cell that is part of the EXCHANGEDATA block, or a cell directly connected to such an Exchange cell. The size of the error depends on the degree of anisotropy around the particular CHD cell and shows up as a discrepancy in the volume budget table reported in the GWF list file. The program has been updated with the correct budget calculation. + \item A memory leak was identified in the way MPI messages were constructed and cached. The memory leak caused problems for large models with many time steps. The message construction and caching was fixed and tests confirm that the memory leak is no longer present. % \item xxx % \item xxx \end{itemize} + diff --git a/doc/ReleaseNotes/previous/v6.4.3.tex b/doc/ReleaseNotes/previous/v6.4.3.tex new file mode 100644 index 00000000000..3140ebd5b07 --- /dev/null +++ b/doc/ReleaseNotes/previous/v6.4.3.tex @@ -0,0 +1,65 @@ +% Use this template for starting initializing the release notes +% after a release has just been made. + + %\item \currentmodflowversion + \subsection{Version mf6.4.3--February 7, 2024} + + \underline{NEW FUNCTIONALITY} + \begin{itemize} + \item The Input Data Processor (IDP), first released in version 6.4.2, is a general utility for reading user-provided input files. Package-specific routines for reading input files continue to be replaced by the IDP approach. For packages that use IDP for input, logging information is written to the simulation list file (mfsim.lst). Additional information on the IDP and the list of supported packages is contained in the MODFLOW 6 Description of Input and Output (mf6io.pdf) under a section titled ``Processing of Program Input.'' + \item The source code was refactored to support compilation of a parallel version of MODFLOW 6 based on the Message Passing Interface (MPI) and the Portable, Extensible Toolkit for Scientific Computation (PETSc) libraries. The parallel version of MODFLOW is considered preliminary. Limited testing of the parallel version has been performed on laptops, desktops, and supercomputers, but significant changes are expected in future releases. User support for the parallel version of MODFLOW 6 may be provided in the future. + % \item xxx + \end{itemize} + + \underline{EXAMPLES} + \begin{itemize} + \item A new exampled called ex-gwf-rad-disu was added. This new example uses a DISU grid to represent radial groundwater flow to a pumping well. + \item A new exampled called ex-gwf-curv-90 was added. This new example demonstrates use of a DISV grid to represent a curvilinear spatial discretization. For this example, the curvilinear grid is applied to one quarter of a radial groundwater flow system. + \item A new exampled called ex-gwf-curvilin was added. This new example uses a curvilinear grid, represented with the DISV Package, to simulate groundwater flow through a multi-region aquifer with bends in the domain boundaries. + \end{itemize} + + \textbf{\underline{BUG FIXES AND OTHER CHANGES TO EXISTING FUNCTIONALITY}} \\ + \underline{BASIC FUNCTIONALITY} + \begin{itemize} + \item Improve error message if the size of data read from a binary array file is inconsistent with READARRAY control line and variable description keywords. + \item The area calculation for cells in the DISV package was inaccurate for some cases with very large cell vertex coordinates. The area calculation was improved by using transformed cell vertex coordinates prior to making the area calculation. + \item Auxiliary variables in RCH and EVT Array-Based input packages are now reset to zero when otherwise not specified in period input data and the auxiliary parameter is not controlled by a time-series. + % \item xxx + % \item xxx + \end{itemize} + + \underline{INTERNAL FLOW PACKAGES} + \begin{itemize} + \item The data header in the binary output file written by the viscosity (VSC) package was printing `` VISCOSI'' instead of `` VISCOSITY''. The viscosity package now prints the full `` VISCOSITY'' header in the binary output file. + \item The CSUB Package did not support output of z-displacement arrays for models using the DISU package. The CSUB package was updated to support calculation of z-displacement arrays for DISU model grids. + % \item xxx + \end{itemize} + + \underline{STRESS PACKAGES} + \begin{itemize} + \item This release contains a fix for a longstanding issue related to the use of AUXMULTNAME and time series. Previous release notes included the following description of a known issue: \textit{``The AUXMULTNAME option can be used to scale input values, such as riverbed conductance, using values in an auxiliary column. When this AUXMULTNAME option is used, the multiplier value in the AUXMULTNAME column should not be represented with a time series unless the value to scale is also represented with a time series.''} With this release, the Input Data Processor (IDP) is now used to read stress package input files, and the limitation with AUXMULTNAME and time series no longer applies. + % \item xxx + % \item xxx + \end{itemize} + + \underline{ADVANCED STRESS PACKAGES} + \begin{itemize} + \item Added functionality to support zero values for each grid dimension when specifying the CELLID for SFR reaches that are not connected to an underlying groundwater grid cell. For example, for a DIS grid a CELLID of 0 0 0 should be specified for reaches with no connection to a groundwater cell. Warning messages will be issued if NONE is specified for the CELLID of an unconnected reach. Specifying a CELLID of NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error. + \item Added functionality to support specification of a DNODATA (3.0E+30) BEDLEAK value for LAK package connections. This DNODATA value is used to identify lake-GWF connections where conductance is solely a function of aquifer properties in the connected GWF cell. In this case, the lakebed sediments are assumed to be absent and all resistance to flow is assumed to be within the GWF cell. Warning messages are now issued if NONE is specified for LAK package connections. Specifying a BEDLEAK value equal to NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error. + \item SFR diversion would not be updated if the outflow of its upstream reach is zero. If diversion was not zero in the previous stress period, it would report mass balance error in the SFR budget. This bug was fixed by always updating the diversion. + % \item xxx + \end{itemize} + + %\underline{SOLUTION} + %\begin{itemize} + % \item xxx + % \item xxx + % \item xxx + %\end{itemize} + + \underline{EXCHANGES} + \begin{itemize} + \item A model budget error would occur when a constant-head (CHD) cell in one model had a direct connection to an active cell in another model. For the model budget to be calculated correctly a new term called ``FLOW-JA-FACE-CHD'' was added to the GWF model budget. This term is only included in the budget table when the GWF Model is connected to another GWF Model using a GWF-GWF Exchange. Additionally, the CHD budget calculation for a very specific (and rare) configuration was also incorrect. The incorrect budget calculation occurred when the following conditions were met: (1) a GWF model was connected to another GWF model with a GWF-GWF Exchange; (2) the model as well as the Exchange had the XT3D option enabled, and (3) the model was configured with a CHD cell that is either an Exchange cell, i.e. a cell that is part of the EXCHANGEDATA block, or a cell directly connected to such an Exchange cell. The size of the error depends on the degree of anisotropy around the particular CHD cell and shows up as a discrepancy in the volume budget table reported in the GWF list file. The program has been updated with the correct budget calculation. + % \item xxx + % \item xxx + \end{itemize} diff --git a/doc/ReleaseNotes/vx.x.x-template.tex b/doc/ReleaseNotes/vx.x.x-template.tex index 7771001dfba..467f4f458ca 100644 --- a/doc/ReleaseNotes/vx.x.x-template.tex +++ b/doc/ReleaseNotes/vx.x.x-template.tex @@ -60,3 +60,10 @@ % \item xxx %\end{itemize} + %\underline{PARALLEL} + %\begin{itemize} + % \item xxx + % \item xxx + % \item xxx + %\end{itemize} + diff --git a/doc/mf6io/mf6ivar/dfn/exg-gwfgwf.dfn b/doc/mf6io/mf6ivar/dfn/exg-gwfgwf.dfn index 2232bfbdf90..e27195cb393 100644 --- a/doc/mf6io/mf6ivar/dfn/exg-gwfgwf.dfn +++ b/doc/mf6io/mf6ivar/dfn/exg-gwfgwf.dfn @@ -8,7 +8,8 @@ shape (naux) reader urword optional true longname keyword to specify aux variables -description an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used as the straight-line connection distance, including the vertical component, between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. +description an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' and when it is required is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used in the calculation of specific discharge within model cells connected by the exchange. For a horizontal connection, CDIST should be specified as the horizontal distance between the cell centers, and should not include the vertical component. For vertical connections, CDIST should be specified as the difference in elevation between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. + block options name boundnames diff --git a/doc/mf6io/mf6ivar/dfn/exg-gwtgwt.dfn b/doc/mf6io/mf6ivar/dfn/exg-gwtgwt.dfn index f5aeb241996..7d8b419ddd1 100644 --- a/doc/mf6io/mf6ivar/dfn/exg-gwtgwt.dfn +++ b/doc/mf6io/mf6ivar/dfn/exg-gwtgwt.dfn @@ -24,7 +24,7 @@ shape (naux) reader urword optional true longname keyword to specify aux variables -description an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used as the straight-line connection distance, including the vertical component, between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. +description an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWT-GWT Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. ANGLDEGX must be specified if dispersion is simulated in the connected GWT models. block options name boundnames diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index a1a6b95f5ac..011c23b068b 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -29,7 +29,7 @@ | SIM | TDIS | PERIODDATA | PERLEN | DOUBLE PRECISION | is the length of a stress period. | | SIM | TDIS | PERIODDATA | NSTP | INTEGER | is the number of time steps in a stress period. | | SIM | TDIS | PERIODDATA | TSMULT | DOUBLE PRECISION | is the multiplier for the length of successive time steps. The length of a time step is calculated by multiplying the length of the previous time step by TSMULT. The length of the first time step, $\Delta t_1$, is related to PERLEN, NSTP, and TSMULT by the relation $\Delta t_1= perlen \frac{tsmult - 1}{tsmult^{nstp}-1}$. | -| EXG | GWFGWF | OPTIONS | AUXILIARY | STRING (NAUX) | an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used as the straight-line connection distance, including the vertical component, between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. | +| EXG | GWFGWF | OPTIONS | AUXILIARY | STRING (NAUX) | an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' and when it is required is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used in the calculation of specific discharge within model cells connected by the exchange. For a horizontal connection, CDIST should be specified as the horizontal distance between the cell centers, and should not include the vertical component. For vertical connections, CDIST should be specified as the difference in elevation between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. | | EXG | GWFGWF | OPTIONS | BOUNDNAMES | KEYWORD | keyword to indicate that boundary names may be provided with the list of GWF Exchange cells. | | EXG | GWFGWF | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of exchange entries will be echoed to the listing file immediately after it is read. | | EXG | GWFGWF | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of exchange flow rates will be printed to the listing file for every stress period in which ``SAVE BUDGET'' is specified in Output Control. | @@ -58,7 +58,7 @@ | EXG | GWFGWF | EXCHANGEDATA | BOUNDNAME | STRING | name of the GWF Exchange cell. BOUNDNAME is an ASCII character variable that can contain as many as 40 characters. If BOUNDNAME contains spaces in it, then the entire name must be enclosed within single quotes. | | EXG | GWTGWT | OPTIONS | GWFMODELNAME1 | STRING | keyword to specify name of first corresponding GWF Model. In the simulation name file, the GWT6-GWT6 entry contains names for GWT Models (exgmnamea and exgmnameb). The GWT Model with the name exgmnamea must correspond to the GWF Model with the name gwfmodelname1. | | EXG | GWTGWT | OPTIONS | GWFMODELNAME2 | STRING | keyword to specify name of second corresponding GWF Model. In the simulation name file, the GWT6-GWT6 entry contains names for GWT Models (exgmnamea and exgmnameb). The GWT Model with the name exgmnameb must correspond to the GWF Model with the name gwfmodelname2. | -| EXG | GWTGWT | OPTIONS | AUXILIARY | STRING (NAUX) | an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used as the straight-line connection distance, including the vertical component, between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. | +| EXG | GWTGWT | OPTIONS | AUXILIARY | STRING (NAUX) | an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWT-GWT Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. ANGLDEGX must be specified if dispersion is simulated in the connected GWT models. | | EXG | GWTGWT | OPTIONS | BOUNDNAMES | KEYWORD | keyword to indicate that boundary names may be provided with the list of GWT Exchange cells. | | EXG | GWTGWT | OPTIONS | PRINT_INPUT | KEYWORD | keyword to indicate that the list of exchange entries will be echoed to the listing file immediately after it is read. | | EXG | GWTGWT | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of exchange flow rates will be printed to the listing file for every stress period in which ``SAVE BUDGET'' is specified in Output Control. | diff --git a/doc/mf6io/mf6ivar/tex/exg-gwfgwf-desc.tex b/doc/mf6io/mf6ivar/tex/exg-gwfgwf-desc.tex index 65778b0847b..825d7a89e7a 100644 --- a/doc/mf6io/mf6ivar/tex/exg-gwfgwf-desc.tex +++ b/doc/mf6io/mf6ivar/tex/exg-gwfgwf-desc.tex @@ -3,7 +3,7 @@ \item \textbf{Block: OPTIONS} \begin{description} -\item \texttt{auxiliary}---an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used as the straight-line connection distance, including the vertical component, between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. +\item \texttt{auxiliary}---an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' and when it is required is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used in the calculation of specific discharge within model cells connected by the exchange. For a horizontal connection, CDIST should be specified as the horizontal distance between the cell centers, and should not include the vertical component. For vertical connections, CDIST should be specified as the difference in elevation between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. \item \texttt{BOUNDNAMES}---keyword to indicate that boundary names may be provided with the list of GWF Exchange cells. diff --git a/doc/mf6io/mf6ivar/tex/exg-gwtgwt-desc.tex b/doc/mf6io/mf6ivar/tex/exg-gwtgwt-desc.tex index 85c90ef8aea..d58eb0416e7 100644 --- a/doc/mf6io/mf6ivar/tex/exg-gwtgwt-desc.tex +++ b/doc/mf6io/mf6ivar/tex/exg-gwtgwt-desc.tex @@ -7,7 +7,7 @@ \item \texttt{gwfmodelname2}---keyword to specify name of second corresponding GWF Model. In the simulation name file, the GWT6-GWT6 entry contains names for GWT Models (exgmnamea and exgmnameb). The GWT Model with the name exgmnameb must correspond to the GWF Model with the name gwfmodelname2. -\item \texttt{auxiliary}---an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWF-GWF Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. If an auxiliary variable with the name ``CDIST'' is found, then this information will be used as the straight-line connection distance, including the vertical component, between the two cell centers. Both ANGLDEGX and CDIST are required if specific discharge is calculated for either of the groundwater models. +\item \texttt{auxiliary}---an array of auxiliary variable names. There is no limit on the number of auxiliary variables that can be provided. Most auxiliary variables will not be used by the GWT-GWT Exchange, but they will be available for use by other parts of the program. If an auxiliary variable with the name ``ANGLDEGX'' is found, then this information will be used as the angle (provided in degrees) between the connection face normal and the x axis, where a value of zero indicates that a normal vector points directly along the positive x axis. The connection face normal is a normal vector on the cell face shared between the cell in model 1 and the cell in model 2 pointing away from the model 1 cell. Additional information on ``ANGLDEGX'' is provided in the description of the DISU Package. ANGLDEGX must be specified if dispersion is simulated in the connected GWT models. \item \texttt{BOUNDNAMES}---keyword to indicate that boundary names may be provided with the list of GWT Exchange cells. diff --git a/doc/version.py b/doc/version.py index 5456e79843b..b370e02776c 100644 --- a/doc/version.py +++ b/doc/version.py @@ -1,3 +1,3 @@ # MODFLOW 6 version file automatically created using...update_version.py -# created on...February 07, 2024 21:53:23 -__version__ = "6.4.3" +# created on...February 13, 2024 18:47:21 +__version__ = "6.4.4" diff --git a/doc/version.tex b/doc/version.tex index 88c37fbb614..fea0bc07fb1 100644 --- a/doc/version.tex +++ b/doc/version.tex @@ -1,3 +1,3 @@ -\newcommand{\modflowversion}{mf6.4.3} -\newcommand{\modflowdate}{February 07, 2024} +\newcommand{\modflowversion}{mf6.4.4} +\newcommand{\modflowdate}{February 13, 2024} \newcommand{\currentmodflowversion}{Version \modflowversion---\modflowdate} diff --git a/make/makefile b/make/makefile index 885d0663ab9..20d30e6d03f 100644 --- a/make/makefile +++ b/make/makefile @@ -5,36 +5,36 @@ include ./makedefaults # Define the source file directories SOURCEDIR1=../src -SOURCEDIR2=../src/Exchange -SOURCEDIR3=../src/Model -SOURCEDIR4=../src/Model/Geometry -SOURCEDIR5=../src/Model/TransportModel +SOURCEDIR2=../src/Model +SOURCEDIR3=../src/Model/TransportModel +SOURCEDIR4=../src/Model/GroundWaterFlow +SOURCEDIR5=../src/Model/Geometry SOURCEDIR6=../src/Model/ModelUtilities -SOURCEDIR7=../src/Model/Connection -SOURCEDIR8=../src/Model/GroundWaterTransport -SOURCEDIR9=../src/Model/GroundWaterFlow -SOURCEDIR10=../src/Distributed -SOURCEDIR11=../src/Solution -SOURCEDIR12=../src/Solution/PETSc -SOURCEDIR13=../src/Solution/LinearMethods -SOURCEDIR14=../src/Timing -SOURCEDIR15=../src/Utilities -SOURCEDIR16=../src/Utilities/TimeSeries -SOURCEDIR17=../src/Utilities/Libraries -SOURCEDIR18=../src/Utilities/Libraries/rcm -SOURCEDIR19=../src/Utilities/Libraries/sparsekit -SOURCEDIR20=../src/Utilities/Libraries/sparskit2 -SOURCEDIR21=../src/Utilities/Libraries/blas -SOURCEDIR22=../src/Utilities/Libraries/daglib -SOURCEDIR23=../src/Utilities/Idm -SOURCEDIR24=../src/Utilities/Idm/selector -SOURCEDIR25=../src/Utilities/Idm/mf6blockfile -SOURCEDIR26=../src/Utilities/Matrix -SOURCEDIR27=../src/Utilities/Vector -SOURCEDIR28=../src/Utilities/Observation -SOURCEDIR29=../src/Utilities/OutputControl -SOURCEDIR30=../src/Utilities/Memory -SOURCEDIR31=../src/Utilities/ArrayRead +SOURCEDIR7=../src/Model/GroundWaterTransport +SOURCEDIR8=../src/Model/Connection +SOURCEDIR9=../src/Distributed +SOURCEDIR10=../src/Utilities +SOURCEDIR11=../src/Utilities/Idm +SOURCEDIR12=../src/Utilities/Idm/mf6blockfile +SOURCEDIR13=../src/Utilities/Idm/selector +SOURCEDIR14=../src/Utilities/Vector +SOURCEDIR15=../src/Utilities/Matrix +SOURCEDIR16=../src/Utilities/Observation +SOURCEDIR17=../src/Utilities/ArrayRead +SOURCEDIR18=../src/Utilities/OutputControl +SOURCEDIR19=../src/Utilities/Libraries +SOURCEDIR20=../src/Utilities/Libraries/blas +SOURCEDIR21=../src/Utilities/Libraries/rcm +SOURCEDIR22=../src/Utilities/Libraries/sparsekit +SOURCEDIR23=../src/Utilities/Libraries/sparskit2 +SOURCEDIR24=../src/Utilities/Libraries/daglib +SOURCEDIR25=../src/Utilities/Memory +SOURCEDIR26=../src/Utilities/TimeSeries +SOURCEDIR27=../src/Timing +SOURCEDIR28=../src/Solution +SOURCEDIR29=../src/Solution/PETSc +SOURCEDIR30=../src/Solution/LinearMethods +SOURCEDIR31=../src/Exchange VPATH = \ ${SOURCEDIR1} \ @@ -237,6 +237,7 @@ $(OBJDIR)/SparseMatrix.o \ $(OBJDIR)/LinearSolverBase.o \ $(OBJDIR)/ims8reordering.o \ $(OBJDIR)/ModflowInput.o \ +$(OBJDIR)/IdmLogger.o \ $(OBJDIR)/Integer2dReader.o \ $(OBJDIR)/VirtualExchange.o \ $(OBJDIR)/GridSorting.o \ @@ -268,7 +269,6 @@ $(OBJDIR)/RouterBase.o \ $(OBJDIR)/ImsLinearSolver.o \ $(OBJDIR)/ims8base.o \ $(OBJDIR)/StructVector.o \ -$(OBJDIR)/IdmLogger.o \ $(OBJDIR)/DefinitionSelect.o \ $(OBJDIR)/InputLoadType.o \ $(OBJDIR)/Integer1dReader.o \ diff --git a/meson.build b/meson.build index f1bca2b04f0..026dbf29972 100644 --- a/meson.build +++ b/meson.build @@ -1,7 +1,7 @@ project( 'MODFLOW 6', 'fortran', - version: '6.4.3', + version: '6.4.4', license: 'CC0', meson_version: '>= 1.1.0', default_options : [ diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index 04ba48ba0dd..f0a03e6f2da 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -54,6 +54,11 @@ + + + + + @@ -64,6 +69,11 @@ + + + + + diff --git a/pymake/excludefiles.txt b/pymake/excludefiles.txt index 33694f4ea5d..c2677d02b3d 100644 --- a/pymake/excludefiles.txt +++ b/pymake/excludefiles.txt @@ -5,6 +5,8 @@ ../src/Utilities/Matrix/PetscMatrix.F90 ../src/Utilities/Vector/PetscVector.F90 ../src/Distributed/MpiMessageBuilder.f90 +../src/Distributed/MpiMessageCache.f90 ../src/Distributed/MpiRouter.f90 ../src/Distributed/MpiRunControl.F90 +../src/Distributed/MpiUnitCache.f90 ../src/Distributed/MpiWorld.f90 diff --git a/src/Distributed/MpiMessageBuilder.f90 b/src/Distributed/MpiMessageBuilder.f90 index b92165923f8..8c9be363329 100644 --- a/src/Distributed/MpiMessageBuilder.f90 +++ b/src/Distributed/MpiMessageBuilder.f90 @@ -42,7 +42,6 @@ module MpiMessageBuilderModule procedure, private :: create_vdc_snd_map procedure, private :: create_vdc_snd_body procedure, private :: create_vdc_rcv_body - procedure, private :: create_element_map end type contains @@ -114,7 +113,7 @@ subroutine create_header_snd(this, rank, stage, hdrs_snd_type) class(MpiMessageBuilderType) :: this integer(I4B) :: rank integer(I4B) :: stage - integer :: hdrs_snd_type + integer, intent(out) :: hdrs_snd_type ! local integer(I4B) :: i, offset, nr_types class(VirtualDataContainerType), pointer :: vdc @@ -186,7 +185,7 @@ end subroutine create_header_snd subroutine create_header_rcv(this, hdr_rcv_type) class(MpiMessageBuilderType) :: this - integer :: hdr_rcv_type + integer, intent(out) :: hdr_rcv_type ! local integer :: ierr @@ -203,7 +202,7 @@ subroutine create_map_snd(this, rank, stage, map_snd_type) class(MpiMessageBuilderType) :: this integer(I4B) :: rank integer(I4B) :: stage - integer :: map_snd_type + integer, intent(out) :: map_snd_type ! local integer(I4B) :: i, offset, nr_types class(VirtualDataContainerType), pointer :: vdc @@ -280,7 +279,7 @@ subroutine create_map_rcv(this, rcv_map, nr_headers, map_rcv_type) class(MpiMessageBuilderType) :: this type(VdcReceiverMapsType), dimension(:) :: rcv_map integer(I4B) :: nr_headers - integer :: map_rcv_type + integer, intent(out) :: map_rcv_type ! local integer(I4B) :: i, j, nr_elems, type_cnt integer :: ierr, max_nr_maps @@ -323,7 +322,7 @@ subroutine create_body_rcv(this, rank, stage, body_rcv_type) class(MpiMessageBuilderType) :: this integer(I4B) :: rank integer(I4B) :: stage - integer :: body_rcv_type + integer, intent(out) :: body_rcv_type ! local integer(I4B) :: i, nr_types, offset class(VirtualDataContainerType), pointer :: vdc @@ -400,7 +399,7 @@ subroutine create_body_snd(this, rank, stage, headers, maps, body_snd_type) integer(I4B) :: stage type(VdcHeaderType), dimension(:) :: headers type(VdcReceiverMapsType), dimension(:) :: maps - integer :: body_snd_type + integer, intent(out) :: body_snd_type ! local integer(I4B) :: i, nr_headers class(VirtualDataContainerType), pointer :: vdc @@ -627,33 +626,6 @@ function create_vdc_snd_body(this, vdc, vdc_maps, rank, stage) result(new_type) end function create_vdc_snd_body - !> @brief Temp. function to generate a dummy (complete) map - !< - function create_element_map(this, rank, vdc, vd) result(el_map) - use MemoryManagerModule, only: get_mem_shape, get_mem_rank - use ConstantsModule, only: MAXMEMRANK - class(MpiMessageBuilderType) :: this - integer(I4B) :: rank - class(VirtualDataContainerType), pointer :: vdc - class(VirtualDataType), pointer :: vd - integer(I4B), dimension(:), pointer, contiguous :: el_map - ! local - integer(I4B), dimension(MAXMEMRANK) :: mem_shp - integer(I4B) :: i, nrow, mem_rank - - el_map => null() - call get_mem_rank(vd%virtual_mt%name, vd%virtual_mt%path, mem_rank) - call get_mem_shape(vd%virtual_mt%name, vd%virtual_mt%path, mem_shp) - if (mem_rank > 0) then - nrow = mem_shp(mem_rank) - allocate (el_map(nrow)) - do i = 1, nrow - el_map(i) = i - 1 - end do - end if - - end function create_element_map - function get_vdc_from_hdr(this, header) result(vdc) class(MpiMessageBuilderType) :: this type(VdcHeaderType) :: header diff --git a/src/Distributed/MpiMessageCache.f90 b/src/Distributed/MpiMessageCache.f90 new file mode 100644 index 00000000000..222a4f2ac27 --- /dev/null +++ b/src/Distributed/MpiMessageCache.f90 @@ -0,0 +1,130 @@ +module MpiMessageCacheModule + use KindModule, only: I4B + use SimStagesModule, only: NR_SIM_STAGES + use ListModule + use STLVecIntModule + use MpiUnitCacheModule + implicit none + private + + ! the message types for caching during a simulation stage: + integer(I4B), public, parameter :: MPI_BDY_RCV = 1 !< receiving data (body) from ranks + integer(I4B), public, parameter :: MPI_BDY_SND = 2 !< sending data (body) to ranks + integer(I4B), public, parameter :: NR_MSG_TYPES = 2 !< the total number of message types to be cached + + ! expose this from the unit cache module + public :: NO_CACHED_VALUE + + !> @brief Facility to cache the constructed MPI datatypes. + !! This will avoid having to construct them over and over + !! again for the communication inside the timestep loop. + !! This class deals with separate caches for different + !! units (solutions or global) and for different types of + !< messages within the communication stage. + type, public :: MpiMessageCacheType + type(STLVecInt) :: cached_ids !< a vector with ids for the cached units (solution ids) + type(ListType) :: unit_caches !< a list with caches per unit + contains + procedure :: init => mmc_init + procedure :: get => mmc_get + procedure :: put => mmc_put + procedure :: destroy => mmc_destroy + end type MpiMessageCacheType + +contains + + !< @brief Initialize the MPI type cache system. + !< + subroutine mmc_init(this) + class(MpiMessageCacheType) :: this !< the message cache + + call this%cached_ids%init() + + end subroutine mmc_init + + !< @brief Get the cached mpi datatype for the given + !! unit, rank, stage, and message element. Returns + !< NO_CACHED_VALUE when not in cache. + function mmc_get(this, unit, rank, stage, msg_id) result(mpi_type) + class(MpiMessageCacheType) :: this !< the message cache + integer(I4B) :: unit !< the unit (solution or global) + integer(I4B) :: rank !< the rank of the MPI process to communicate with + integer(I4B) :: stage !< the simulation stage at which the message is sent + integer(I4B) :: msg_id !< the message type as an integer between 1 and NR_MSG_TYPES (see above for predefined values) + integer :: mpi_type !< the resulting mpi datatype + ! local + integer(I4B) :: unit_idx + class(*), pointer :: obj_ptr + + mpi_type = NO_CACHED_VALUE + + unit_idx = this%cached_ids%get_index(unit) + if (unit_idx == -1) return ! not cached + + obj_ptr => this%unit_caches%GetItem(unit_idx) + select type (obj_ptr) + class is (MpiUnitCacheType) + mpi_type = obj_ptr%get_cached(rank, stage, msg_id) + end select + + end function mmc_get + + !> @brief Put the mpi datatype for this particular unit, + !! rank, and stage in cache. The datatype should be + !< committed to the type database externally. + subroutine mmc_put(this, unit, rank, stage, msg_id, mpi_type) + class(MpiMessageCacheType) :: this !< the message cache + integer(I4B) :: unit !< the unit (solution or global) + integer(I4B) :: rank !< the rank of the MPI process to communicate with + integer(I4B) :: stage !< the simulation stage at which the message is sent + integer(I4B) :: msg_id !< the message type as an integer between 1 and NR_MSG_TYPES (see above for predefined values) + integer :: mpi_type !< the mpi datatype to cache + ! local + integer(I4B) :: unit_idx + type(MpiUnitCacheType), pointer :: new_cache + class(*), pointer :: obj_ptr + + unit_idx = this%cached_ids%get_index(unit) + if (unit_idx == -1) then + ! add to vector with cached unit ids + call this%cached_ids%push_back(unit) + ! create and add unit cache + allocate (new_cache) + call new_cache%init(NR_SIM_STAGES, NR_MSG_TYPES) + obj_ptr => new_cache + call this%unit_caches%Add(obj_ptr) + unit_idx = this%cached_ids%size + end if + + ! get the cache for this unit + obj_ptr => this%unit_caches%GetItem(unit_idx) + select type (obj_ptr) + class is (MpiUnitCacheType) + call obj_ptr%cache(rank, stage, msg_id, mpi_type) + end select + + end subroutine mmc_put + + !< @brief Destroy the MPI type cache system. + !< + subroutine mmc_destroy(this) + class(MpiMessageCacheType) :: this !< the message cache + ! local + integer(I4B) :: i + class(*), pointer :: obj_ptr + + ! clear caches + do i = 1, this%cached_ids%size + obj_ptr => this%unit_caches%GetItem(i) + select type (obj_ptr) + class is (MpiUnitCacheType) + call obj_ptr%destroy() + end select + end do + call this%unit_caches%Clear(destroy=.true.) + + call this%cached_ids%destroy() + + end subroutine mmc_destroy + +end module diff --git a/src/Distributed/MpiRouter.f90 b/src/Distributed/MpiRouter.f90 index e3cc403b832..b4fcd6c1eff 100644 --- a/src/Distributed/MpiRouter.f90 +++ b/src/Distributed/MpiRouter.f90 @@ -11,6 +11,7 @@ module MpiRouterModule use VirtualExchangeModule, only: VirtualExchangeType use VirtualSolutionModule use MpiMessageBuilderModule + use MpiMessageCacheModule use MpiWorldModule use mpi implicit none @@ -27,6 +28,7 @@ module MpiRouterModule type(VdcPtrType), dimension(:), pointer :: rte_models => null() !< the currently active models to be routed type(VdcPtrType), dimension(:), pointer :: rte_exchanges => null() !< the currently active exchanges to be routed type(MpiMessageBuilderType) :: message_builder + type(MpiMessageCacheType) :: msg_cache type(MpiWorldType), pointer :: mpi_world => null() integer(I4B) :: imon !< the output file unit for the mpi monitor logical(LGP) :: enable_monitor !< when true, log diagnostics @@ -38,11 +40,14 @@ module MpiRouterModule ! private procedure, private :: activate procedure, private :: deactivate - procedure, private :: mr_update_senders - procedure, private :: mr_update_senders_sln - procedure, private :: mr_update_receivers - procedure, private :: mr_update_receivers_sln - procedure, private :: mr_route_active + procedure, private :: update_senders + procedure, private :: update_senders_sln + procedure, private :: update_receivers + procedure, private :: update_receivers_sln + procedure, private :: route_active + procedure, private :: is_cached + procedure, private :: compose_messages + procedure, private :: load_messages end type MpiRouterType contains @@ -73,8 +78,9 @@ subroutine mr_initialize(this) ! to log or not to log this%enable_monitor = .false. - ! initialize the MPI message builder + ! initialize the MPI message builder and cache call this%message_builder%init() + call this%msg_cache%init() ! get mpi world for our process this%mpi_world => get_mpi_world() @@ -102,6 +108,7 @@ subroutine mr_initialize(this) call MPI_Allreduce(MPI_IN_PLACE, this%model_proc_ids, nr_models, & MPI_INTEGER, MPI_SUM, this%mpi_world%comm, ierr) + call CHECK_MPI(ierr) ! set the process id to the models and exchanges do i = 1, nr_models @@ -181,7 +188,7 @@ subroutine mr_route_all(this, stage) ! route all call this%activate(this%all_models, this%all_exchanges) - call this%mr_route_active(stage) + call this%route_active(0, stage) call this%deactivate() if (this%enable_monitor) then @@ -206,7 +213,7 @@ subroutine mr_route_sln(this, virtual_sol, stage) ! route for this solution call this%activate(virtual_sol%models, virtual_sol%exchanges) - call this%mr_route_active(stage) + call this%route_active(virtual_sol%solution_id, stage) call this%deactivate() if (this%enable_monitor) then @@ -215,39 +222,50 @@ subroutine mr_route_sln(this, virtual_sol, stage) end subroutine mr_route_sln - !> @brief Routes the models and exchanges. This is the - !< workhorse routine - subroutine mr_route_active(this, stage) - class(MpiRouterType) :: this - integer(I4B) :: stage + !> @brief Routes the models and exchanges over MPI, + !! either constructing the message bodies in a sequence + !! of communication steps, or by loading from cache + !< for efficiency. + subroutine route_active(this, unit, stage) + use SimModule, only: store_error + class(MpiRouterType) :: this !< this mpi router + integer(I4B) :: unit !< the solution id, or equal to 0 when global + integer(I4B) :: stage !< the stage to route ! local - integer(I4B) :: i, j, k + integer(I4B) :: i integer(I4B) :: rnk integer :: ierr, msg_size + logical(LGP) :: from_cache ! mpi handles integer, dimension(:), allocatable :: rcv_req integer, dimension(:), allocatable :: snd_req integer, dimension(:, :), allocatable :: rcv_stat - integer, dimension(:, :), allocatable :: snd_stat - ! message header - integer(I4B) :: max_headers - type(VdcHeaderType), dimension(:, :), allocatable :: headers - integer, dimension(:), allocatable :: hdr_rcv_t - integer, dimension(:), allocatable :: hdr_snd_t - integer, dimension(:), allocatable :: hdr_rcv_cnt - - ! maps - type(VdcReceiverMapsType), dimension(:, :), allocatable :: rcv_maps - integer, dimension(:), allocatable :: map_rcv_t - integer, dimension(:), allocatable :: map_snd_t ! message body integer, dimension(:), allocatable :: body_rcv_t integer, dimension(:), allocatable :: body_snd_t ! update adress list - call this%mr_update_senders() - call this%mr_update_receivers() + call this%update_senders() + call this%update_receivers() + + if (this%senders%size /= this%receivers%size) then + call store_error("Internal error: receivers should equal senders", & + terminate=.true.) + end if + + ! allocate body data + allocate (body_rcv_t(this%senders%size)) + allocate (body_snd_t(this%receivers%size)) + + ! allocate handles + allocate (rcv_req(this%receivers%size)) + allocate (snd_req(this%senders%size)) + allocate (rcv_stat(MPI_STATUS_SIZE, this%receivers%size)) + + ! always initialize request handles + rcv_req = MPI_REQUEST_NULL + snd_req = MPI_REQUEST_NULL if (this%enable_monitor) then write (this%imon, '(2x,a,*(i3))') "process ids sending data: ", & @@ -256,11 +274,132 @@ subroutine mr_route_active(this, stage) this%receivers%get_values() end if + ! from cache or build + from_cache = this%is_cached(unit, stage) + if (.not. from_cache) then + call this%compose_messages(unit, stage, body_snd_t, body_rcv_t) + else + call this%load_messages(unit, stage, body_snd_t, body_rcv_t) + end if + + if (this%enable_monitor) then + write (this%imon, '(2x,a)') "== communicating bodies ==" + end if + + ! recv bodies + do i = 1, this%senders%size + rnk = this%senders%at(i) + if (this%enable_monitor) then + write (this%imon, '(4x,a,i0)') "receiving from process: ", rnk + end if + + call MPI_Type_size(body_rcv_t(i), msg_size, ierr) + if (msg_size > 0) then + call MPI_Irecv(MPI_BOTTOM, 1, body_rcv_t(i), rnk, stage, & + this%mpi_world%comm, rcv_req(i), ierr) + call CHECK_MPI(ierr) + end if + + if (this%enable_monitor) then + write (this%imon, '(6x,a,i0)') "message body size: ", msg_size + end if + end do + + ! send bodies + do i = 1, this%receivers%size + rnk = this%receivers%at(i) + if (this%enable_monitor) then + write (this%imon, '(4x,a,i0)') "sending to process: ", rnk + end if + + call MPI_Type_size(body_snd_t(i), msg_size, ierr) + if (msg_size > 0) then + call MPI_Isend(MPI_Bottom, 1, body_snd_t(i), rnk, stage, & + this%mpi_world%comm, snd_req(i), ierr) + call CHECK_MPI(ierr) + end if + + if (this%enable_monitor) then + write (this%imon, '(6x,a,i0)') "message body size: ", msg_size + end if + call flush (this%imon) + end do + + ! wait for exchange of all messages + call MPI_WaitAll(this%senders%size, rcv_req, rcv_stat, ierr) + call CHECK_MPI(ierr) + + deallocate (rcv_req, snd_req, rcv_stat) + deallocate (body_rcv_t, body_snd_t) + + end subroutine route_active + + !> @brief Constructs the message bodies' MPI datatypes. + !! + !! Constructs the message bodies' MPI datatypes for a + !! unit (a solution) and a given stage. This is done in + !! a sequence of 6 steps (distributed over 3 phases): + !! + !! == synchronizing headers (VdcHeaderType) == + !! + !! step 1: + !! Receive headers from remote addresses requesting data + !! from virtual data containers (models, exchanges, ...) local to this process + !! step 2: + !! Send headers to remote addresses to indicate for which + !! virtual data containers (models, exchanges, ...) data is requested + !! + !! == synchronizing maps (VdcReceiverMapsType) == + !! + !! step 3: + !! Based on the received headers, receive element maps (which elements are + !! to be sent from a contiguous array) for outgoing data + !! step 4: + !! Send element maps to remote addresses to specify incoming data + !! + !! == construct message body data types (VirtualDataContainerType) == + !! + !! step 5: + !! Construct the message receive body based on the virtual + !! data items in the virtual data containers, and cache + !! + !! step 6: + !! Construct the message send body, based on received header data and + !! and maps, from the virtual data containers, and cache + !< + subroutine compose_messages(this, unit, stage, body_snd_t, body_rcv_t) + class(MpiRouterType) :: this + integer(I4B) :: unit + integer(I4B) :: stage + integer, dimension(:) :: body_snd_t + integer, dimension(:) :: body_rcv_t + ! local + integer(I4B) :: i, j, k + integer(I4B) :: rnk + integer :: ierr + ! mpi handles + integer, dimension(:), allocatable :: rcv_req + integer, dimension(:), allocatable :: snd_req + integer, dimension(:, :), allocatable :: rcv_stat + ! message header + integer(I4B) :: max_headers + type(VdcHeaderType), dimension(:, :), allocatable :: headers + integer, dimension(:), allocatable :: hdr_rcv_t + integer, dimension(:), allocatable :: hdr_snd_t + integer, dimension(:), allocatable :: hdr_rcv_cnt + ! maps + type(VdcReceiverMapsType), dimension(:, :), allocatable :: rcv_maps + integer, dimension(:), allocatable :: map_rcv_t + integer, dimension(:), allocatable :: map_snd_t + ! allocate handles allocate (rcv_req(this%receivers%size)) allocate (snd_req(this%senders%size)) allocate (rcv_stat(MPI_STATUS_SIZE, this%receivers%size)) - allocate (snd_stat(MPI_STATUS_SIZE, this%senders%size)) + + ! init handles + rcv_req = MPI_REQUEST_NULL + snd_req = MPI_REQUEST_NULL ! allocate header data max_headers = size(this%rte_models) + size(this%rte_exchanges) @@ -274,10 +413,6 @@ subroutine mr_route_active(this, stage) allocate (map_rcv_t(this%receivers%size)) allocate (rcv_maps(max_headers, this%receivers%size)) ! for every header, we potentially need the maps - ! allocate body data - allocate (body_rcv_t(this%senders%size)) - allocate (body_snd_t(this%receivers%size)) - if (this%enable_monitor) then write (this%imon, '(2x,a)') "== communicating headers ==" end if @@ -291,7 +426,7 @@ subroutine mr_route_active(this, stage) call this%message_builder%create_header_rcv(hdr_rcv_t(i)) call MPI_Irecv(headers(:, i), max_headers, hdr_rcv_t(i), rnk, stage, & this%mpi_world%comm, rcv_req(i), ierr) - ! don't free mpi datatype, we need the count below + call CHECK_MPI(ierr) end do ! send header for incoming data @@ -303,11 +438,16 @@ subroutine mr_route_active(this, stage) call this%message_builder%create_header_snd(rnk, stage, hdr_snd_t(i)) call MPI_Isend(MPI_BOTTOM, 1, hdr_snd_t(i), rnk, stage, & this%mpi_world%comm, snd_req(i), ierr) - call MPI_Type_free(hdr_snd_t(i), ierr) + call CHECK_MPI(ierr) end do ! wait for exchange of all headers call MPI_WaitAll(this%receivers%size, rcv_req, rcv_stat, ierr) + call CHECK_MPI(ierr) + + ! reinit handles + rcv_req = MPI_REQUEST_NULL + snd_req = MPI_REQUEST_NULL ! after WaitAll we can count incoming headers from statuses do i = 1, this%receivers%size @@ -323,9 +463,15 @@ subroutine mr_route_active(this, stage) write (this%imon, '(6x,a,99i6)') "map sizes: ", headers(j, i)%map_sizes end do end if + end do + ! clean up types + do i = 1, this%receivers%size call MPI_Type_free(hdr_rcv_t(i), ierr) end do + do i = 1, this%senders%size + call MPI_Type_free(hdr_snd_t(i), ierr) + end do if (this%enable_monitor) then write (this%imon, '(2x,a)') "== communicating maps ==" @@ -347,9 +493,9 @@ subroutine mr_route_active(this, stage) call this%message_builder%create_map_rcv(rcv_maps(:, i), hdr_rcv_cnt(i), & map_rcv_t(i)) - call MPI_Irecv(MPI_BOTTOM, 1, map_rcv_t(i), rnk, stage, & this%mpi_world%comm, rcv_req(i), ierr) + call CHECK_MPI(ierr) end do ! send maps @@ -358,13 +504,16 @@ subroutine mr_route_active(this, stage) if (this%enable_monitor) then write (this%imon, '(4x,a,i0)') "send map to process: ", rnk end if + call this%message_builder%create_map_snd(rnk, stage, map_snd_t(i)) call MPI_Isend(MPI_BOTTOM, 1, map_snd_t(i), rnk, stage, & this%mpi_world%comm, snd_req(i), ierr) + call CHECK_MPI(ierr) end do ! wait on receiving maps call MPI_WaitAll(this%receivers%size, rcv_req, rcv_stat, ierr) + call CHECK_MPI(ierr) ! print maps if (this%enable_monitor) then @@ -395,77 +544,75 @@ subroutine mr_route_active(this, stage) end do if (this%enable_monitor) then - write (this%imon, '(2x,a)') "== communicating bodies ==" + write (this%imon, '(2x,a)') "== composing message bodies ==" end if - ! recv bodies + ! construct recv bodies and cache do i = 1, this%senders%size rnk = this%senders%at(i) if (this%enable_monitor) then - write (this%imon, '(4x,a,i0)') "receiving from process: ", rnk + write (this%imon, '(4x,a,i0)') "build recv body for process: ", rnk end if call this%message_builder%create_body_rcv(rnk, stage, body_rcv_t(i)) - call MPI_Type_size(body_rcv_t(i), msg_size, ierr) - if (msg_size > 0) then - call MPI_Irecv(MPI_BOTTOM, 1, body_rcv_t(i), rnk, stage, & - this%mpi_world%comm, rcv_req(i), ierr) - end if - - if (this%enable_monitor) then - write (this%imon, '(6x,a,i0)') "message body size: ", msg_size - end if + call this%msg_cache%put(unit, rnk, stage, MPI_BDY_RCV, body_rcv_t(i)) end do - ! send bodies + ! construct send bodies and cache do i = 1, this%receivers%size rnk = this%receivers%at(i) if (this%enable_monitor) then - write (this%imon, '(4x,a,i0)') "sending to process: ", rnk + write (this%imon, '(4x,a,i0)') "build send body for process: ", rnk end if + call this%message_builder%create_body_snd( & rnk, stage, headers(1:hdr_rcv_cnt(i), i), & rcv_maps(:, i), body_snd_t(i)) - call MPI_Type_size(body_snd_t(i), msg_size, ierr) - if (msg_size > 0) then - call MPI_Isend(MPI_Bottom, 1, body_snd_t(i), rnk, stage, & - this%mpi_world%comm, snd_req(i), ierr) - end if - - if (this%enable_monitor) then - write (this%imon, '(6x,a,i0)') "message body size: ", msg_size - end if - call flush (this%imon) - end do - - ! wait for exchange of all messages - call MPI_WaitAll(this%senders%size, rcv_req, snd_stat, ierr) - - ! clean up types - do i = 1, this%senders%size - call MPI_Type_free(body_rcv_t(i), ierr) - end do - do i = 1, this%receivers%size - call MPI_Type_free(body_snd_t(i), ierr) + call this%msg_cache%put(unit, rnk, stage, MPI_BDY_SND, body_snd_t(i)) end do - ! done sending, clean up element maps + ! clean up element maps do i = 1, this%receivers%size do j = 1, hdr_rcv_cnt(i) call rcv_maps(j, i)%destroy() end do end do - deallocate (rcv_req, snd_req, rcv_stat, snd_stat) + deallocate (rcv_req, snd_req, rcv_stat) deallocate (hdr_rcv_t, hdr_snd_t, hdr_rcv_cnt) deallocate (headers) deallocate (map_rcv_t, map_snd_t) deallocate (rcv_maps) - deallocate (body_rcv_t, body_snd_t) - end subroutine mr_route_active + end subroutine compose_messages - subroutine mr_update_senders(this) + !> @brief Load the message body MPI datatypes from cache + !< + subroutine load_messages(this, unit, stage, body_snd_t, body_rcv_t) + class(MpiRouterType) :: this + integer(I4B) :: unit + integer(I4B) :: stage + integer, dimension(:), allocatable :: body_snd_t + integer, dimension(:), allocatable :: body_rcv_t + ! local + integer(I4B) :: i, rnk + + if (this%enable_monitor) then + write (this%imon, '(2x,a)') "... running from cache ..." + end if + + do i = 1, this%receivers%size + rnk = this%receivers%at(i) + body_snd_t(i) = this%msg_cache%get(unit, rnk, stage, MPI_BDY_SND) + end do + do i = 1, this%senders%size + rnk = this%senders%at(i) + body_rcv_t(i) = this%msg_cache%get(unit, rnk, stage, MPI_BDY_RCV) + end do + + end subroutine load_messages + + subroutine update_senders(this) class(MpiRouterType) :: this ! local integer(I4B) :: i @@ -486,9 +633,9 @@ subroutine mr_update_senders(this) end if end do - end subroutine mr_update_senders + end subroutine update_senders - subroutine mr_update_senders_sln(this, virtual_sol) + subroutine update_senders_sln(this, virtual_sol) class(MpiRouterType) :: this type(VirtualSolutionType) :: virtual_sol ! local @@ -510,9 +657,9 @@ subroutine mr_update_senders_sln(this, virtual_sol) end if end do - end subroutine mr_update_senders_sln + end subroutine update_senders_sln - subroutine mr_update_receivers(this) + subroutine update_receivers(this) class(MpiRouterType) :: this ! local integer(I4B) :: i @@ -524,9 +671,9 @@ subroutine mr_update_receivers(this) call this%receivers%push_back(this%senders%at(i)) end do - end subroutine mr_update_receivers + end subroutine update_receivers - subroutine mr_update_receivers_sln(this, virtual_sol) + subroutine update_receivers_sln(this, virtual_sol) class(MpiRouterType) :: this type(VirtualSolutionType) :: virtual_sol ! local @@ -539,11 +686,51 @@ subroutine mr_update_receivers_sln(this, virtual_sol) call this%receivers%push_back(this%senders%at(i)) end do - end subroutine mr_update_receivers_sln + end subroutine update_receivers_sln + + !> @brief Check if this stage is cached + !< + function is_cached(this, unit, stage) result(in_cache) + use SimModule, only: ustop + class(MpiRouterType) :: this + integer(I4B) :: unit + integer(I4B) :: stage + logical(LGP) :: in_cache + ! local + integer(I4B) :: i, rnk + integer(I4B) :: no_cache_cnt + integer :: cached_type + + in_cache = .false. + no_cache_cnt = 0 + + do i = 1, this%receivers%size + rnk = this%receivers%at(i) + cached_type = this%msg_cache%get(unit, rnk, stage, MPI_BDY_SND) + if (cached_type == NO_CACHED_VALUE) no_cache_cnt = no_cache_cnt + 1 + end do + do i = 1, this%senders%size + rnk = this%senders%at(i) + cached_type = this%msg_cache%get(unit, rnk, stage, MPI_BDY_RCV) + if (cached_type == NO_CACHED_VALUE) no_cache_cnt = no_cache_cnt + 1 + end do + + ! it should be all or nothing + if (no_cache_cnt == this%receivers%size + this%senders%size) then + in_cache = .false. + else if (no_cache_cnt == 0) then + in_cache = .true. + else + call ustop("Internal error: MPI message cache corrupt...") + end if + + end function is_cached subroutine mr_destroy(this) class(MpiRouterType) :: this + call this%msg_cache%destroy() + call this%senders%destroy() call this%receivers%destroy() diff --git a/src/Distributed/MpiRunControl.F90 b/src/Distributed/MpiRunControl.F90 index e37f50c174a..796b59bf9b6 100644 --- a/src/Distributed/MpiRunControl.F90 +++ b/src/Distributed/MpiRunControl.F90 @@ -79,6 +79,7 @@ subroutine mpi_ctrl_start(this) #else if (.not. mpi_world%has_comm()) then call MPI_Init(ierr) + call CHECK_MPI(ierr) call mpi_world%set_comm(MPI_COMM_WORLD) end if #endif @@ -130,6 +131,7 @@ subroutine mpi_ctrl_finish(this) CHKERRQ(ierr) #else call MPI_Finalize(ierr) + call CHECK_MPI(ierr) #endif ! finish everything else by calling parent diff --git a/src/Distributed/MpiUnitCache.f90 b/src/Distributed/MpiUnitCache.f90 new file mode 100644 index 00000000000..51bf4a546f1 --- /dev/null +++ b/src/Distributed/MpiUnitCache.f90 @@ -0,0 +1,166 @@ +module MpiUnitCacheModule + use KindModule, only: I4B, LGP + use ListModule + use STLVecIntModule + use SimStagesModule, only: NR_SIM_STAGES + use mpi + implicit none + private + + integer(I4B), public, parameter :: NO_CACHED_VALUE = -1 + + type, public :: MpiUnitCacheType + ! private + type(STLVecInt), private :: cached_ranks + type(STLVecInt), private :: cached_messages + integer(I4B), private :: nr_stages + integer(I4B), private :: nr_msg_types + contains + procedure :: init => cc_init + procedure :: get_cached => cc_get_cached + procedure :: cache => mc_cache + procedure :: destroy => cc_destroy + ! private + procedure, private :: is_rank_cached + procedure, private :: add_rank_cache + procedure, private :: get_rank_index + procedure, private :: get_msg_index + end type MpiUnitCacheType + +contains + + !> @brief Initialize the unit cache. + !< + subroutine cc_init(this, nr_stages, nr_msg_types) + class(MpiUnitCacheType) :: this + integer(I4B) :: nr_stages !< number of (simulation) stages + integer(I4B) :: nr_msg_types !< number of message types to be cached during a stage + + this%nr_stages = nr_stages + this%nr_msg_types = nr_msg_types + call this%cached_ranks%init() + call this%cached_messages%init() + + end subroutine cc_init + + !> @brief Get the cached mpi type for this rank and + !< stage. Equal to NO_CACHED_VALUE when not present. + function cc_get_cached(this, rank, stage, msg_id) result(mpi_type) + class(MpiUnitCacheType) :: this + integer(I4B) :: rank + integer(I4B) :: stage + integer(I4B) :: msg_id + integer :: mpi_type + ! local + integer(I4B) :: msg_idx + + mpi_type = NO_CACHED_VALUE + msg_idx = this%get_msg_index(rank, stage, msg_id) + if (msg_idx > 0) then + mpi_type = this%cached_messages%at(msg_idx) + end if + + end function cc_get_cached + + !> @brief Cache the mpi datatype for this particular + !! rank and stage. The datatype should be committed + !< to the type database externally. + subroutine mc_cache(this, rank, stage, msg_id, mpi_type) + class(MpiUnitCacheType) :: this + integer(I4B) :: rank + integer(I4B) :: stage + integer(I4B) :: msg_id + integer :: mpi_type + ! local + integer(I4B) :: msg_idx + + ! add if rank not present in cache yet + if (.not. this%is_rank_cached(rank)) then + call this%add_rank_cache(rank) + end if + + ! rank has been added to cache, now set + ! mpi datatype for this stage's message: + msg_idx = this%get_msg_index(rank, stage, msg_id) + call this%cached_messages%set(msg_idx, mpi_type) + + end subroutine mc_cache + + function is_rank_cached(this, rank) result(in_cache) + class(MpiUnitCacheType) :: this + integer(I4B) :: rank + logical(LGP) :: in_cache + + in_cache = this%cached_ranks%contains(rank) + + end function is_rank_cached + + subroutine add_rank_cache(this, rank) + class(MpiUnitCacheType) :: this + integer(I4B) :: rank + ! local + integer(I4B) :: i, j + + call this%cached_ranks%push_back(rank) + do i = 1, this%nr_stages + do j = 1, this%nr_msg_types + call this%cached_messages%push_back(NO_CACHED_VALUE) + end do + end do + + end subroutine add_rank_cache + + !> @Brief returns -1 when not present + !< + function get_rank_index(this, rank) result(rank_index) + class(MpiUnitCacheType) :: this + integer(I4B) :: rank + integer(I4B) :: rank_index + + rank_index = this%cached_ranks%get_index(rank) + + end function get_rank_index + + !> @Brief returns -1 when not present + !< + function get_msg_index(this, rank, stage, msg_id) result(msg_index) + class(MpiUnitCacheType) :: this + integer(I4B) :: rank + integer(I4B) :: stage + integer(I4B) :: msg_id + integer(I4B) :: msg_index + ! local + integer(I4B) :: rank_idx + integer(I4B) :: rank_offset, stage_offset + + msg_index = -1 + rank_idx = this%get_rank_index(rank) + if (rank_idx < 1) return + + rank_offset = (rank_idx - 1) * (this%nr_stages * this%nr_msg_types) + stage_offset = (stage - 1) * this%nr_msg_types + msg_index = rank_offset + stage_offset + msg_id + + end function get_msg_index + + !> @brief Clean up the unit cache. + !< + subroutine cc_destroy(this) + class(MpiUnitCacheType) :: this + ! local + integer(I4B) :: i + integer :: mpi_type, ierr + + do i = 1, this%cached_messages%size + mpi_type = this%cached_messages%at(i) + if (mpi_type /= NO_CACHED_VALUE) then + call MPI_Type_free(mpi_type, ierr) + end if + end do + + call this%cached_ranks%destroy() + call this%cached_messages%destroy() + + end subroutine cc_destroy + +end module diff --git a/src/Distributed/MpiWorld.f90 b/src/Distributed/MpiWorld.f90 index 8962c99d757..7ee3498467e 100644 --- a/src/Distributed/MpiWorld.f90 +++ b/src/Distributed/MpiWorld.f90 @@ -6,6 +6,7 @@ module MpiWorldModule private public :: get_mpi_world + public :: CHECK_MPI type, public :: MpiWorldType integer(I4B) :: mpi_rank !< the id for this process @@ -108,4 +109,21 @@ subroutine mpiw_destroy(this) end subroutine mpiw_destroy + !> @brief Check the MPI error code, report, and + !< terminate when not MPI_SUCCESS + subroutine CHECK_MPI(mpi_error_code) + use SimModule, only: store_error + integer :: mpi_error_code + ! local + character(len=1024) :: mpi_err_msg + integer :: err_len + integer :: ierr + + if (mpi_error_code /= MPI_SUCCESS) then + call MPI_Error_string(mpi_error_code, mpi_err_msg, err_len, ierr) + call store_error("Internal error: "//trim(mpi_err_msg), terminate=.true.) + end if + + end subroutine CHECK_MPI + end module MpiWorldModule diff --git a/src/Distributed/VirtualDataContainer.f90 b/src/Distributed/VirtualDataContainer.f90 index 1a49f66b8d7..d8ca55e4fed 100644 --- a/src/Distributed/VirtualDataContainer.f90 +++ b/src/Distributed/VirtualDataContainer.f90 @@ -326,6 +326,9 @@ subroutine print_items(this, imon, items) vdi => get_virtual_data_from_list(this%virtual_data_list, items%at(i)) write (imon, *) vdi%var_name, ":", vdi%mem_path end do + if (items%size == 0) then + write (imon, *) "... empty ...", this%name + end if write (imon, *) "<===== items" end subroutine print_items diff --git a/src/Model/Connection/GwfGwfConnection.f90 b/src/Model/Connection/GwfGwfConnection.f90 index 625c41d8525..bbfdb6961dc 100644 --- a/src/Model/Connection/GwfGwfConnection.f90 +++ b/src/Model/Connection/GwfGwfConnection.f90 @@ -13,7 +13,7 @@ module GwfGwfConnectionModule use DisConnExchangeModule use GwfGwfExchangeModule, only: GwfExchangeType, GetGwfExchangeFromList, & CastAsGwfExchange - use GwfNpfModule, only: GwfNpfType, hcond, vcond + use GwfNpfModule, only: GwfNpfType use GwfBuyModule, only: GwfBuyType use BaseDisModule, only: DisBaseType use ConnectionsModule, only: ConnectionsType diff --git a/src/Model/GroundWaterFlow/gwf3lak8.f90 b/src/Model/GroundWaterFlow/gwf3lak8.f90 index 3834f640263..bda281b18a6 100644 --- a/src/Model/GroundWaterFlow/gwf3lak8.f90 +++ b/src/Model/GroundWaterFlow/gwf3lak8.f90 @@ -829,7 +829,7 @@ subroutine lak_read_lake_connections(this) 'DNODATA (', DNODATA, ') value.' ! ! -- create deprecation warning - call deprecation_warning('CONNECTIONDATA', 'bedleak=NONE', '6.5.0', & + call deprecation_warning('CONNECTIONDATA', 'bedleak=NONE', '6.4.3', & warnmsg, this%parser%GetUnit()) case default read (keyword, *) rval diff --git a/src/Model/GroundWaterFlow/gwf3mvr8.f90 b/src/Model/GroundWaterFlow/gwf3mvr8.f90 index 5488b31e0e9..cf536f3c4e1 100644 --- a/src/Model/GroundWaterFlow/gwf3mvr8.f90 +++ b/src/Model/GroundWaterFlow/gwf3mvr8.f90 @@ -102,7 +102,7 @@ ! if(this%inmvr > 0) call this%mvr%mvr_ot() ! module GwfMvrModule - use KindModule, only: DP, I4B + use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LENMEMPATH, LENPACKAGENAME, LENMODELNAME, & LENBUDTXT, LENAUXNAME, LENPAKLOC, & DZERO, DNODATA, MAXCHARLEN, TABCENTER, & @@ -123,6 +123,7 @@ module GwfMvrModule public :: GwfMvrType, mvr_cr type, extends(NumericalPackageType) :: GwfMvrType + logical(LGP), pointer :: reset_mapped_id ! flag to indicate mapped ids must be reset; true when movers change integer(I4B), pointer :: ibudgetout => null() !< binary budget output file integer(I4B), pointer :: ibudcsv => null() !< unit number for csv budget output file integer(I4B), pointer :: maxmvr => null() !< max number of movers to be specified @@ -172,6 +173,7 @@ module GwfMvrModule procedure, private :: mvr_setup_budobj procedure, private :: mvr_setup_outputtab procedure, private :: mvr_print_outputtab + procedure, private :: set_mapped_id end type GwfMvrType @@ -337,6 +339,7 @@ subroutine mvr_rp(this) write (this%iout, '(/,2x,a,i0)') 'READING WATER MOVERS FOR PERIOD ', kper nlist = -1 i = 1 + this%reset_mapped_id = .true. ! ! -- set mname to '' if this is an exchange mover, or to the model name if (this%iexgmvr == 0) then @@ -492,19 +495,16 @@ subroutine mvr_bd(this) ! -- dummy class(GwfMvrType) :: this ! -- locals - integer(I4B) :: i, mapped_id - class(PackageMoverType), pointer :: pkg_mvr ! -- formats -! ------------------------------------------------------------------------------ ! - ! -- set the feature maps - allocate (pkg_mvr) - do i = 1, this%nmvr - call set_packagemover_pointer(pkg_mvr, this%mvr(i)%mem_path_src) - mapped_id = pkg_mvr%iprmap(this%mvr(i)%iRchNrSrc) - this%mvr(i)%iRchNrSrcMapped = mapped_id - end do - deallocate (pkg_mvr) + ! -- set the feature maps; for performance reasons, + ! this should only be called for the first time + ! step of a stress period in which a new set of + ! movers was provided in a period block. + if (this%reset_mapped_id) then + call this%set_mapped_id() + this%reset_mapped_id = .false. + end if ! ! -- fill the budget object call this%fill_budobj() @@ -699,6 +699,7 @@ subroutine mvr_da(this) end if ! ! -- Scalars + call mem_deallocate(this%reset_mapped_id) call mem_deallocate(this%ibudgetout) call mem_deallocate(this%ibudcsv) call mem_deallocate(this%maxmvr) @@ -1045,6 +1046,7 @@ subroutine allocate_scalars(this) call this%NumericalPackageType%allocate_scalars() ! ! -- Allocate + call mem_allocate(this%reset_mapped_id, 'RESET_MAPPED_ID', this%memoryPath) call mem_allocate(this%ibudgetout, 'IBUDGETOUT', this%memoryPath) call mem_allocate(this%ibudcsv, 'IBUDCSV', this%memoryPath) call mem_allocate(this%maxmvr, 'MAXMVR', this%memoryPath) @@ -1055,6 +1057,7 @@ subroutine allocate_scalars(this) call mem_allocate(this%imodelnames, 'IMODELNAMES', this%memoryPath) ! ! -- Initialize + this%reset_mapped_id = .false. this%ibudgetout = 0 this%ibudcsv = 0 this%maxmvr = -1 @@ -1345,4 +1348,41 @@ subroutine mvr_print_outputtab(this) return end subroutine mvr_print_outputtab + !> @brief Set mapped id + !! + !! For the budget output, we don't write outlet number, + !! instead we write the lake number. Normally the receiver + !! number is the same as the feature number provided by the + !! user. For moving water from a lake, the user specifies the + !! outlet number, not the lake number, in the mover package. + !! The iRchNrSrcMapped variable contains the lake number, not + !! the outlet number, and is written to the budget files. For + !! other packages, the iRchNrSrcMapped value is simply the well + !! number, the stream reach, or the uzf cell number. + !! This routine needs to be called each time a new set of movers + !! is read. It can't be called from within mvr_rp because the + !! iprmap isn't updated by lake until lak_rp, which is called + !! after mvr_rp. + !< + subroutine set_mapped_id(this) + ! -- dummy + class(GwfMvrType) :: this + ! -- locals + integer(I4B) :: i, mapped_id + class(PackageMoverType), pointer :: pkg_mvr + ! -- formats + ! + ! -- set the feature maps + allocate (pkg_mvr) + do i = 1, this%nmvr + call set_packagemover_pointer(pkg_mvr, this%mvr(i)%mem_path_src) + mapped_id = pkg_mvr%iprmap(this%mvr(i)%iRchNrSrc) + this%mvr(i)%iRchNrSrcMapped = mapped_id + end do + deallocate (pkg_mvr) + ! + ! -- Return + return + end subroutine set_mapped_id + end module diff --git a/src/Model/GroundWaterFlow/gwf3sfr8.f90 b/src/Model/GroundWaterFlow/gwf3sfr8.f90 index c9a4141fc9d..dea9390e846 100644 --- a/src/Model/GroundWaterFlow/gwf3sfr8.f90 +++ b/src/Model/GroundWaterFlow/gwf3sfr8.f90 @@ -923,7 +923,7 @@ subroutine sfr_read_packagedata(this) 'of 0 0 0 should be specified for unconnected reaches' ! ! -- create deprecation warning - call deprecation_warning('PACKAGEDATA', 'CELLID=NONE', '6.5.0', & + call deprecation_warning('PACKAGEDATA', 'CELLID=NONE', '6.4.3', & warnmsg, this%parser%GetUnit()) else diff --git a/src/Solution/PETSc/PetscConvergence.F90 b/src/Solution/PETSc/PetscConvergence.F90 index 7c4dcd2cf81..0928d7aa017 100644 --- a/src/Solution/PETSc/PetscConvergence.F90 +++ b/src/Solution/PETSc/PetscConvergence.F90 @@ -92,10 +92,11 @@ subroutine petsc_check_convergence(ksp, n, rnorm, flag, context, ierr) CHKERRQ(ierr) call VecCopy(res, context%res_old, ierr) CHKERRQ(ierr) - call VecDestroy(res, ierr) - CHKERRQ(ierr) flag = KSP_CONVERGED_ITERATING end if + + call VecDestroy(res, ierr) + CHKERRQ(ierr) return end if @@ -131,10 +132,10 @@ subroutine petsc_check_convergence(ksp, n, rnorm, flag, context, ierr) call VecDestroy(res, ierr) CHKERRQ(ierr) - ! get dv and dr per local model - call VecGetArrayF90(context%delta_x, local_dx, ierr) + ! get dv and dr per local model (readonly!) + call VecGetArrayReadF90(context%delta_x, local_dx, ierr) CHKERRQ(ierr) - call VecGetArrayF90(context%delta_res, local_dr, ierr) + call VecGetArrayReadF90(context%delta_res, local_dr, ierr) CHKERRQ(ierr) do i = 1, summary%convnmod ! reset @@ -162,9 +163,9 @@ subroutine petsc_check_convergence(ksp, n, rnorm, flag, context, ierr) summary%convlocdr(i, iter_cnt) = idx_dr end if end do - call VecRestoreArrayF90(x, local_dx, ierr) + call VecRestoreArrayF90(context%delta_x, local_dx, ierr) CHKERRQ(ierr) - call VecRestoreArrayF90(x, local_dr, ierr) + call VecRestoreArrayF90(context%delta_res, local_dr, ierr) CHKERRQ(ierr) if (norm < context%dvclose) then diff --git a/src/Solution/ParallelSolution.f90 b/src/Solution/ParallelSolution.f90 index 687755e71ae..09c00392ef9 100644 --- a/src/Solution/ParallelSolution.f90 +++ b/src/Solution/ParallelSolution.f90 @@ -43,6 +43,7 @@ function par_has_converged(this, max_dvc) result(has_converged) abs_max_dvc = abs(max_dvc) call MPI_Allreduce(abs_max_dvc, global_max_dvc, 1, MPI_DOUBLE_PRECISION, & MPI_MAX, mpi_world%comm, ierr) + call CHECK_MPI(ierr) if (global_max_dvc <= this%dvclose) then has_converged = .true. end if @@ -68,6 +69,7 @@ function par_package_convergence(this, dpak, cpakout, iend) & call MPI_Allreduce(icnvg_local, icnvg_global, 1, MPI_INTEGER, & MPI_MIN, mpi_world%comm, ierr) + call CHECK_MPI(ierr) end function par_package_convergence @@ -82,6 +84,7 @@ function par_sync_newtonur_flag(this, inewtonur) result(ivalue) mpi_world => get_mpi_world() call MPI_Allreduce(inewtonur, ivalue, 1, MPI_INTEGER, & MPI_MAX, mpi_world%comm, ierr) + call CHECK_MPI(ierr) end function par_sync_newtonur_flag @@ -108,6 +111,7 @@ function par_nur_has_converged(this, dxold_max, hncg) & call MPI_Allreduce(icnvg_local, icnvg_global, 1, MPI_INTEGER, & MPI_MIN, mpi_world%comm, ierr) + call CHECK_MPI(ierr) if (icnvg_global == 1) has_converged = .true. end function par_nur_has_converged @@ -131,6 +135,7 @@ subroutine par_calc_ptc(this, iptc, ptcf) ! now reduce call MPI_Allreduce(ptcf_loc, ptcf_glo_max, 1, MPI_DOUBLE_PRECISION, & MPI_MAX, mpi_world%comm, ierr) + call CHECK_MPI(ierr) iptc = 0 ptcf = DZERO @@ -161,8 +166,10 @@ subroutine par_underrelax(this, kiter, bigch, neq, active, x, xtemp) ! first reduce largest change over all processes call MPI_Allreduce(bigch, dvc_global_max, 1, MPI_DOUBLE_PRECISION, & MPI_MAX, mpi_world%comm, ierr) + call CHECK_MPI(ierr) call MPI_Allreduce(bigch, dvc_global_min, 1, MPI_DOUBLE_PRECISION, & MPI_MIN, mpi_world%comm, ierr) + call CHECK_MPI(ierr) if (abs(dvc_global_min) > abs(dvc_global_max)) then dvc_global_max = dvc_global_min @@ -192,6 +199,7 @@ subroutine par_backtracking_xupdate(this, btflag) call MPI_Allreduce(btflag_local, btflag, 1, MPI_INTEGER, & MPI_MAX, mpi_world%comm, ierr) + call CHECK_MPI(ierr) end subroutine par_backtracking_xupdate diff --git a/src/Utilities/STLVecInt.f90 b/src/Utilities/STLVecInt.f90 index 6551e697705..7431da49026 100644 --- a/src/Utilities/STLVecInt.f90 +++ b/src/Utilities/STLVecInt.f90 @@ -21,10 +21,12 @@ module STLVecIntModule procedure, pass(this) :: add_array_unique !< adds elements of array at the end of the vector, if not present yet procedure, pass(this) :: at !< random access, unsafe, no bounds checking procedure, pass(this) :: at_safe !< random access with bounds checking + procedure, pass(this) :: set !< set value at index, no bounds checking procedure, pass(this) :: clear !< empties the vector, leaves memory unchanged procedure, pass(this) :: shrink_to_fit !< reduces the allocated memory to fit the actual vector size procedure, pass(this) :: destroy !< deletes the memory procedure, pass(this) :: contains !< true when element already present + procedure, pass(this) :: get_index !< return index of first occurrence of value in array, -1 when not present procedure, pass(this) :: get_values !< returns a copy of the values ! private procedure, private, pass(this) :: expand @@ -119,6 +121,15 @@ function at_safe(this, idx) result(value) end function at_safe + subroutine set(this, idx, value) + class(STLVecInt), intent(inout) :: this + integer(I4B), intent(in) :: idx + integer(I4B) :: value + + this%values(idx) = value + + end subroutine set + subroutine clear(this) class(STLVecInt), intent(inout) :: this @@ -200,6 +211,25 @@ function contains(this, val) result(res) end function contains + !> @brief Return index of first occurrence, + !< returns -1 when not present + function get_index(this, val) result(idx) + class(STLVecInt), intent(inout) :: this + integer(I4B) :: val + integer(I4B) :: idx + ! local + integer(I4B) :: i + + idx = -1 + do i = 1, this%size + if (this%at(i) == val) then + idx = i + return + end if + end do + + end function get_index + function get_values(this) result(values) class(STLVecInt), intent(in) :: this integer(I4B), dimension(:), allocatable :: values diff --git a/src/Utilities/SimStages.f90 b/src/Utilities/SimStages.f90 index 9168f05b9a8..01bd7d1e312 100644 --- a/src/Utilities/SimStages.f90 +++ b/src/Utilities/SimStages.f90 @@ -21,6 +21,7 @@ module SimStagesModule integer(I4B), public, parameter :: STG_BFR_EXG_AD = 12 !< before exchange advance (per solution) integer(I4B), public, parameter :: STG_BFR_EXG_CF = 13 !< before exchange calculate (per solution) integer(I4B), public, parameter :: STG_BFR_EXG_FC = 14 !< before exchange formulate (per solution) + integer(I4B), public, parameter :: NR_SIM_STAGES = 14 !< before exchange formulate (per solution) contains diff --git a/src/Utilities/version.f90 b/src/Utilities/version.f90 index 011f2ec6917..809c8f5e3f4 100644 --- a/src/Utilities/version.f90 +++ b/src/Utilities/version.f90 @@ -17,8 +17,8 @@ module VersionModule public ! -- modflow 6 version integer(I4B), parameter :: IDEVELOPMODE = 1 - character(len=*), parameter :: VERSIONNUMBER = '6.4.3' - character(len=*), parameter :: VERSIONTAG = ' 02/07/2024' + character(len=*), parameter :: VERSIONNUMBER = '6.4.4' + character(len=*), parameter :: VERSIONTAG = ' 02/13/2024' character(len=40), parameter :: VERSION = VERSIONNUMBER//VERSIONTAG character(len=2), parameter :: MFVNAM = ' 6' character(len=*), parameter :: MFTITLE = & diff --git a/src/meson.build b/src/meson.build index b25a391edb3..8bd2cb66d93 100644 --- a/src/meson.build +++ b/src/meson.build @@ -275,6 +275,8 @@ modflow_petsc_sources = files( modflow_mpi_sources = files( 'Distributed' / 'MpiMessageBuilder.f90', + 'Distributed' / 'MpiMessageCache.f90', + 'Distributed' / 'MpiUnitCache.f90', 'Distributed' / 'MpiRouter.f90', 'Distributed' / 'MpiRunControl.F90', 'Distributed' / 'MpiWorld.f90', diff --git a/utils/mf5to6/make/makefile b/utils/mf5to6/make/makefile index ccf5741d638..0e3b414c298 100644 --- a/utils/mf5to6/make/makefile +++ b/utils/mf5to6/make/makefile @@ -5,10 +5,10 @@ include ./makedefaults # Define the source file directories SOURCEDIR1=../src -SOURCEDIR2=../src/LGR -SOURCEDIR3=../src/Preproc -SOURCEDIR4=../src/MF2005 -SOURCEDIR5=../src/NWT +SOURCEDIR2=../src/MF2005 +SOURCEDIR3=../src/NWT +SOURCEDIR4=../src/LGR +SOURCEDIR5=../src/Preproc SOURCEDIR6=../../../src/Utilities/Memory SOURCEDIR7=../../../src/Utilities/TimeSeries SOURCEDIR8=../../../src/Utilities diff --git a/version.txt b/version.txt index 8e0a2119a7d..8b38d69c522 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -6.4.3 \ No newline at end of file +6.4.4 \ No newline at end of file