From e60ab0fe83e832433617d849d720b44fe2a3597a Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Mon, 24 Jul 2023 09:53:14 -0400 Subject: [PATCH 01/13] docs: describe building and testing makefiles in DEVELOPER.md (#1307) Co-authored-by: Eric Morway Co-authored-by: langevin-usgs --- DEVELOPER.md | 78 ++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 66 insertions(+), 12 deletions(-) diff --git a/DEVELOPER.md b/DEVELOPER.md index 4b30e2f8174..1cfbba5524a 100644 --- a/DEVELOPER.md +++ b/DEVELOPER.md @@ -7,7 +7,6 @@ To build and test a parallel version of the program, first read the instructions - - [Prerequisites](#prerequisites) - [Git](#git) - [Fortran compiler](#fortran-compiler) @@ -43,11 +42,16 @@ To build and test a parallel version of the program, first read the instructions - [Selecting tests with markers](#selecting-tests-with-markers) - [External model tests](#external-model-tests) - [Writing tests](#writing-tests) +- [Generating makefiles](#generating-makefiles) + - [Updating extra and excluded files](#updating-extra-and-excluded-files) + - [Testing makefiles](#testing-makefiles) + - [Installing `make` on Windows](#installing-make-on-windows) + - [Using Conda from Git Bash](#using-conda-from-git-bash) - [Git Strategy for Managing Long-Lived Branches](#git-strategy-for-managing-long-lived-branches) - - [Create a Backup](#create-a-backup) - - [Squash Feature Branch Commits](#squash-feature-branch-commits) - - [Rebase Feature Branch with Develop](#rebase-feature-branch-with-develop) - - [Cleanup](#cleanup) + - [Create a Backup](#create-a-backup) + - [Squash Feature Branch Commits](#squash-feature-branch-commits) + - [Rebase Feature Branch with Develop](#rebase-feature-branch-with-develop) + - [Cleanup](#cleanup) @@ -432,11 +436,59 @@ Tests should ideally follow a few conventions for easier maintenance: **Note:** If all three external model repositories are not installed as described above, some tests will be skipped. The full test suite includes >750 cases. All must pass before changes can be merged into this repository. +## Generating makefiles + +Run `build_makefiles.py` in the `distribution/` directory after adding, removing, or renaming source files. This script uses [Pymake](https://github.com/modflowpy/pymake) to regenerate makefiles. For instance: + +```shell +python build_makefiles.py +``` + +### Updating extra and excluded files + +If the utilities located in the `utils` directory (e.g., `mf5to6` and `zbud6`) are affected by changes to the modflow6 `src/` directory (such as new or refactored source files), then the new module source file should also be added to the utility's `utils//pymake/extrafiles.txt` file. This file informs Pymake of source files living outside the main source directory, so they can be included in generated makefiles. + +Module dependencies for features still under development should be added to `excludefiles.txt`. Source files listed in this file will be excluded from makefiles generated by Pymake. Makefiles should only include the source files needed to the build officially released/supported features. + +### Testing makefiles + +Makefile generation and usage can be tested from the `distribution` directory by running the `build_makefiles.py` script with Pytest: + +```shell +pytest -v build_makefiles.py +``` + +**Note**: `make` is required to test compiling MODFLOW 6 with makefiles. If `make` is not discovered on the system path, compile tests will be skipped. + +Makefiles may also be tested manually by changing to the appropriate `make` subdirectory (of the project root for MODFLOW 6, or inside the corresponding `utils` subdirectory for the zonebudget or converter utilities) and invoking `make` (`make clean` may first be necessary to remove previously created object files). + +### Installing `make` on Windows + +On Windows, it is recommended to generate and test makefiles from a Unix-like shell rather than PowerShell or Command Prompt. Make can be installed via [Conda](https://anaconda.org/conda-forge/make) or [Chocolatey](https://community.chocolatey.org/packages/make). Alternatively, it is included with [mingw](https://sourceforge.net/projects/mingw/), which is also available from [Chocolatey](https://community.chocolatey.org/packages/mingw). + +#### Using Conda from Git Bash + +To use Conda from Git Bash on Windows, first run the `conda.sh` script located in your Conda installation's `/etc/profile.d` subdirectory. For instance, with Anaconda3: + +```shell +. /c/Anaconda3/etc/profile.d/conda.sh +``` + +Or Miniconda3: + +```shell +. /c/ProgramData/miniconda3/etc/profile.d/conda.sh +``` + +After this, `conda` commands should be available. + +This command may be added to a `.bashrc` or `.bash_profile` file in your home directory to permanently configure Git Bash for Conda. + ## Git Strategy for Managing Long-Lived Branches When a feature branch takes a long time to develop, it is easy to become out of sync with the develop branch. Depending on the situation, it may be advisable to periodically squash the commits on the feature branch and rebase the change set with develop. The following approach for updating a long-lived feature branch has proven robust. -In the example below, the feature branch is assumed to be called `feat-xyz` +In the example below, the feature branch is assumed to be called `feat-xyz`. ### Create a Backup @@ -452,22 +504,24 @@ git checkout feat-xyz Next, consider squashing commits on the feature branch. If there are many commits, it is beneficial to squash them before trying to rebase with develop. There is a nice article on [squashing commits into one using git](https://www.internalpointers.com/post/squash-commits-into-one-git), which has been very useful for consolidating commits on a long-lived modflow6 feature branch. -A quick and dirty way to squash without interactive rebase (as an alternative to the approach described in the article mentioned in the preceding paragraph) is a soft reset followed by an ammended commit. The following commands will also squash the feature branch commits, provided a backup of the branch has been made as described above (use caution; accidentally typing `--hard` will wipe out all your work). +A quick and dirty way to squash without interactive rebase (as an alternative to the approach described in the article mentioned in the preceding paragraph) is a soft reset followed by an ammended commit. First making a backup of the feature branch is strongly recommended before using this approach, as accidentally typing `--hard` instead of `--soft` will wipe out all your work. ``` -git reset --soft +git reset --soft git commit --amend -m "consolidated commit message" ``` -Once the commits on the feature branch have been consolidated, a force push to origin is recommended. The force push to origin is not strictly required, but it can serve an intermediate backup/checkpoint so the squashed branch state can be retrieved if rebasing fails. The following command will push `feat-xyz` to origin. +Once the commits on the feature branch have been consolidated, a force push to origin is recommended. This is not strictly required, but it can serve as an intermediate backup/checkpoint so the squashed branch state can be retrieved if rebasing fails. The following command will push `feat-xyz` to origin. ``` git push origin feat-xyz --force ``` +The `--force` flag's short form is `-f`. + ### Rebase Feature Branch with Develop -Now that the commits on `feat-xyz` have been consolidated, it is time to rebase with develop. If there are multiple commits in `feat-xyz` that make changes and undo them or rename files and move things around in subsequent commits, then there may be multiple sets of merge conflicts that will need to be resolved as the rebase works its way through the commit change sets. This is why it is beneficial to squash the feature commits before rebasing with develop. +Now that the commits on `feat-xyz` have been consolidated, it is time to rebase with develop. If there are multiple commits in `feat-xyz` that make changes, undo them, rename files, and/or move things around in subsequent commits, then there may be multiple sets of merge conflicts that will need to be resolved as the rebase works its way through the commit change sets. This is why it is beneficial to squash the feature commits before rebasing with develop. To rebase with develop, make sure the feature branch is checked out and then type: @@ -475,7 +529,7 @@ To rebase with develop, make sure the feature branch is checked out and then typ git rebase develop ``` -If anything goes wrong during a rebase, there is the `rebase --abort` command to unwind a rebase that's gotten out of hand. +If anything goes wrong during a rebase, there is the `rebase --abort` command to unwind it. If there are merge conflicts, they will need to be resolved before going forward. Once any conflicts are resolved, it may be worthwhile to rebuild the MODFLOW 6 program and run the smoke tests to ensure nothing is broken. @@ -493,4 +547,4 @@ Lastly, if you are satisfied with the results and confident the procedure went w git branch -d feat-xyz-backup ``` -This process can be repeated periodically to stay in sync with the develop branch and keep a clean commit history. \ No newline at end of file +This process can be repeated periodically to stay in sync with the develop branch and keep a clean commit history. From d6b4350d0ba7ab54fb8747d3cc76d12af1d04126 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Wed, 26 Jul 2023 11:20:33 -0500 Subject: [PATCH 02/13] fix(backspace): remove Fortran backspace (#1310) * remove backspace so that ifx can be used to compile * update meson for ifx * fix initialization error in mf5to6 ghost node writer * introduce LongLineReader as a way to emulate backspace functionality * update vscode tasks.json for consistent build group definition --- .vscode/tasks.json | 40 ++- autotest/test_cli.py | 9 +- make/makedefaults | 2 +- make/makefile | 35 +-- meson.build | 18 +- msvs/mf6core.vfproj | 127 +++++---- src/Model/GroundWaterFlow/gwf3evt8.f90 | 3 +- src/Model/GroundWaterFlow/gwf3rch8.f90 | 3 +- src/Model/ModelUtilities/BoundaryPackage.f90 | 3 +- .../ModelUtilities/DiscretizationBase.f90 | 13 +- src/Utilities/ArrayReaders.f90 | 4 +- src/Utilities/BlockParser.f90 | 257 +++++++++++++++++- src/Utilities/InputOutput.f90 | 244 +---------------- src/Utilities/ListReader.f90 | 31 ++- src/Utilities/LongLineReader.f90 | 117 ++++++++ src/meson.build | 1 + utils/mf5to6/make/makedefaults | 2 +- utils/mf5to6/make/makefile | 3 +- utils/mf5to6/msvs/mf5to6.vfproj | 54 ++-- utils/mf5to6/pymake/extrafiles.txt | 1 + utils/mf5to6/src/Connection.f90 | 30 ++ utils/mf5to6/src/Preproc/ObsBlock.f90 | 2 +- utils/mf5to6/src/Preproc/Preproc.f90 | 7 +- utils/mf5to6/src/Preproc/Utilities.f90 | 64 +---- utils/mf5to6/src/mf5to6.f90 | 3 +- utils/zonebudget/make/makedefaults | 2 +- utils/zonebudget/make/makefile | 3 +- utils/zonebudget/msvs/zonebudget.vfproj | 53 ++-- utils/zonebudget/pymake/extrafiles.txt | 1 + 29 files changed, 656 insertions(+), 476 deletions(-) create mode 100644 src/Utilities/LongLineReader.f90 diff --git a/.vscode/tasks.json b/.vscode/tasks.json index 490baff6416..76f6cb2fe2d 100644 --- a/.vscode/tasks.json +++ b/.vscode/tasks.json @@ -26,7 +26,10 @@ "release", "build", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -53,7 +56,10 @@ "release", "rebuild", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -81,7 +87,10 @@ "release", "build", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -108,7 +117,10 @@ "release", "rebuild", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -136,7 +148,10 @@ "debug", "build", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -163,7 +178,10 @@ "debug", "rebuild", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -191,7 +209,10 @@ "debug", "build", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } @@ -218,7 +239,10 @@ "debug", "rebuild", ], - "group": "build", + "group": { + "kind": "build", + "isDefault": true + }, "presentation": { "clear": true } diff --git a/autotest/test_cli.py b/autotest/test_cli.py index 6aa73052988..d1f033bb0e9 100644 --- a/autotest/test_cli.py +++ b/autotest/test_cli.py @@ -1,16 +1,21 @@ import subprocess +import platform from conftest import project_root_path bin_path = project_root_path / "bin" +app = "mf6" +ext = ".exe" if platform.system() == "Windows" else "" +app = f"{app}{ext}" + def test_cli_version(): output = " ".join( - subprocess.check_output([str(bin_path / "mf6"), "-v"]).decode().split() + subprocess.check_output([str(bin_path / app), "-v"]).decode().split() ) print(output) - assert output.startswith("mf6:") + assert output.startswith(f"{app}:"), f"found: {output}" version = ( output.lower().split(' ')[1] diff --git a/make/makedefaults b/make/makedefaults index 7bbd4dd3293..a9174e9746d 100644 --- a/make/makedefaults +++ b/make/makedefaults @@ -1,4 +1,4 @@ -# makedefaults created by pymake (version 1.2.7) for the 'mf6' executable. +# makedefaults created by pymake (version 1.2.9.dev0) for the 'mf6' executable. # determine OS ifeq ($(OS), Windows_NT) diff --git a/make/makefile b/make/makefile index 552d7c70a1d..1fc009719ea 100644 --- a/make/makefile +++ b/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.7) for the 'mf6' executable. +# makefile created by pymake (version 1.2.9.dev0) for the 'mf6' executable. include ./makedefaults @@ -88,6 +88,7 @@ $(OBJDIR)/MemoryHelper.o \ $(OBJDIR)/CharString.o \ $(OBJDIR)/Memory.o \ $(OBJDIR)/List.o \ +$(OBJDIR)/LongLineReader.o \ $(OBJDIR)/MemoryList.o \ $(OBJDIR)/TimeSeriesRecord.o \ $(OBJDIR)/BlockParser.o \ @@ -122,28 +123,21 @@ $(OBJDIR)/PackageMover.o \ $(OBJDIR)/Obs3.o \ $(OBJDIR)/NumericalPackage.o \ $(OBJDIR)/Budget.o \ -$(OBJDIR)/SeqVector.o \ $(OBJDIR)/sort.o \ $(OBJDIR)/SfrCrossSectionUtils.o \ $(OBJDIR)/BudgetTerm.o \ -$(OBJDIR)/BoundaryPackage.o \ -$(OBJDIR)/BaseModel.o \ -$(OBJDIR)/SparseMatrix.o \ -$(OBJDIR)/LinearSolverBase.o \ -$(OBJDIR)/ims8reordering.o \ $(OBJDIR)/VirtualBase.o \ $(OBJDIR)/STLVecInt.o \ +$(OBJDIR)/BoundaryPackage.o \ +$(OBJDIR)/BaseModel.o \ $(OBJDIR)/InputDefinition.o \ $(OBJDIR)/SfrCrossSectionManager.o \ $(OBJDIR)/dag_module.o \ $(OBJDIR)/BudgetObject.o \ -$(OBJDIR)/NumericalModel.o \ -$(OBJDIR)/BaseExchange.o \ -$(OBJDIR)/ImsLinearSolver.o \ -$(OBJDIR)/ims8base.o \ $(OBJDIR)/VirtualDataLists.o \ $(OBJDIR)/VirtualDataContainer.o \ $(OBJDIR)/SimStages.o \ +$(OBJDIR)/NumericalModel.o \ $(OBJDIR)/simnamidm.o \ $(OBJDIR)/gwt1idm.o \ $(OBJDIR)/gwt1dsp1idm.o \ @@ -168,13 +162,9 @@ $(OBJDIR)/gwf3lak8.o \ $(OBJDIR)/GwfVscInputData.o \ $(OBJDIR)/gwf3ghb8.o \ $(OBJDIR)/gwf3drn8.o \ -$(OBJDIR)/Timer.o \ -$(OBJDIR)/NumericalExchange.o \ -$(OBJDIR)/LinearSolverFactory.o \ -$(OBJDIR)/ims8linear.o \ -$(OBJDIR)/BaseSolution.o \ $(OBJDIR)/IndexMap.o \ $(OBJDIR)/VirtualModel.o \ +$(OBJDIR)/BaseExchange.o \ $(OBJDIR)/IdmSimDfnSelector.o \ $(OBJDIR)/IdmGwtDfnSelector.o \ $(OBJDIR)/IdmGwfDfnSelector.o \ @@ -187,9 +177,10 @@ $(OBJDIR)/gwf3tvk8.o \ $(OBJDIR)/MemoryManagerExt.o \ $(OBJDIR)/gwf3vsc8.o \ $(OBJDIR)/GwfNpfOptions.o \ -$(OBJDIR)/NumericalSolution.o \ $(OBJDIR)/InterfaceMap.o \ +$(OBJDIR)/SeqVector.o \ $(OBJDIR)/CellWithNbrs.o \ +$(OBJDIR)/NumericalExchange.o \ $(OBJDIR)/IdmDfnSelector.o \ $(OBJDIR)/gwf3uzf8.o \ $(OBJDIR)/gwt1apt1.o \ @@ -207,6 +198,9 @@ $(OBJDIR)/GwfMvrPeriodData.o \ $(OBJDIR)/ims8misc.o \ $(OBJDIR)/GwfBuyInputData.o \ $(OBJDIR)/VirtualSolution.o \ +$(OBJDIR)/SparseMatrix.o \ +$(OBJDIR)/LinearSolverBase.o \ +$(OBJDIR)/ims8reordering.o \ $(OBJDIR)/ArrayReaderBase.o \ $(OBJDIR)/VirtualExchange.o \ $(OBJDIR)/gwf3disu8.o \ @@ -244,12 +238,18 @@ $(OBJDIR)/GhostNode.o \ $(OBJDIR)/gwf3evt8.o \ $(OBJDIR)/gwf3chd8.o \ $(OBJDIR)/RouterBase.o \ +$(OBJDIR)/ImsLinearSolver.o \ +$(OBJDIR)/ims8base.o \ $(OBJDIR)/Integer2dReader.o \ $(OBJDIR)/GridConnection.o \ $(OBJDIR)/DistributedVariable.o \ $(OBJDIR)/gwt1.o \ $(OBJDIR)/gwf3.o \ $(OBJDIR)/SerialRouter.o \ +$(OBJDIR)/Timer.o \ +$(OBJDIR)/LinearSolverFactory.o \ +$(OBJDIR)/ims8linear.o \ +$(OBJDIR)/BaseSolution.o \ $(OBJDIR)/StructVector.o \ $(OBJDIR)/IdmLogger.o \ $(OBJDIR)/Integer1dReader.o \ @@ -262,6 +262,7 @@ $(OBJDIR)/GwtGwtExchange.o \ $(OBJDIR)/GwfInterfaceModel.o \ $(OBJDIR)/GwfGwfExchange.o \ $(OBJDIR)/RouterFactory.o \ +$(OBJDIR)/NumericalSolution.o \ $(OBJDIR)/MappedMemory.o \ $(OBJDIR)/StructArray.o \ $(OBJDIR)/ModflowInput.o \ diff --git a/meson.build b/meson.build index 879bf5d4c3a..adc2078da00 100644 --- a/meson.build +++ b/meson.build @@ -25,6 +25,7 @@ message('The used profile is:', profile) # parse compiler options fc = meson.get_compiler('fortran') fc_id = fc.get_id() +message('The fc_id is:', fc_id) compile_args = [] link_args = [] @@ -85,8 +86,23 @@ elif fc_id == 'intel' '-diag-disable:5268', # Line too long ] link_args += '-static-intel' -endif + +# Command line options for ifx +elif fc_id == 'intel-llvm-cl' + # windows + compile_args += ['/fpe:0', # Activate all floating point exceptions + '/heap-arrays:0', + '/traceback', + '/fpp', # Activate preprocessing + '/Qdiag-disable:7416', # f2008 warning + '/Qdiag-disable:7025', # f2008 warning + '/Qdiag-disable:5268', # Line too long + ] + link_args += ['/ignore:4217', # access through ddlimport might be inefficient + '/ignore:4286' # same as 4217, but more general + ] +endif # parallel build options is_parallel_build = get_option('parallel') diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index e2be9f2e7b6..b2fe810728e 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -2,44 +2,50 @@ - + + - - - - - - - - + + + + + + + + + - - - - - - - - + + + + + + + + + - - - - - - - - + + + + + + + + + - - - - - - - - + + + + + + + + + + @@ -72,13 +78,13 @@ - + - + - + - + @@ -224,13 +230,13 @@ - + - + - + - + @@ -239,13 +245,13 @@ - + - + - + - + @@ -336,29 +342,30 @@ - + - + - + - + - + - + - + - + + @@ -381,13 +388,15 @@ - + - + - + - + - - + + + + diff --git a/src/Model/GroundWaterFlow/gwf3evt8.f90 b/src/Model/GroundWaterFlow/gwf3evt8.f90 index 34d687cfbdf..37652b31151 100644 --- a/src/Model/GroundWaterFlow/gwf3evt8.f90 +++ b/src/Model/GroundWaterFlow/gwf3evt8.f90 @@ -1093,7 +1093,8 @@ subroutine evt_rp_list(this, inrate) ! nlist = -1 maxboundorig = this%maxbound - call this%dis%read_list(this%parser%iuactive, this%iout, this%iprpak, & + call this%dis%read_list(this%parser%line_reader, & + this%parser%iuactive, this%iout, this%iprpak, & nlist, this%inamedbound, this%iauxmultcol, & this%nodelist, this%bound, this%auxvar, & this%auxname, this%boundname, this%listlabel, & diff --git a/src/Model/GroundWaterFlow/gwf3rch8.f90 b/src/Model/GroundWaterFlow/gwf3rch8.f90 index e50ad4f279b..d0add5b9f79 100644 --- a/src/Model/GroundWaterFlow/gwf3rch8.f90 +++ b/src/Model/GroundWaterFlow/gwf3rch8.f90 @@ -591,7 +591,8 @@ subroutine rch_rp_list(this, inrech) ! ! -- read the list of recharge values; scale the recharge by auxmultcol ! if it is specified. - call this%dis%read_list(this%parser%iuactive, this%iout, this%iprpak, & + call this%dis%read_list(this%parser%line_reader, & + this%parser%iuactive, this%iout, this%iprpak, & nlist, this%inamedbound, this%iauxmultcol, & this%nodelist, this%bound, this%auxvar, & this%auxname, this%boundname, this%listlabel, & diff --git a/src/Model/ModelUtilities/BoundaryPackage.f90 b/src/Model/ModelUtilities/BoundaryPackage.f90 index 1ec95293e16..8fd71401f45 100644 --- a/src/Model/ModelUtilities/BoundaryPackage.f90 +++ b/src/Model/ModelUtilities/BoundaryPackage.f90 @@ -364,7 +364,8 @@ subroutine bnd_rp(this) call this%TasManager%Reset(this%packName) ! ! -- Read data as a list - call this%dis%read_list(this%parser%iuactive, this%iout, & + call this%dis%read_list(this%parser%line_reader, & + this%parser%iuactive, this%iout, & this%iprpak, nlist, this%inamedbound, & this%iauxmultcol, this%nodelist, & this%bound, this%auxvar, this%auxname, & diff --git a/src/Model/ModelUtilities/DiscretizationBase.f90 b/src/Model/ModelUtilities/DiscretizationBase.f90 index dd03ac8faa3..6c43a150d01 100644 --- a/src/Model/ModelUtilities/DiscretizationBase.f90 +++ b/src/Model/ModelUtilities/DiscretizationBase.f90 @@ -1016,9 +1016,9 @@ subroutine fill_dbl_array(this, buff1, buff2) return end subroutine fill_dbl_array - subroutine read_list(this, in, iout, iprpak, nlist, inamedbound, & - iauxmultcol, nodelist, rlist, auxvar, auxname, & - boundname, label, pkgname, tsManager, iscloc, & + subroutine read_list(this, line_reader, in, iout, iprpak, nlist, & + inamedbound, iauxmultcol, nodelist, rlist, auxvar, & + auxname, boundname, label, pkgname, tsManager, iscloc, & indxconvertflux) ! ****************************************************************************** ! read_list -- Read a list using the list reader object. @@ -1032,6 +1032,7 @@ subroutine read_list(this, in, iout, iprpak, nlist, inamedbound, & ! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LENBOUNDNAME, LINELENGTH + use LongLineReaderModule, only: LongLineReaderType use ListReaderModule, only: ListReaderType use SimModule, only: store_error, store_error_unit, count_errors use InputOutputModule, only: urword @@ -1039,6 +1040,7 @@ subroutine read_list(this, in, iout, iprpak, nlist, inamedbound, & use TimeSeriesManagerModule, only: read_value_or_time_series ! -- dummy class(DisBaseType) :: this + type(LongLineReaderType), intent(inout) :: line_reader integer(I4B), intent(in) :: in integer(I4B), intent(in) :: iout integer(I4B), intent(in) :: iprpak @@ -1070,8 +1072,9 @@ subroutine read_list(this, in, iout, iprpak, nlist, inamedbound, & ! ------------------------------------------------------------------------------ ! ! -- Read the list - call lstrdobj%read_list(in, iout, nlist, inamedbound, this%mshape, & - nodelist, rlist, auxvar, auxname, boundname, label) + call lstrdobj%read_list(line_reader, in, iout, nlist, inamedbound, & + this%mshape, nodelist, rlist, auxvar, auxname, & + boundname, label) ! ! -- Go through all locations where a text string was found instead of ! a double precision value and make time-series links to rlist diff --git a/src/Utilities/ArrayReaders.f90 b/src/Utilities/ArrayReaders.f90 index 37d0141b369..32e49130faa 100644 --- a/src/Utilities/ArrayReaders.f90 +++ b/src/Utilities/ArrayReaders.f90 @@ -674,6 +674,7 @@ subroutine read_control_dbl(iu, iout, aname, locat, cnstnt, & end subroutine read_control_dbl subroutine read_control_1(iu, iout, aname, locat, iclose, line, icol, fname) + use SimModule, only: ustop ! -- Read CONSTANT, INTERNAL, or OPEN/CLOSE from array control record. ! -- dummy integer(I4B), intent(in) :: iu @@ -690,7 +691,8 @@ subroutine read_control_1(iu, iout, aname, locat, iclose, line, icol, fname) integer(I4B) :: ierr real(DP) :: r ! - ! -- Read array control record. + ! -- Read array control record. Any future refactoring + ! should use the LongLineReader here instead of u9rdcom call u9rdcom(iu, iout, line, ierr) ! iclose = 0 diff --git a/src/Utilities/BlockParser.f90 b/src/Utilities/BlockParser.f90 index 2cc5159db58..9b6e015a3bb 100644 --- a/src/Utilities/BlockParser.f90 +++ b/src/Utilities/BlockParser.f90 @@ -7,17 +7,18 @@ module BlockParserModule use KindModule, only: DP, I4B - use ConstantsModule, only: LENHUGELINE, LINELENGTH, MAXCHARLEN + use ConstantsModule, only: LENBIGLINE, LENHUGELINE, LINELENGTH, MAXCHARLEN use VersionModule, only: IDEVELOPMODE - use InputOutputModule, only: uget_block, uget_any_block, uterminate_block, & - u9rdcom, urword, upcase + use InputOutputModule, only: urword, upcase, openfile, & + io_getunit => GetUnit use SimModule, only: store_error, store_error_unit use SimVariablesModule, only: errmsg + use LongLineReaderModule, only: LongLineReaderType implicit none private - public :: BlockParserType + public :: BlockParserType, uget_block, uget_any_block, uterminate_block type :: BlockParserType integer(I4B), public :: iuactive !< flag indicating if a file unit is active, variable is not used internally @@ -30,6 +31,7 @@ module BlockParserModule character(len=LINELENGTH), private :: blockNameFound !< block name found character(len=LENHUGELINE), private :: laststring !< last string read character(len=:), allocatable, private :: line !< current line + type(LongLineReaderType) :: line_reader contains procedure, public :: Initialize procedure, public :: Clear @@ -155,8 +157,9 @@ subroutine GetBlock(this, blockName, isFound, ierr, supportOpenClose, & this%blockNameFound = '' ! if (blockName == '*') then - call uget_any_block(this%inunit, this%iout, isFound, this%lloc, & - this%line, blockNameFound, this%iuext) + call uget_any_block(this%line_reader, this%inunit, this%iout, & + isFound, this%lloc, this%line, blockNameFound, & + this%iuext) if (isFound) then this%blockNameFound = blockNameFound ierr = 0 @@ -164,7 +167,8 @@ subroutine GetBlock(this, blockName, isFound, ierr, supportOpenClose, & ierr = 1 end if else - call uget_block(this%inunit, this%iout, this%blockName, ierr, isFound, & + call uget_block(this%line_reader, this%inunit, this%iout, & + this%blockName, ierr, isFound, & this%lloc, this%line, this%iuext, continueRead, & supportOpenCloseLocal) if (isFound) this%blockNameFound = this%blockName @@ -202,7 +206,7 @@ subroutine GetNextLine(this, endOfBlock) ! -- read next line loop1: do if (lineread) exit loop1 - call u9rdcom(this%iuext, this%iout, this%line, ierr) + call this%line_reader%rdcom(this%iuext, this%iout, this%line, ierr) this%lloc = 1 call urword(this%line, this%lloc, istart, istop, 0, ival, rval, & this%iout, this%iuext) @@ -587,4 +591,241 @@ subroutine DevOpt(this) return end subroutine DevOpt + ! -- static methods previously in InputOutput + !> @brief Find a block in a file + !! + !! Subroutine to read from a file until the tag (ctag) for a block is + !! is found. Return isfound with true, if found. + !! + !< + subroutine uget_block(line_reader, iin, iout, ctag, ierr, isfound, & + lloc, line, iuext, blockRequired, supportopenclose) + implicit none + ! -- dummy variables + type(LongLineReaderType), intent(inout) :: line_reader + integer(I4B), intent(in) :: iin !< file unit + integer(I4B), intent(in) :: iout !< output listing file unit + character(len=*), intent(in) :: ctag !< block tag + integer(I4B), intent(out) :: ierr !< error + logical, intent(inout) :: isfound !< boolean indicating if the block was found + integer(I4B), intent(inout) :: lloc !< position in line + character(len=:), allocatable, intent(inout) :: line !< line + integer(I4B), intent(inout) :: iuext !< external file unit number + logical, optional, intent(in) :: blockRequired !< boolean indicating if the block is required + logical, optional, intent(in) :: supportopenclose !< boolean indicating if the block supports open/close + ! -- local variables + integer(I4B) :: istart + integer(I4B) :: istop + integer(I4B) :: ival + integer(I4B) :: lloc2 + real(DP) :: rval + character(len=:), allocatable :: line2 + character(len=LINELENGTH) :: fname + character(len=MAXCHARLEN) :: ermsg + logical :: supportoc, blockRequiredLocal + ! + ! -- code + if (present(blockRequired)) then + blockRequiredLocal = blockRequired + else + blockRequiredLocal = .true. + end if + supportoc = .false. + if (present(supportopenclose)) then + supportoc = supportopenclose + end if + iuext = iin + isfound = .false. + mainloop: do + lloc = 1 + call line_reader%rdcom(iin, iout, line, ierr) + if (ierr < 0) then + if (blockRequiredLocal) then + ermsg = 'Required block "'//trim(ctag)// & + '" not found. Found end of file instead.' + call store_error(ermsg) + call store_error_unit(iuext) + end if + ! block not found so exit + exit + end if + call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) + if (line(istart:istop) == 'BEGIN') then + call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) + if (line(istart:istop) == ctag) then + isfound = .true. + if (supportoc) then + ! Look for OPEN/CLOSE on 1st line after line starting with BEGIN + call line_reader%rdcom(iin, iout, line2, ierr) + if (ierr < 0) exit + lloc2 = 1 + call urword(line2, lloc2, istart, istop, 1, ival, rval, iin, iout) + if (line2(istart:istop) == 'OPEN/CLOSE') then + ! -- Get filename and preserve case + call urword(line2, lloc2, istart, istop, 0, ival, rval, iin, iout) + fname = line2(istart:istop) + ! If line contains '(BINARY)' or 'SFAC', handle this block elsewhere + chk: do + call urword(line2, lloc2, istart, istop, 1, ival, rval, iin, iout) + if (line2(istart:istop) == '') exit chk + if (line2(istart:istop) == '(BINARY)' .or. & + line2(istart:istop) == 'SFAC') then + call line_reader%bkspc(iin) + exit mainloop + end if + end do chk + iuext = io_getunit() + call openfile(iuext, iout, fname, 'OPEN/CLOSE') + else + call line_reader%bkspc(iin) + end if + end if + else + if (blockRequiredLocal) then + ermsg = 'Error: Required block "'//trim(ctag)// & + '" not found. Found block "'//line(istart:istop)// & + '" instead.' + call store_error(ermsg) + call store_error_unit(iuext) + else + call line_reader%bkspc(iin) + end if + end if + exit mainloop + else if (line(istart:istop) == 'END') then + call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) + if (line(istart:istop) == ctag) then + ermsg = 'Error: Looking for BEGIN '//trim(ctag)// & + ' but found END '//line(istart:istop)// & + ' instead.' + call store_error(ermsg) + call store_error_unit(iuext) + end if + end if + end do mainloop + ! + ! -- return + return + end subroutine uget_block + + !> @brief Find the next block in a file + !! + !! Subroutine to read from a file until next block is found. + !! Return isfound with true, if found, and return the block name. + !! + !< + subroutine uget_any_block(line_reader, iin, iout, isfound, & + lloc, line, ctagfound, iuext) + implicit none + ! -- dummy variables + type(LongLineReaderType), intent(inout) :: line_reader + integer(I4B), intent(in) :: iin !< file unit number + integer(I4B), intent(in) :: iout !< output listing file unit + logical, intent(inout) :: isfound !< boolean indicating if a block was found + integer(I4B), intent(inout) :: lloc !< position in line + character(len=:), allocatable, intent(inout) :: line !< line + character(len=*), intent(out) :: ctagfound !< block name + integer(I4B), intent(inout) :: iuext !< external file unit number + ! -- local variables + integer(I4B) :: ierr, istart, istop + integer(I4B) :: ival, lloc2 + real(DP) :: rval + character(len=100) :: ermsg + character(len=:), allocatable :: line2 + character(len=LINELENGTH) :: fname + ! + ! -- code + isfound = .false. + ctagfound = '' + iuext = iin + do + lloc = 1 + call line_reader%rdcom(iin, iout, line, ierr) + if (ierr < 0) exit + call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) + if (line(istart:istop) == 'BEGIN') then + call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) + if (line(istart:istop) /= '') then + isfound = .true. + ctagfound = line(istart:istop) + call line_reader%rdcom(iin, iout, line2, ierr) + if (ierr < 0) exit + lloc2 = 1 + call urword(line2, lloc2, istart, istop, 1, ival, rval, iout, iin) + if (line2(istart:istop) == 'OPEN/CLOSE') then + iuext = io_getunit() + call urword(line2, lloc2, istart, istop, 0, ival, rval, iout, iin) + fname = line2(istart:istop) + call openfile(iuext, iout, fname, 'OPEN/CLOSE') + else + call line_reader%bkspc(iin) + end if + else + ermsg = 'Block name missing in file.' + call store_error(ermsg) + call store_error_unit(iin) + end if + exit + end if + end do + return + end subroutine uget_any_block + + !> @brief Evaluate if the end of a block has been found + !! + !! Subroutine to evaluate if the end of a block has been found. Abnormal + !! termination if 'begin' is found or if 'end' encountered with + !! incorrect tag. + !! + !< + subroutine uterminate_block(iin, iout, key, ctag, lloc, line, ierr, iuext) + implicit none + ! -- dummy variables + integer(I4B), intent(in) :: iin !< file unit number + integer(I4B), intent(in) :: iout !< output listing file unit number + character(len=*), intent(in) :: key !< keyword in block + character(len=*), intent(in) :: ctag !< block name + integer(I4B), intent(inout) :: lloc !< position in line + character(len=*), intent(inout) :: line !< line + integer(I4B), intent(inout) :: ierr !< error + integer(I4B), intent(inout) :: iuext !< external file unit number + ! -- local variables + character(len=LENBIGLINE) :: ermsg + integer(I4B) :: istart + integer(I4B) :: istop + integer(I4B) :: ival + real(DP) :: rval + ! -- format +1 format('ERROR. "', A, '" DETECTED WITHOUT "', A, '". ', '"END', 1X, A, & + '" MUST BE USED TO END ', A, '.') +2 format('ERROR. "', A, '" DETECTED BEFORE "END', 1X, A, '". ', '"END', 1X, A, & + '" MUST BE USED TO END ', A, '.') + ! + ! -- code + ierr = 1 + select case (key) + case ('END') + call urword(line, lloc, istart, istop, 1, ival, rval, iout, iin) + if (line(istart:istop) /= ctag) then + write (ermsg, 1) trim(key), trim(ctag), trim(ctag), trim(ctag) + call store_error(ermsg) + call store_error_unit(iin) + else + ierr = 0 + if (iuext /= iin) then + ! -- close external file + close (iuext) + iuext = iin + end if + end if + case ('BEGIN') + write (ermsg, 2) trim(key), trim(ctag), trim(ctag), trim(ctag) + call store_error(ermsg) + call store_error_unit(iin) + end select + ! + ! -- return + return + end subroutine uterminate_block + end module BlockParserModule diff --git a/src/Utilities/InputOutput.f90 b/src/Utilities/InputOutput.f90 index 8d12c1bd21b..ed27ffd90e6 100644 --- a/src/Utilities/InputOutput.f90 +++ b/src/Utilities/InputOutput.f90 @@ -13,14 +13,14 @@ module InputOutputModule DZERO use GenericUtilitiesModule, only: is_same, sim_message private - public :: GetUnit, uget_block, & - uterminate_block, UPCASE, URWORD, ULSTLB, UBDSV4, & + public :: GetUnit, & + UPCASE, URWORD, ULSTLB, UBDSV4, & ubdsv06, UBDSVB, UCOLNO, ULAPRW, & ULASAV, ubdsv1, ubdsvc, ubdsvd, UWWORD, & same_word, get_node, get_ijk, unitinquire, & ParseLine, ulaprufw, openfile, & linear_interpolate, lowcase, & - read_line, uget_any_block, & + read_line, & GetFileFromPath, extract_idnum_or_bndname, urdaux, & get_jk, print_format, BuildFixedFormat, & BuildFloatFormat, BuildIntFormat, fseek_stream, & @@ -201,239 +201,6 @@ function getunit() return end function getunit - !> @brief Find a block in a file - !! - !! Subroutine to read from a file until the tag (ctag) for a block is - !! is found. Return isfound with true, if found. - !! - !< - subroutine uget_block(iin, iout, ctag, ierr, isfound, lloc, line, iuext, & - blockRequired, supportopenclose) - implicit none - ! -- dummy variables - integer(I4B), intent(in) :: iin !< file unit - integer(I4B), intent(in) :: iout !< output listing file unit - character (len=*), intent(in) :: ctag !< block tag - integer(I4B), intent(out) :: ierr !< error - logical, intent(inout) :: isfound !< boolean indicating if the block was found - integer(I4B), intent(inout) :: lloc !< position in line - character (len=:), allocatable, intent(inout) :: line !< line - integer(I4B), intent(inout) :: iuext !< external file unit number - logical, optional, intent(in) :: blockRequired !< boolean indicating if the block is required - logical, optional, intent(in) :: supportopenclose !< boolean indicating if the block supports open/close - ! -- local variables - integer(I4B) :: istart - integer(I4B) :: istop - integer(I4B) :: ival - integer(I4B) :: lloc2 - real(DP) :: rval - character (len=:), allocatable :: line2 - character(len=LINELENGTH) :: fname - character(len=MAXCHARLEN) :: ermsg - logical :: supportoc, blockRequiredLocal - ! - ! -- code - if (present(blockRequired)) then - blockRequiredLocal = blockRequired - else - blockRequiredLocal = .true. - endif - supportoc = .false. - if (present(supportopenclose)) then - supportoc = supportopenclose - endif - iuext = iin - isfound = .false. - mainloop: do - lloc = 1 - call u9rdcom(iin, iout, line, ierr) - if (ierr < 0) then - if (blockRequiredLocal) then - ermsg = 'Required block "' // trim(ctag) // & - '" not found. Found end of file instead.' - call store_error(ermsg) - call store_error_unit(iuext) - end if - ! block not found so exit - exit - end if - call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) - if (line(istart:istop) == 'BEGIN') then - call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) - if (line(istart:istop) == ctag) then - isfound = .true. - if (supportoc) then - ! Look for OPEN/CLOSE on 1st line after line starting with BEGIN - call u9rdcom(iin, iout, line2, ierr) - if (ierr < 0) exit - lloc2 = 1 - call urword(line2, lloc2, istart, istop, 1, ival, rval, iin, iout) - if (line2(istart:istop) == 'OPEN/CLOSE') then - ! -- Get filename and preserve case - call urword(line2, lloc2, istart, istop, 0, ival, rval, iin, iout) - fname = line2(istart:istop) - ! If line contains '(BINARY)' or 'SFAC', handle this block elsewhere - chk: do - call urword(line2, lloc2, istart, istop, 1, ival, rval, iin, iout) - if (line2(istart:istop) == '') exit chk - if (line2(istart:istop) == '(BINARY)' .or. & - line2(istart:istop) == 'SFAC') then - backspace(iin) - exit mainloop - end if - end do chk - iuext = GetUnit() - call openfile(iuext,iout,fname,'OPEN/CLOSE') - else - backspace(iin) - end if - end if - else - if (blockRequiredLocal) then - ermsg = 'Error: Required block "' // trim(ctag) // & - '" not found. Found block "' // line(istart:istop) // & - '" instead.' - call store_error(ermsg) - call store_error_unit(iuext) - else - backspace(iin) - endif - end if - exit mainloop - else if (line(istart:istop) == 'END') then - call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) - if (line(istart:istop) == ctag) then - ermsg = 'Error: Looking for BEGIN ' // trim(ctag) // & - ' but found END ' // line(istart:istop) // & - ' instead.' - call store_error(ermsg) - call store_error_unit(iuext) - endif - end if - end do mainloop - ! - ! -- return - return - end subroutine uget_block - - !> @brief Find the next block in a file - !! - !! Subroutine to read from a file until next block is found. - !! Return isfound with true, if found, and return the block name. - !! - !< - subroutine uget_any_block(iin,iout,isfound,lloc,line,ctagfound,iuext) - implicit none - ! -- dummy variables - integer(I4B), intent(in) :: iin !< file unit number - integer(I4B), intent(in) :: iout !< output listing file unit - logical, intent(inout) :: isfound !< boolean indicating if a block was found - integer(I4B), intent(inout) :: lloc !< position in line - character (len=:), allocatable, intent(inout) :: line !< line - character(len=*), intent(out) :: ctagfound !< block name - integer(I4B), intent(inout) :: iuext !< external file unit number - ! -- local variables - integer(I4B) :: ierr, istart, istop - integer(I4B) :: ival, lloc2 - real(DP) :: rval - character(len=100) :: ermsg - character (len=:), allocatable :: line2 - character(len=LINELENGTH) :: fname - ! - ! -- code - isfound = .false. - ctagfound = '' - iuext = iin - do - lloc = 1 - call u9rdcom(iin,iout,line,ierr) - if (ierr < 0) exit - call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) - if (line(istart:istop) == 'BEGIN') then - call urword(line, lloc, istart, istop, 1, ival, rval, iin, iout) - if (line(istart:istop) /= '') then - isfound = .true. - ctagfound = line(istart:istop) - call u9rdcom(iin,iout,line2,ierr) - if (ierr < 0) exit - lloc2 = 1 - call urword(line2,lloc2,istart,istop,1,ival,rval,iout,iin) - if (line2(istart:istop) == 'OPEN/CLOSE') then - iuext = GetUnit() - call urword(line2,lloc2,istart,istop,0,ival,rval,iout,iin) - fname = line2(istart:istop) - call openfile(iuext,iout,fname,'OPEN/CLOSE') - else - backspace(iin) - endif - else - ermsg = 'Block name missing in file.' - call store_error(ermsg) - call store_error_unit(iin) - end if - exit - end if - end do - return - end subroutine uget_any_block - - !> @brief Evaluate if the end of a block has been found - !! - !! Subroutine to evaluate if the end of a block has been found. Abnormal - !! termination if 'begin' is found or if 'end' encountered with - !! incorrect tag. - !! - !< - subroutine uterminate_block(iin,iout,key,ctag,lloc,line,ierr,iuext) - implicit none - ! -- dummy variables - integer(I4B), intent(in) :: iin !< file unit number - integer(I4B), intent(in) :: iout !< output listing file unit number - character (len=*), intent(in) :: key !< keyword in block - character (len=*), intent(in) :: ctag !< block name - integer(I4B), intent(inout) :: lloc !< position in line - character (len=*), intent(inout) :: line !< line - integer(I4B), intent(inout) :: ierr !< error - integer(I4B), intent(inout) :: iuext !< external file unit number - ! -- local variables - character(len=LENBIGLINE) :: ermsg - integer(I4B) :: istart - integer(I4B) :: istop - integer(I4B) :: ival - real(DP) :: rval - ! -- format -1 format('ERROR. "',A,'" DETECTED WITHOUT "',A,'". ','"END',1X,A, & - '" MUST BE USED TO END ',A,'.') -2 format('ERROR. "',A,'" DETECTED BEFORE "END',1X,A,'". ','"END',1X,A, & - '" MUST BE USED TO END ',A,'.') - ! - ! -- code - ierr = 1 - select case(key) - case ('END') - call urword(line, lloc, istart, istop, 1, ival, rval, iout, iin) - if (line(istart:istop) /= ctag) then - write(ermsg, 1) trim(key), trim(ctag), trim(ctag), trim(ctag) - call store_error(ermsg) - call store_error_unit(iin) - else - ierr = 0 - if (iuext /= iin) then - ! -- close external file - close(iuext) - iuext = iin - endif - end if - case ('BEGIN') - write(ermsg, 2) trim(key), trim(ctag), trim(ctag), trim(ctag) - call store_error(ermsg) - call store_error_unit(iin) - end select - ! - ! -- return - return - end subroutine uterminate_block - !> @brief Convert to upper case !! !! Subroutine to convert a character string to upper case. @@ -2128,9 +1895,8 @@ subroutine u9rdcom(iin, iout, line, ierr) pcomments: do call get_line(iin, line, ierr) if (ierr == IOSTAT_END) then - ! -- End of file reached. - ! -- Backspace is needed for gfortran. - backspace(iin) + ! -- End of file reached. Return with ierr = IOSTAT_END + ! and line as an empty string line = ' ' exit pcomments elseif (ierr /= 0) then diff --git a/src/Utilities/ListReader.f90 b/src/Utilities/ListReader.f90 index fc5150937ac..959b1cdd0e2 100644 --- a/src/Utilities/ListReader.f90 +++ b/src/Utilities/ListReader.f90 @@ -6,6 +6,7 @@ module ListReaderModule LENAUXNAME, LENLISTLABEL, DONE use SimVariablesModule, only: errmsg use SimModule, only: store_error, count_errors, store_error_unit + use LongLineReaderModule, only: LongLineReaderType implicit none private @@ -41,6 +42,7 @@ module ListReaderModule integer(I4B), dimension(:), allocatable :: idxtxtauxcol ! col locations of text in auxvar character(len=LENTIMESERIESNAME), dimension(:), allocatable :: txtrlist ! text found in rlist character(len=LENTIMESERIESNAME), dimension(:), allocatable :: txtauxvar ! text found in auxvar + type(LongLineReaderType), pointer :: line_reader => null() contains procedure :: read_list procedure :: write_list @@ -53,8 +55,9 @@ module ListReaderModule contains - subroutine read_list(this, in, iout, nlist, inamedbound, mshape, nodelist, & - rlist, auxvar, auxname, boundname, label) + subroutine read_list(this, line_reader, in, iout, nlist, inamedbound, & + mshape, nodelist, rlist, auxvar, auxname, boundname, & + label) ! ****************************************************************************** ! init -- Initialize the reader ! ****************************************************************************** @@ -65,6 +68,7 @@ subroutine read_list(this, in, iout, nlist, inamedbound, mshape, nodelist, & use ConstantsModule, only: LENBOUNDNAME ! -- dummy class(ListReaderType) :: this + type(LongLineReaderType), intent(inout), target :: line_reader integer(I4B), intent(in) :: in integer(I4B), intent(in) :: iout integer(I4B), intent(inout) :: nlist @@ -95,6 +99,7 @@ subroutine read_list(this, in, iout, nlist, inamedbound, mshape, nodelist, & this%auxvar => auxvar this%auxname => auxname this%boundname => boundname + this%line_reader => line_reader ! ! -- Allocate arrays for storing text and text locations if (.not. allocated(this%idxtxtrow)) allocate (this%idxtxtrow(0)) @@ -125,7 +130,7 @@ subroutine read_control_record(this) ! SPECIFICATIONS: ! ------------------------------------------------------------------------------ ! -- modules - use InputOutputModule, only: u9rdcom, urword + use InputOutputModule, only: urword ! -- dummy class(ListReaderType) :: this ! -- local @@ -142,7 +147,7 @@ subroutine read_control_record(this) this%ibinary = 0 ! ! -- Read to the first non-commented line - call u9rdcom(this%in, this%iout, this%line, this%ierr) + call this%line_reader%rdcom(this%in, this%iout, this%line, this%ierr) this%lloc = 1 call urword(this%line, this%lloc, this%istart, this%istop, 1, idum, r, & this%iout, this%in) @@ -167,7 +172,7 @@ subroutine set_openclose(this) ! SPECIFICATIONS: ! ------------------------------------------------------------------------------ ! -- modules - use InputOutputModule, only: u9rdcom, urword, openfile + use InputOutputModule, only: urword, openfile use OpenSpecModule, only: form, access use ConstantsModule, only: LINELENGTH ! -- dummy @@ -237,8 +242,9 @@ subroutine set_openclose(this) ! ! -- Read the first line from inlist to be consistent with how the list is ! read when it is included in the package input file - if (this%ibinary /= 1) call u9rdcom(this%inlist, this%iout, this%line, & - this%ierr) + if (this%ibinary /= 1) & + call this%line_reader%rdcom(this%inlist, this%iout, this%line, & + this%ierr) ! ! -- return return @@ -394,7 +400,7 @@ subroutine read_ascii(this) ! ------------------------------------------------------------------------------ ! -- modules use ConstantsModule, only: LENBOUNDNAME, LINELENGTH, DZERO - use InputOutputModule, only: u9rdcom, urword, get_node + use InputOutputModule, only: urword, get_node use ArrayHandlersModule, only: ExpandArray use TdisModule, only: kper ! -- dummy @@ -427,7 +433,8 @@ subroutine read_ascii(this) readloop: do ! ! -- First line was already read, so don't read again - if (ii /= 1) call u9rdcom(this%inlist, 0, this%line, this%ierr) + if (ii /= 1) & + call this%line_reader%rdcom(this%inlist, 0, this%line, this%ierr) ! ! -- If this is an unknown-length list, then check for END. ! If found, then backspace, set nlist, and exit readloop. @@ -436,10 +443,10 @@ subroutine read_ascii(this) call urword(this%line, this%lloc, this%istart, this%istop, 1, idum, r, & this%iout, this%inlist) if (this%line(this%istart:this%istop) == 'END' .or. this%ierr < 0) then - ! If ierr < 0, backspace was already performed in u9rdcom, so only - ! need to backspace if END was found. + ! If END was found then call line_reader backspace + ! emulator so that caller can proceed with reading END. if (this%ierr == 0) then - backspace (this%inlist) + call this%line_reader%bkspc(this%inlist) end if this%nlist = ii - 1 exit readloop diff --git a/src/Utilities/LongLineReader.f90 b/src/Utilities/LongLineReader.f90 new file mode 100644 index 00000000000..bce57afd45d --- /dev/null +++ b/src/Utilities/LongLineReader.f90 @@ -0,0 +1,117 @@ +!> @brief This module contains the LongLineReaderType +!! +!! The LongLineReader is a utility for reading text lines +!! from mf6 input files. It calls u9rdcom (which calls +!! get_line) to read the first non-commented line of an +!! input file. The LongLineReader can emulate the Fortran +!! backspace command by calling the bkspc method, which stores +!! the current line in last_line, and will return last_line +!! upon the next call to rdcom. The LongLineReader was +!! implemented to replace all Fortran backspace calls, due +!! to a bug in ifort and ifx that prevented the backspace +!! command from working properly with non-advancing IO. +!! +!< +module LongLineReaderModule + + use, intrinsic :: iso_fortran_env, only: IOSTAT_END + use KindModule, only: I4B + use SimModule, only: store_error + use InputOutputModule, only: u9rdcom + + implicit none + + private + public :: LongLineReaderType + + !> @brief LongLineReaderType + !! + !! Object for reading input from mf6 input files + !! + !< + type :: LongLineReaderType + + character(len=:), allocatable :: line + character(len=:), allocatable :: last_line + integer(I4B) :: nbackspace = 0 + integer(I4B) :: iostat = 0 + integer(I4B) :: last_unit = 0 + + contains + + procedure :: bkspc + procedure :: rdcom + + end type LongLineReaderType + +contains + + !> @brief Return the first non-comment line + !! + !! Skip through any comments and return the first + !! non-commented line. If an end of file was + !! encountered previously, then return a blank line. + !! If a backspace was called prior to this call, + !! then do not read a new line and return last_line + !! instead. + !! + !< + subroutine rdcom(this, iu, iout, line, ierr) + class(LongLineReaderType) :: this + integer(I4B), intent(in) :: iu + integer(I4B), intent(in) :: iout + character(len=:), intent(inout), allocatable :: line + integer(I4B), intent(inout) :: ierr + + ierr = 0 + + ! If using this reader to read from a new file + ! then reset state + if (iu /= this%last_unit) then + this%nbackspace = 0 + this%iostat = 0 + end if + + if (this%nbackspace == 1) then + ! If backspace was called, then return last line + if (allocated(line)) deallocate (line) + allocate (character(len=len(this%last_line) + 1) :: line) + line(:) = this%last_line(:) + this%nbackspace = 0 + else + ! if end of file was reached previously, then return a + ! blank line and return ierr as IOSTAT_END + if (this%iostat == IOSTAT_END) then + line = ' ' + ierr = IOSTAT_END + else + call u9rdcom(iu, iout, line, ierr) + end if + this%last_line = line + this%iostat = ierr + end if + this%last_unit = iu + return + end subroutine rdcom + + !> @brief Emulate a Fortan backspace + !! + !! Emulate a fortran backspace call by storing + !! the current line in long_line + !! + !< + subroutine bkspc(this, iin) + class(LongLineReaderType) :: this + integer(I4B), intent(in) :: iin + if (this%nbackspace > 0) then + call store_error( & + "Programming error in LongLineReaderType%bkspc(). Backspace & + & called more than once for an open file.", & + terminate=.true.) + else + this%nbackspace = 1 + end if + return + end subroutine bkspc + +end module LongLineReaderModule diff --git a/src/meson.build b/src/meson.build index 5f57e065d7b..eea4e2a0618 100644 --- a/src/meson.build +++ b/src/meson.build @@ -209,6 +209,7 @@ modflow_sources = files( 'Utilities' / 'kind.f90', 'Utilities' / 'List.f90', 'Utilities' / 'ListReader.f90', + 'Utilities' / 'LongLineReader.f90', 'Utilities' / 'Message.f90', 'Utilities' / 'OpenSpec.f90', 'Utilities' / 'PackageBudget.f90', diff --git a/utils/mf5to6/make/makedefaults b/utils/mf5to6/make/makedefaults index fdc2fdb51cb..bcea1425b67 100644 --- a/utils/mf5to6/make/makedefaults +++ b/utils/mf5to6/make/makedefaults @@ -1,4 +1,4 @@ -# makedefaults created by pymake (version 1.2.7) for the 'mf5to6' executable. +# makedefaults created by pymake (version 1.2.9.dev0) for the 'mf5to6' executable. # determine OS ifeq ($(OS), Windows_NT) diff --git a/utils/mf5to6/make/makefile b/utils/mf5to6/make/makefile index d8eb7e6963d..c455f6d76b0 100644 --- a/utils/mf5to6/make/makefile +++ b/utils/mf5to6/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.7) for the 'mf5to6' executable. +# makefile created by pymake (version 1.2.9.dev0) for the 'mf5to6' executable. include ./makedefaults @@ -46,6 +46,7 @@ $(OBJDIR)/CharString.o \ $(OBJDIR)/Memory.o \ $(OBJDIR)/List.o \ $(OBJDIR)/MemoryList.o \ +$(OBJDIR)/LongLineReader.o \ $(OBJDIR)/Utilities.o \ $(OBJDIR)/ConstantsPHMF.o \ $(OBJDIR)/MemoryManager.o \ diff --git a/utils/mf5to6/msvs/mf5to6.vfproj b/utils/mf5to6/msvs/mf5to6.vfproj index 0fdd1b971ce..a027537f04b 100644 --- a/utils/mf5to6/msvs/mf5to6.vfproj +++ b/utils/mf5to6/msvs/mf5to6.vfproj @@ -1,28 +1,32 @@ - + + - - - - - - - - - + + + + + + + + + + - - - - - - - - - + + + + + + + + + + + @@ -93,12 +97,16 @@ - + + + + + @@ -190,5 +198,7 @@ - - + + + + diff --git a/utils/mf5to6/pymake/extrafiles.txt b/utils/mf5to6/pymake/extrafiles.txt index 1699e51ed46..89a298d0b93 100644 --- a/utils/mf5to6/pymake/extrafiles.txt +++ b/utils/mf5to6/pymake/extrafiles.txt @@ -13,6 +13,7 @@ ../../../src/Utilities/InputOutput.f90 ../../../src/Utilities/kind.f90 ../../../src/Utilities/List.f90 +../../../src/Utilities/LongLineReader.f90 ../../../src/Utilities/OpenSpec.f90 ../../../src/Utilities/defmacro.F90 ../../../src/Utilities/version.f90 diff --git a/utils/mf5to6/src/Connection.f90 b/utils/mf5to6/src/Connection.f90 index 561d737a049..8c0d8fec34b 100644 --- a/utils/mf5to6/src/Connection.f90 +++ b/utils/mf5to6/src/Connection.f90 @@ -67,13 +67,43 @@ subroutine WriteGhostNodeCorrection(this, iu, numalphaj) ! select case (numalphaj) case (1) + if (this%alphaj1 == 0) then + this%k1 = 0 + this%i1 = 0 + this%j1 = 0 + end if write(iu,10)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, & this%k1, this%i1, this%j1, this%alphaj1 case (2) + if (this%alphaj1 == 0.) then + this%k1 = 0 + this%i1 = 0 + this%j1 = 0 + end if + if (this%alphaj2 == 0.) then + this%k2 = 0 + this%i2 = 0 + this%j2 = 0 + end if write(iu,20)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, & this%k1, this%i1, this%j1, this%k2, this%i2, this%j2, & this%alphaj1, this%alphaj2 case (3) + if (this%alphaj1 == 0.) then + this%k1 = 0 + this%i1 = 0 + this%j1 = 0 + end if + if (this%alphaj2 == 0.) then + this%k2 = 0 + this%i2 = 0 + this%j2 = 0 + end if + if (this%alphaj12 == 0.) then + this%k12 = 0 + this%i12 = 0 + this%j12 = 0 + end if write(iu,30)this%kp, this%ip, this%jp, this%kc, this%ic, this%jc, & this%k1, this%i1, this%j1, this%k2, this%i2, this%j2, & this%k12, this%i12, this%j12, this%alphaj1, this%alphaj2, & diff --git a/utils/mf5to6/src/Preproc/ObsBlock.f90 b/utils/mf5to6/src/Preproc/ObsBlock.f90 index 9a16e3bfd50..61f3e5bad78 100644 --- a/utils/mf5to6/src/Preproc/ObsBlock.f90 +++ b/utils/mf5to6/src/Preproc/ObsBlock.f90 @@ -7,7 +7,7 @@ module ObsBlockModule use ConstantsPHMFModule, only: CONTINUOUS, SINGLE, LENOBSNAMENEW use DnmDis3dModule, only: Dis3dType use GlobalVariablesPHMFModule, only: verbose - use InputOutputModule, only: UPCASE, URWORD, uterminate_block + use InputOutputModule, only: UPCASE, URWORD use ListModule, only: ListType use ObserveModule, only: ObserveType, AddObserveToList, & GetObserveFromList, ConstructObservation diff --git a/utils/mf5to6/src/Preproc/Preproc.f90 b/utils/mf5to6/src/Preproc/Preproc.f90 index c320028b428..df58ae136e1 100644 --- a/utils/mf5to6/src/Preproc/Preproc.f90 +++ b/utils/mf5to6/src/Preproc/Preproc.f90 @@ -14,9 +14,7 @@ module PreprocModule use GLOBAL, only: NCOL, NROW, DELC, DELR use globalPHMF, only: ioutPHMF, outfile use GlobalVariablesPHMFModule, only: prognamPHMF, verbose, vnam - use InputOutputModule, only: GetUnit, uget_block, urword, & - uterminate_block, GetUnit, openfile, & - uget_any_block + use InputOutputModule, only: GetUnit, urword, GetUnit, openfile use ListModule, only: ListType use ObsBlockModule, only: ObsBlockType, ConstructObsBlockType, & AddObsBlockToList, GetObsBlockFromList @@ -203,8 +201,6 @@ subroutine read_options(this) ierr = 0 ! ! -- get BEGIN line of OPTIONS block -! call uget_block(iin, 0, blockTypeWanted, ierr, found, & -! lloc, line, iuext, continueread) call this%parser%GetBlock('OPTIONS', found, ierr, supportOpenClose=.true.) if (ierr /= 0) then ! end of file @@ -733,7 +729,6 @@ subroutine read_any_block(this, iu ,k, eof, dis3d, WriteBeginEnd) ! ! -- Read any block as long as it's SINGLE or CONTINUOUS. lloc = 1 -! call uget_any_block(iu, this%iout, isfound, lloc, line, ctagfound, iuext) call this%parser%GetBlock('*', isfound, ierr, .true., & .false., ctagfound) if (.not. isfound) then diff --git a/utils/mf5to6/src/Preproc/Utilities.f90 b/utils/mf5to6/src/Preproc/Utilities.f90 index 3d97b4c0618..13e7ec9c3f4 100644 --- a/utils/mf5to6/src/Preproc/Utilities.f90 +++ b/utils/mf5to6/src/Preproc/Utilities.f90 @@ -3,8 +3,7 @@ module UtilitiesModule use ConstantsModule, only: MAXCHARLEN, DZERO, MAXCHARLEN use GlobalVariablesModule, only: optfile, PathToPostObsMf, ScriptType, & verbose, echo - use InputOutputModule, only: GetUnit, openfile, UPCASE, URWORD, & - uget_block, uterminate_block, u9rdcom + use InputOutputModule, only: GetUnit, openfile, UPCASE, URWORD use SimModule, onlY: store_error, store_note, store_warning, ustop private @@ -15,7 +14,7 @@ module UtilitiesModule BuildArrayFormat, Write1dValues, & Write2dValues, Write3dValues, findcell, & close_file, GreaterOf, GreatestOf, RemoveElement, & - get_extension, ReadMf5to6Options, count_file_records, & + get_extension, count_file_records, & CalcContribFactors, PhmfOption interface RemoveElement @@ -876,65 +875,6 @@ subroutine get_extension(name, ext) return end subroutine get_extension - subroutine ReadMf5to6Options() - implicit none - ! local - integer :: ierr, istart, istop, iu, idum, icol - double precision :: rdum - character(len=MAXCHARLEN) :: ermsg - character(len=:), allocatable :: line - character(len=10) :: stype - logical :: continueread=.true., found - integer :: iuext - ! - if (optfile /= '') then - iu = GetUnit() - call openfile(iu, 0, optfile, 'OPTIONS', filstat_opt='OLD') - call uget_block(iu, 0, 'OPTIONS', ierr, found, icol, line, iuext, & - continueread) - if (found) then - do - icol = 1 - call u9rdcom(iu, 0, line, ierr) - call urword(line, icol, istart, istop, 1, idum, rdum, 0, iu) - select case (line(istart:istop)) - case ('PATHTOPOSTOBSMF') - call urword(line, icol, istart, istop, 0, idum, rdum, 0, iu) - PathToPostObsMf = line(istart:istop) - case ('SCRIPT') - call urword(line, icol, istart, istop, 1, idum, rdum, 0, iu) - stype = line(istart:istop) - select case (stype) - case ('BATCH') - ScriptType = 'BATCH' - case ('PYTHON') - ScriptType = 'PYTHON' - case default - ermsg = 'Unknown Script option: ' // line(istart:istop) - call store_error(ermsg) - call ustop() - end select - case default - ermsg = 'Unknown Mf5to6 option: ' // line(istart:istop) - call store_error(ermsg) - call ustop() - case ('END','BEGIN') - call uterminate_block(iu,0,line(istart:istop), & - 'OPTIONS', icol,line,ierr, iuext) - if(ierr==0) exit - end select - enddo - close(iu) - else - ermsg = 'Mf5to6 options file not found: ' // trim(optfile) - call store_error(ermsg) - call ustop() - endif - endif - ! - return - end subroutine ReadMf5to6Options - function count_file_records(filename) result(nrecs) ! Open a text file, count the number of records in it, and close the file. ! dummy diff --git a/utils/mf5to6/src/mf5to6.f90 b/utils/mf5to6/src/mf5to6.f90 index 1e7b07dc8da..ecf2131c8de 100644 --- a/utils/mf5to6/src/mf5to6.f90 +++ b/utils/mf5to6/src/mf5to6.f90 @@ -19,7 +19,7 @@ program mf5to6 use SimFileWriterModule, only: SimFileWriterType use SimModule, only: ustop use SimListVariablesModule, only: SimMovers - use UtilitiesModule, only: GetArgs, ReadMf5to6Options, PhmfOption + use UtilitiesModule, only: GetArgs, PhmfOption ! implicit none integer :: iexg, igrid, ispw, iu @@ -69,7 +69,6 @@ program mf5to6 ! provide a command-prompt instruction that will run PostObsMF, or maybe ! generate a batch or python file (could also be an option) that would ! run PostObsMF (twice, if there are multilayer head observations). - if (SupportPreproc) call ReadMf5to6Options() SimFileWriter%BaseName = basnam if (ilgr > 0) then ! LGR is active; read and initialize parent and all children. diff --git a/utils/zonebudget/make/makedefaults b/utils/zonebudget/make/makedefaults index 919f5ed9b06..e27dfca39de 100644 --- a/utils/zonebudget/make/makedefaults +++ b/utils/zonebudget/make/makedefaults @@ -1,4 +1,4 @@ -# makedefaults created by pymake (version 1.2.7) for the 'zbud6' executable. +# makedefaults created by pymake (version 1.2.9.dev0) for the 'zbud6' executable. # determine OS ifeq ($(OS), Windows_NT) diff --git a/utils/zonebudget/make/makefile b/utils/zonebudget/make/makefile index 092dc3b0294..1918bd8fee2 100644 --- a/utils/zonebudget/make/makefile +++ b/utils/zonebudget/make/makefile @@ -1,4 +1,4 @@ -# makefile created by pymake (version 1.2.7) for the 'zbud6' executable. +# makefile created by pymake (version 1.2.9.dev0) for the 'zbud6' executable. include ./makedefaults @@ -26,6 +26,7 @@ $(OBJDIR)/Message.o \ $(OBJDIR)/Sim.o \ $(OBJDIR)/OpenSpec.o \ $(OBJDIR)/InputOutput.o \ +$(OBJDIR)/LongLineReader.o \ $(OBJDIR)/sort.o \ $(OBJDIR)/BlockParser.o \ $(OBJDIR)/ArrayReaders.o \ diff --git a/utils/zonebudget/msvs/zonebudget.vfproj b/utils/zonebudget/msvs/zonebudget.vfproj index 682849e3f29..919ef4c8552 100644 --- a/utils/zonebudget/msvs/zonebudget.vfproj +++ b/utils/zonebudget/msvs/zonebudget.vfproj @@ -1,28 +1,32 @@ - + + - - - - - - - - - + + + + + + + + + + - - - - - - - - - + + + + + + + + + + + @@ -34,14 +38,15 @@ - + - + + @@ -52,5 +57,7 @@ - - + + + + diff --git a/utils/zonebudget/pymake/extrafiles.txt b/utils/zonebudget/pymake/extrafiles.txt index c4cfdf23967..0848ab0c086 100644 --- a/utils/zonebudget/pymake/extrafiles.txt +++ b/utils/zonebudget/pymake/extrafiles.txt @@ -7,6 +7,7 @@ ../../../src/Utilities/genericutils.f90 ../../../src/Utilities/InputOutput.f90 ../../../src/Utilities/kind.f90 +../../../src/Utilities/LongLineReader.f90 ../../../src/Utilities/OpenSpec.f90 ../../../src/Utilities/sort.f90 ../../../src/Utilities/Message.f90 From f84accdb11199685e186645200450b23d005c70f Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Wed, 26 Jul 2023 15:17:28 -0500 Subject: [PATCH 03/13] fix(msvs): msvs project file needs to be updated for backspace change (#1312) --- msvs/mf6core.vfproj | 2 +- msvs/mf6lib.vfproj | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/msvs/mf6core.vfproj b/msvs/mf6core.vfproj index b2fe810728e..daf0fa21fe2 100644 --- a/msvs/mf6core.vfproj +++ b/msvs/mf6core.vfproj @@ -365,9 +365,9 @@ - + diff --git a/msvs/mf6lib.vfproj b/msvs/mf6lib.vfproj index e89b383dde7..4028ac77f17 100644 --- a/msvs/mf6lib.vfproj +++ b/msvs/mf6lib.vfproj @@ -150,6 +150,7 @@ + From 6598100a55c1c7694bc2c7cf4ca0761c6da7e166 Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Fri, 28 Jul 2023 09:29:38 -0500 Subject: [PATCH 04/13] feat(pak_cc): add support for parallel package convergence (#1309) * remove `dpak` check in `sln_nur_has_converged()` --- autotest/test_gwf_uzf_gwet.py | 50 ++++++++-------- autotest/test_par_gwf_pakcc.py | 78 +++++++++++++++++++++++++ src/Solution/NumericalSolution.f90 | 94 +++++++++++++++++------------- src/Solution/ParallelSolution.f90 | 30 ++++++++-- 4 files changed, 182 insertions(+), 70 deletions(-) create mode 100644 autotest/test_par_gwf_pakcc.py diff --git a/autotest/test_gwf_uzf_gwet.py b/autotest/test_gwf_uzf_gwet.py index ab7d74556d8..a06115897d3 100644 --- a/autotest/test_gwf_uzf_gwet.py +++ b/autotest/test_gwf_uzf_gwet.py @@ -3,16 +3,20 @@ import flopy import numpy as np import pytest + from framework import TestFramework from simulation import TestSimulation ex = ["uzf_3lay"] +name = "model" iuz_cell_dict = {} cell_iuz_dict = {} +nouter, ninner = 100, 300 +hclose, rclose, relax = 1e-9, 1e-3, 0.97 -def build_model(idx, dir): +def build_model(idx, dir): nlay, nrow, ncol = 3, 1, 10 nper = 5 perlen = [20.0, 20.0, 20.0, 500.0, 2000.0] @@ -23,15 +27,10 @@ def build_model(idx, dir): delc = 1.0 strt = -25 - nouter, ninner = 100, 300 - hclose, rclose, relax = 1e-9, 1e-3, 0.97 - tdis_rc = [] for i in range(nper): tdis_rc.append((perlen[i], nstp[i], tsmult[i])) - name = ex[idx] - # build MODFLOW 6 files ws = dir sim = flopy.mf6.MFSimulation( @@ -263,11 +262,7 @@ def eval_model(sim): ws = sim.simpath sim = flopy.mf6.MFSimulation.load(sim_ws=ws) - fpth = os.path.join(ws, "uzf_3lay.hds") - hobj = flopy.utils.HeadFile(fpth, precision="double") - hds = hobj.get_alldata() - - bpth = os.path.join(ws, "uzf_3lay.cbc") + bpth = ws / f"{name}.cbc" bobj = flopy.utils.CellBudgetFile(bpth, precision="double") bobj.get_unique_record_names() # ' STO-SS' @@ -281,7 +276,7 @@ def eval_model(sim): gwet = bobj.get_data(text="UZF-GWET") gwet = np.array(gwet) - uzpth = os.path.join(ws, "uzf_3lay.uzf.bud") + uzpth = os.path.join(ws, f"{name}.uzf.bud") uzobj = flopy.utils.CellBudgetFile(uzpth, precision="double") uzobj.get_unique_record_names() # ' FLOW-JA-FACE' @@ -303,18 +298,18 @@ def eval_model(sim): gwet_arr = np.zeros( ( tot_stp, - sim.uzf_3lay.dis.nlay.get_data(), - sim.uzf_3lay.dis.nrow.get_data(), - sim.uzf_3lay.dis.ncol.get_data(), + sim.model.dis.nlay.get_data(), + sim.model.dis.nrow.get_data(), + sim.model.dis.ncol.get_data(), ) ) uzet_arr = np.zeros( ( tot_stp, - sim.uzf_3lay.dis.nlay.get_data(), - sim.uzf_3lay.dis.nrow.get_data(), - sim.uzf_3lay.dis.ncol.get_data(), + sim.model.dis.nlay.get_data(), + sim.model.dis.nrow.get_data(), + sim.model.dis.ncol.get_data(), ) ) @@ -336,7 +331,7 @@ def eval_model(sim): uzet_arr[tm, lay, row, col] = itm[2] - uzf_strsPerDat = sim.uzf_3lay.uzf.perioddata.get_data() + uzf_strsPerDat = sim.model.uzf.perioddata.get_data() pet = 0 for tm in range(tot_stp): nstps = 0 @@ -345,8 +340,8 @@ def eval_model(sim): if tm < nstps: break - for i in range(sim.uzf_3lay.dis.nrow.get_data()): - for j in range(sim.uzf_3lay.dis.ncol.get_data()): + for i in range(sim.model.dis.nrow.get_data()): + for j in range(sim.model.dis.ncol.get_data()): if (0, i, j) in cell_iuz_dict: iuz = cell_iuz_dict[ (0, i, j) @@ -370,14 +365,17 @@ def eval_model(sim): print("Finished running checks") -@pytest.mark.parametrize("name", ex) -def test_mf6model(name, function_tmpdir, targets): - ws = str(function_tmpdir) +@pytest.mark.parametrize("idx,name", enumerate(ex)) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = function_tmpdir test = TestFramework() - test.build(build_model, 0, ws) + test.build(build_model, idx, ws) test.run( TestSimulation( - name=name, exe_dict=targets, exfunc=eval_model, idxsim=0 + name=name, + exe_dict=targets, + exfunc=eval_model, + idxsim=idx, ), ws, ) diff --git a/autotest/test_par_gwf_pakcc.py b/autotest/test_par_gwf_pakcc.py new file mode 100644 index 00000000000..b25fcfc37e5 --- /dev/null +++ b/autotest/test_par_gwf_pakcc.py @@ -0,0 +1,78 @@ +import os +import pathlib as pl +from decimal import Decimal + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +# This tests reuses the simulation data in test_gwf_uzf_gwet +# and runs it in parallel on one and two cpus with +# +# so we can compare the parallel coupling of two models +# with a serial model. +# +# This test also checks that Newton under_relaxation works in parallel. + +ex = ["par_uzf_3lay_1p", "par_uzf_3lay_2p"] + + +def build_petsc_db(idx, exdir): + from test_gwf_uzf_gwet import hclose, ninner + + petsc_db_file = os.path.join(exdir, ".petscrc") + with open(petsc_db_file, "w") as petsc_file: + petsc_file.write("-ksp_type bicg\n") + petsc_file.write("-pc_type bjacobi\n") + petsc_file.write("-sub_pc_type ilu\n") + petsc_file.write("-sub_pc_factor_levels 2\n") + petsc_file.write(f"-dvclose {Decimal(hclose):.2E}\n") + petsc_file.write(f"-nitermax {ninner}\n") + petsc_file.write("-options_left no\n") + + +def build_model(idx, exdir): + from test_gwf_uzf_gwet import build_model as build_model_ext + + build_petsc_db(idx, exdir) + sim, dummy = build_model_ext(idx, exdir) + if idx == 1: + sim.set_sim_path(exdir / "working") + sim.write_simulation(silent=True) + mfsplit = flopy.mf6.utils.Mf6Splitter(sim) + split_array = np.zeros((10), dtype=int) + split_array[5:] = 1 + new_sim = mfsplit.split_model(split_array) + new_sim.set_sim_path(exdir) + mfsplit.save_node_mapping(pl.Path(f"{exdir}/mapping.json")) + return new_sim, None + else: + return sim, dummy + + +@pytest.mark.parallel +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework() + test.build(build_model, idx, function_tmpdir) + if idx == 1: + ncpus = 2 + else: + ncpus = 1 + test.run( + TestSimulation( + name=name, + exe_dict=targets, + idxsim=idx, + make_comparison=False, + parallel=True, + ncpus=ncpus, + ), + str(function_tmpdir), + ) diff --git a/src/Solution/NumericalSolution.f90 b/src/Solution/NumericalSolution.f90 index 9bd7c91b94c..c27ca4df333 100644 --- a/src/Solution/NumericalSolution.f90 +++ b/src/Solution/NumericalSolution.f90 @@ -27,7 +27,7 @@ module NumericalSolutionModule AddNumericalExchangeToList, & GetNumericalExchangeFromList use SparseModule, only: sparsematrix - use SimVariablesModule, only: iout, isim_mode + use SimVariablesModule, only: iout, isim_mode, errmsg use SimStagesModule use BlockParserModule, only: BlockParserType use IMSLinearModule @@ -160,6 +160,7 @@ module NumericalSolutionModule ! 'protected' (this can be overridden) procedure :: sln_has_converged + procedure :: sln_package_convergence procedure :: sln_sync_newtonur_flag procedure :: sln_nur_has_converged procedure :: sln_calc_ptc @@ -531,7 +532,6 @@ subroutine sln_ar(this) ! -- local variables class(NumericalModelType), pointer :: mp => null() class(NumericalExchangeType), pointer :: cp => null() - character(len=linelength) :: errmsg character(len=linelength) :: warnmsg character(len=linelength) :: keyword character(len=linelength) :: fname @@ -1502,7 +1502,6 @@ subroutine solve(this, kiter) class(NumericalExchangeType), pointer :: cp => null() character(len=LINELENGTH) :: title character(len=LINELENGTH) :: tag - character(len=LINELENGTH) :: line character(len=LENPAKLOC) :: cmod character(len=LENPAKLOC) :: cpak character(len=LENPAKLOC) :: cpakout @@ -1692,40 +1691,34 @@ subroutine solve(this, kiter) end if end do ! - ! -- evaluate package convergence - if (abs(dpak) > this%dvclose) then - this%icnvg = 0 - ! -- write message to stdout - if (iend /= 0) then - write (line, '(3a)') & - 'PACKAGE (', trim(cpakout), ') CAUSED CONVERGENCE FAILURE' - call sim_message(line) - end if - end if - ! - ! -- write maximum change in package convergence check - if (this%iprims > 0) then - cval = 'Package' - if (this%icnvg /= 1) then - cmsg = ' ' - else - cmsg = '*' - end if - if (len_trim(cpakout) > 0) then - ! - ! -- add data to outertab - call this%outertab%add_term(cval) - call this%outertab%add_term(kiter) - call this%outertab%add_term(' ') - if (this%numtrack > 0) then - call this%outertab%add_term(' ') - call this%outertab%add_term(' ') - call this%outertab%add_term(' ') + ! -- evaluate package convergence - only done if convergence is achieved + if (this%icnvg == 1) then + this%icnvg = this%sln_package_convergence(dpak, cpakout, iend) + ! + ! -- write maximum change in package convergence check + if (this%iprims > 0) then + cval = 'Package' + if (this%icnvg /= 1) then + cmsg = ' ' + else + cmsg = '*' + end if + if (len_trim(cpakout) > 0) then + ! + ! -- add data to outertab + call this%outertab%add_term(cval) + call this%outertab%add_term(kiter) call this%outertab%add_term(' ') + if (this%numtrack > 0) then + call this%outertab%add_term(' ') + call this%outertab%add_term(' ') + call this%outertab%add_term(' ') + call this%outertab%add_term(' ') + end if + call this%outertab%add_term(dpak) + call this%outertab%add_term(cmsg) + call this%outertab%add_term(cpakout) end if - call this%outertab%add_term(dpak) - call this%outertab%add_term(cmsg) - call this%outertab%add_term(cpakout) end if end if ! @@ -1762,7 +1755,7 @@ subroutine solve(this, kiter) call this%sln_maxval(this%neq, this%dxold, dxold_max) ! ! -- evaluate convergence - if (this%sln_nur_has_converged(dxold_max, this%hncg(kiter), dpak)) then + if (this%sln_nur_has_converged(dxold_max, this%hncg(kiter))) then ! ! -- converged this%icnvg = 1 @@ -3144,6 +3137,29 @@ function sln_has_converged(this, max_dvc) result(has_converged) end function sln_has_converged + !> @brief Check package convergence + !< + function sln_package_convergence(this, dpak, cpakout, iend) result(ivalue) + ! dummy + class(NumericalSolutionType) :: this !< NumericalSolutionType instance + real(DP), intent(in) :: dpak !< Newton Under-relaxation flag + character(len=LENPAKLOC), intent(in) :: cpakout + integer(I4B), intent(in) :: iend + ! local + integer(I4B) :: ivalue !< + ivalue = 1 + if (abs(dpak) > this%dvclose) then + ivalue = 0 + ! -- write message to stdout + if (iend /= 0) then + write (errmsg, '(3a)') & + 'PACKAGE (', trim(cpakout), ') CAUSED CONVERGENCE FAILURE' + call sim_message(errmsg) + end if + end if + + end function sln_package_convergence + !> @brief Syncronize Newton Under-relaxation flag !< function sln_sync_newtonur_flag(this, inewtonur) result(ivalue) @@ -3159,18 +3175,16 @@ end function sln_sync_newtonur_flag !> @brief Custom convergence check for when Newton UR has been applied !< - function sln_nur_has_converged(this, dxold_max, hncg, dpak) & + function sln_nur_has_converged(this, dxold_max, hncg) & result(has_converged) class(NumericalSolutionType) :: this !< NumericalSolutionType instance real(DP), intent(in) :: dxold_max !< the maximum dependent variable change for unrelaxed cells real(DP), intent(in) :: hncg !< largest dep. var. change at end of Picard iteration - real(DP), intent(in) :: dpak !< largest change in advanced packages logical(LGP) :: has_converged !< True, when converged has_converged = .false. if (abs(dxold_max) <= this%dvclose .and. & - abs(hncg) <= this%dvclose .and. & - abs(dpak) <= this%dvclose) then + abs(hncg) <= this%dvclose) then has_converged = .true. end if diff --git a/src/Solution/ParallelSolution.f90 b/src/Solution/ParallelSolution.f90 index b3ac1675623..1029a473ce7 100644 --- a/src/Solution/ParallelSolution.f90 +++ b/src/Solution/ParallelSolution.f90 @@ -1,6 +1,6 @@ module ParallelSolutionModule use KindModule, only: DP, LGP, I4B - use ConstantsModule, only: DONE, DZERO + use ConstantsModule, only: LENPAKLOC, DONE, DZERO use NumericalSolutionModule, only: NumericalSolutionType use mpi use MpiWorldModule @@ -13,6 +13,7 @@ module ParallelSolutionModule contains ! override procedure :: sln_has_converged => par_has_converged + procedure :: sln_package_convergence => par_package_convergence procedure :: sln_sync_newtonur_flag => par_sync_newtonur_flag procedure :: sln_nur_has_converged => par_nur_has_converged procedure :: sln_calc_ptc => par_calc_ptc @@ -47,6 +48,28 @@ function par_has_converged(this, max_dvc) result(has_converged) end function par_has_converged + function par_package_convergence(this, dpak, cpakout, iend) & + result(icnvg_global) + class(ParallelSolutionType) :: this !< parallel solution instance + real(DP), intent(in) :: dpak !< Newton Under-relaxation flag + character(len=LENPAKLOC), intent(in) :: cpakout + integer(I4B), intent(in) :: iend + ! local + integer(I4B) :: icnvg_global + integer(I4B) :: icnvg_local + integer :: ierr + type(MpiWorldType), pointer :: mpi_world + + mpi_world => get_mpi_world() + + icnvg_local = & + this%NumericalSolutionType%sln_package_convergence(dpak, cpakout, iend) + + call MPI_Allreduce(icnvg_local, icnvg_global, 1, MPI_INTEGER, & + MPI_MIN, mpi_world%comm, ierr) + + end function par_package_convergence + function par_sync_newtonur_flag(this, inewtonur) result(ivalue) class(ParallelSolutionType) :: this !< parallel solution instance integer(I4B), intent(in) :: inewtonur !< local Newton Under-relaxation flag @@ -61,12 +84,11 @@ function par_sync_newtonur_flag(this, inewtonur) result(ivalue) end function par_sync_newtonur_flag - function par_nur_has_converged(this, dxold_max, hncg, dpak) & + function par_nur_has_converged(this, dxold_max, hncg) & result(has_converged) class(ParallelSolutionType) :: this !< parallel solution instance real(DP), intent(in) :: dxold_max !< the maximum dependent variable change for cells not adjusted by NUR real(DP), intent(in) :: hncg !< largest dep. var. change at end of Picard iter. - real(DP), intent(in) :: dpak !< largest change in advanced packages logical(LGP) :: has_converged !< True, when converged ! local integer(I4B) :: icnvg_local @@ -79,7 +101,7 @@ function par_nur_has_converged(this, dxold_max, hncg, dpak) & has_converged = .false. icnvg_local = 0 if (this%NumericalSolutionType%sln_nur_has_converged( & - dxold_max, hncg, dpak)) then + dxold_max, hncg)) then icnvg_local = 1 end if From 2000548392718b940a84d4172abd65ffd44d04dc Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Fri, 28 Jul 2023 10:40:46 -0500 Subject: [PATCH 05/13] feat(sfr): add support for CELLIDs with zero values for unconnected reaches (#1311) * add deprecation warning for use of NONE for unconnected reaches * update mf6io with SFR changes * update ReleaseNotes with SFR changes * update error message in gwf3dis8 to use standard errmsg variable * update error message in gwf3dis8 to write a single error message per failure * standardize error messages in gwf3disv.f90 and gwf3disu.f90 * make specific test (test_gwf_sfr_conn) for sfr packagedata cellid for dis, disv, and disu * fix issue in mf6io SFR example * update mf6io SFR RBTH and RHK variable descriptions --- autotest/test_gwf_sfr_gwf_conn.py | 292 ++++++++++++++++++ doc/ReleaseNotes/develop.tex | 8 +- doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn | 6 +- .../mf6ivar/examples/gwf-sfr-example.dat | 76 ++--- doc/mf6io/mf6ivar/md/mf6ivar.md | 10 +- doc/mf6io/mf6ivar/tex/gwf-sfr-desc.tex | 6 +- src/Model/GroundWaterFlow/gwf3dis8.f90 | 97 +++--- src/Model/GroundWaterFlow/gwf3disu8.f90 | 27 +- src/Model/GroundWaterFlow/gwf3disv8.f90 | 130 ++++---- src/Model/GroundWaterFlow/gwf3sfr8.f90 | 28 +- 10 files changed, 503 insertions(+), 177 deletions(-) create mode 100644 autotest/test_gwf_sfr_gwf_conn.py diff --git a/autotest/test_gwf_sfr_gwf_conn.py b/autotest/test_gwf_sfr_gwf_conn.py new file mode 100644 index 00000000000..14f41968cb8 --- /dev/null +++ b/autotest/test_gwf_sfr_gwf_conn.py @@ -0,0 +1,292 @@ +import os + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation + +paktest = "sfr" +ex = [ + "sfr_dis", + "sfr_dis_fail", + "sfr_dis_none", + "sfr_disv", + "sfr_disv_fail", + "sfr_disv_none", + "sfr_disu", + "sfr_disu_fail", + "sfr_disu_none", +] +dis_types = [ + "dis", + "dis", + "dis", + "disv", + "disv", + "disv", + "disu", + "disu", + "disu", +] + +# spatial discretization data +nlay, nrow, ncol = 1, 1, 1 +delr, delc = 100.0, 100.0 +top = 0.0 +botm = -10.0 +strt = 0.0 + +# spatial discretization data for disv and disu +vertices = [(0, 0.0, 0.0), (1, 0.0, delc), (2, delr, delc), (3, delr, 0.0)] +cell2d = [(0, delr / 2.0, delc / 2.0, 4, 0, 1, 2, 3)] + +# sfr data +nreaches = 10 +rlen = 10.0 +rwid = 10.0 +roughness = 0.001 +rbth = 1.0 +rhk = 0.0 +slope = 0.001 +ustrf = 1.0 +ndv = 0 + + +def build_model(idx, ws): + # static model data + # temporal discretization + nper = 1 + tdis_rc = [(11.0, 11, 1.0)] + ts_times = np.arange(0.0, 12.0, 1.0, dtype=float) + ts_flows = np.array([1000.0] + [float(q) for q in range(1000, -100, -100)]) + + # build MODFLOW 6 files + name = ex[idx] + dis_type = dis_types[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, + version="mf6", + exe_name="mf6", + sim_ws=ws, + ) + sim.simulation_data.verify_data = False + + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, + time_units="seconds", + nper=nper, + perioddata=tdis_rc, + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="ALL", + ) + + # create gwf model + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=name, + ) + + if dis_type == "dis": + dis = flopy.mf6.ModflowGwfdis( + gwf, + length_units="meters", + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=1, + ) + elif dis_type == "disv": + dis = flopy.mf6.ModflowGwfdisv( + gwf, + length_units="meters", + nlay=nlay, + ncpl=1, + nvert=4, + vertices=vertices, + cell2d=cell2d, + top=top, + botm=botm, + idomain=1, + ) + else: + disukwargs = flopy.utils.gridutil.get_disu_kwargs( + nlay, + nrow, + ncol, + [delr], + [delc], + top, + [botm], + ) + dis = flopy.mf6.ModflowGwfdisu( + gwf, + vertices=vertices, + cell2d=cell2d, + **disukwargs, + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf(gwf) + + # chd files + # chd data + if dis_type == "dis": + spd = [ + [(0, 0, 0), 0.0], + ] + elif dis_type == "disv": + spd = [ + [(0, 0), 0.0], + ] + else: + spd = [ + [(0,), 0.0], + ] + + chd = flopy.mf6.modflow.ModflowGwfchd( + gwf, maxbound=1, stress_period_data=spd, pname="chd-1" + ) + + # sfr file + if dis_type == "dis": + if "fail" in name: + cellid = (2, 2, 2) + elif "none" in name: + cellid = "none" + else: + cellid = (-1, -1, -1) + elif dis_type == "disv": + if "fail" in name: + cellid = (2, 2) + elif "none" in name: + cellid = "none" + else: + cellid = (-1, -1) + else: + if "fail" in name: + cellid = (2,) + elif "none" in name: + cellid = "none" + else: + cellid = (-1,) + packagedata = [] + for irch in range(nreaches): + nconn = 1 + if 0 < irch < nreaches - 1: + nconn += 1 + rp = [ + irch, + cellid, + rlen, + rwid, + slope, + top, + rbth, + rhk, + roughness, + nconn, + ustrf, + ndv, + ] + packagedata.append(rp) + + if not ws.endswith("mf6"): + packagedata = packagedata[::-1] + + connectiondata = [] + inflow_loc = 0 + ioutflow_loc = nreaches - 1 + for irch in range(nreaches): + rc = [irch] + if irch > 0: + rc.append(irch - 1) + if irch < nreaches - 1: + rc.append(-(irch + 1)) + connectiondata.append(rc) + + ts_names = ["inflow"] + perioddata = [ + [inflow_loc, "inflow", "inflow"], + ] + ts_methods = ["linearend"] * len(ts_names) + ts_data = [] + for t, q in zip(ts_times, ts_flows): + ts_data.append((t, q)) + + sfr = flopy.mf6.ModflowGwfsfr( + gwf, + print_stage=True, + print_flows=True, + print_input=True, + mover=True, + nreaches=nreaches, + packagedata=packagedata, + connectiondata=connectiondata, + perioddata=perioddata, + pname="sfr-1", + ) + fname = f"{name}.sfr.ts" + sfr.ts.initialize( + filename=fname, + timeseries=ts_data, + time_series_namerecord=ts_names, + interpolation_methodrecord=ts_methods, + ) + fname = f"{name}.sfr.obs" + sfr_obs = { + f"{fname}.csv": [ + ("inflow", "ext-inflow", (inflow_loc,)), + ("outflow", "ext-outflow", (ioutflow_loc,)), + ] + } + sfr.obs.initialize(filename=fname, print_input=True, continuous=sfr_obs) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + printrecord=[ + ("BUDGET", "ALL"), + ], + ) + + return sim + + +def build_models(idx, base_ws): + sim = build_model(idx, base_ws) + return sim, None + + +@pytest.mark.parametrize("idx, name", enumerate(ex)) +def test_mf6model(idx, name, function_tmpdir, targets): + if "fail" in name: + require_failure = True + else: + require_failure = False + ws = str(function_tmpdir) + test = TestFramework() + test.build(build_models, idx, ws) + test.run( + TestSimulation( + name=name, + exe_dict=targets, + idxsim=idx, + mf6_regression=False, + make_comparison=False, + require_failure=require_failure, + ), + ws, + ) diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index 9dab6e9899a..b1fcec5c8c8 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -39,12 +39,12 @@ % \item xxx %\end{itemize} - %\underline{ADVANCED STRESS PACKAGES} - %\begin{itemize} - % \item xxx + \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 unconnected reaches. Warning messages will be issued if NONE is specified for unconnected reaches. Specifying a CELLID of NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error. % \item xxx % \item xxx - %\end{itemize} + \end{itemize} %\underline{SOLUTION} %\begin{itemize} diff --git a/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn b/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn index 5a23b919e45..d5d5313db3b 100644 --- a/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwf-sfr.dfn @@ -367,7 +367,7 @@ tagged false in_record true reader urword longname cell identifier -description The keyword `NONE' must be specified for reaches that are not connected to an underlying GWF cell. The keyword `NONE' is used for reaches that are in cells that have IDOMAIN values less than one or are in areas not covered by the GWF model grid. Reach-aquifer flow is not calculated if the keyword `NONE' is specified. +description is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. For reaches that are not connected to an underlying GWF cell, a zero should be specified for each grid dimension. For example, for a DIS grid a CELLID of 0 0 0 should be specified. Reach-aquifer flow is not calculated for unconnected reaches. The keyword NONE can be still be specified to identify unconnected reaches for backward compatibility with previous versions of MODFLOW 6 but eventually NONE will be deprecated and will cause MODFLOW 6 to terminate with an error. block packagedata name rlen @@ -417,7 +417,7 @@ tagged false in_record true reader urword longname streambed thickness -description real value that defines the thickness of the reach streambed. RBTH can be any value if CELLID is `NONE'. Otherwise, RBTH must be greater than zero. +description real value that defines the thickness of the reach streambed. RBTH can be any value if the reach is not connected to an underlying GWF cell. Otherwise, RBTH must be greater than zero. block packagedata name rhk @@ -427,7 +427,7 @@ tagged false in_record true reader urword longname -description real value that defines the hydraulic conductivity of the reach streambed. RHK can be any positive value if CELLID is `NONE'. Otherwise, RHK must be greater than zero. +description real value that defines the hydraulic conductivity of the reach streambed. RHK can be any positive value if the reach is not connected to an underlying GWF cell. Otherwise, RHK must be greater than zero. block packagedata name man diff --git a/doc/mf6io/mf6ivar/examples/gwf-sfr-example.dat b/doc/mf6io/mf6ivar/examples/gwf-sfr-example.dat index 18d885dfd4c..eb7791bac41 100644 --- a/doc/mf6io/mf6ivar/examples/gwf-sfr-example.dat +++ b/doc/mf6io/mf6ivar/examples/gwf-sfr-example.dat @@ -13,44 +13,44 @@ BEGIN DIMENSIONS END DIMENSIONS BEGIN PACKAGEDATA -#rno k i j rlen rwid rgrd rtp rbth rhk man ncon ustrf ndv boundname - 1 1 1 1 4500. 12 8.67E-04 1093.048 3.0 0.00003 0.03 1 1.0 0 reach1 - 2 1 2 2 7000. 12 8.67E-04 1088.059 3.0 0.00003 0.03 2 1.0 0 reach2 - 3 1 3 3 6000. 12 8.67E-04 1082.419 3.0 0.00003 0.03 2 1.0 0 reach3 - 4 1 3 4 5550. 12 8.67E-04 1077.408 3.0 0.00003 0.03 3 1.0 1 reach4 - 5 1 4 5 6500. 12 9.43E-04 1071.934 3.0 0.00003 0.03 2 1.0 0 - 6 1 5 6 5000. 12 9.43E-04 1066.509 3.0 0.00003 0.03 2 1.0 0 - 7 1 6 6 5000. 12 9.43E-04 1061.792 3.0 0.00003 0.03 2 1.0 0 - 8 1 7 6 5000. 12 9.43E-04 1057.075 3.0 0.00003 0.03 2 1.0 0 - 9 1 8 6 5000. 12 9.43E-04 1052.359 3.0 0.00003 0.03 2 1.0 0 - 10 1 3 5 5000. 10 5.45E-04 1073.636 2.0 0.00003 0.03 2 0.0 0 canal - 11 1 3 6 5000. 10 5.45E-04 1070.909 2.0 0.00003 0.03 2 1.0 0 canal - 12 1 3 7 4500. 10 5.45E-04 1068.318 2.0 0.00003 0.03 2 1.0 0 canal - 13 1 4 8 6000. 10 5.45E-04 1065.455 2.0 0.00003 0.03 2 1.0 0 canal - 14 1 5 8 5000. 10 5.45E-04 1062.455 2.0 0.00003 0.03 2 1.0 0 canal - 15 1 6 8 2000. 10 5.45E-04 1060.545 2.0 0.00003 0.03 2 1.0 0 canal - 16 1 510 2500. 10 1.81E-03 1077.727 3.0 0.00003 0.03 1 1.0 0 - 17 1 5 9 5000. 10 1.81E-03 1070.909 3.0 0.00003 0.03 2 1.0 0 - 18 1 6 8 3500. 10 1.81E-03 1063.182 3.0 0.00003 0.03 2 1.0 0 - 19 1 6 8 4000. 15 1.00E-03 1058.000 3.0 0.00003 0.03 3 1.0 0 - 20 1 7 7 5000. 15 1.00E-03 1053.500 3.0 0.00003 0.03 2 1.0 0 - 21 1 8 7 3500. 15 1.00E-03 1049.250 3.0 0.00003 0.03 2 1.0 0 - 22 1 8 6 2500. 15 1.00E-03 1046.250 3.0 0.00003 0.03 2 1.0 0 - 23 1 9 6 5000. 12 9.09E-04 1042.727 3.0 0.00003 0.03 3 1.0 0 - 24 1 10 7 5000. 12 9.09E-04 1038.182 3.0 0.00003 0.03 2 1.0 0 - 25 1 11 7 5000. 12 9.09E-04 1033.636 3.0 0.00003 0.03 2 1.0 0 - 26 1 12 7 5000. 12 9.09E-04 1029.091 3.0 0.00003 0.03 2 1.0 0 - 27 1 13 7 2000. 12 9.09E-04 1025.909 3.0 0.00003 0.03 2 1.0 0 - 28 1 14 9 5000. 55 9.67E-04 1037.581 3.0 0.00006 0.025 1 1.0 0 - 29 1 13 8 5500. 55 9.67E-04 1032.500 3.0 0.00006 0.025 2 1.0 0 - 30 1 13 7 5000. 55 9.67E-04 1027.419 3.0 0.00006 0.025 2 1.0 0 - 31 1 13 6 5000. 40 1.25E-03 1021.875 3.0 0.00006 0.025 3 1.0 0 - 32 1 13 5 5000. 40 1.25E-03 1015.625 3.0 0.00006 0.025 2 1.0 0 - 33 1 13 4 5000. 40 1.25E-03 1009.375 3.0 0.00006 0.025 2 1.0 0 - 34 1 13 3 5000. 40 1.25E-03 1003.125 3.0 0.00006 0.025 2 1.0 0 - 35 1 13 2 5000. 40 1.25E-03 996.8750 3.0 0.00006 0.025 2 1.0 0 - 36 1 13 1 3000. 40 1.25E-03 991.8750 3.0 0.00006 0.025 2 1.0 0 - 37 none 5000. 40 1.25E-03 985.6250 3.0 0.00006 0.025 1 1.0 0 +#rno k i j rlen rwid rgrd rtp rbth rhk man ncon ustrf ndv boundname + 1 1 1 1 4500. 12 8.67E-04 1093.048 3.0 0.00003 0.03 1 1.0 0 reach1 + 2 1 2 2 7000. 12 8.67E-04 1088.059 3.0 0.00003 0.03 2 1.0 0 reach2 + 3 1 3 3 6000. 12 8.67E-04 1082.419 3.0 0.00003 0.03 2 1.0 0 reach3 + 4 1 3 4 5550. 12 8.67E-04 1077.408 3.0 0.00003 0.03 3 1.0 1 reach4 + 5 1 4 5 6500. 12 9.43E-04 1071.934 3.0 0.00003 0.03 2 1.0 0 + 6 1 5 6 5000. 12 9.43E-04 1066.509 3.0 0.00003 0.03 2 1.0 0 + 7 1 6 6 5000. 12 9.43E-04 1061.792 3.0 0.00003 0.03 2 1.0 0 + 8 1 7 6 5000. 12 9.43E-04 1057.075 3.0 0.00003 0.03 2 1.0 0 + 9 1 8 6 5000. 12 9.43E-04 1052.359 3.0 0.00003 0.03 2 1.0 0 + 10 1 3 5 5000. 10 5.45E-04 1073.636 2.0 0.00003 0.03 2 0.0 0 canal + 11 1 3 6 5000. 10 5.45E-04 1070.909 2.0 0.00003 0.03 2 1.0 0 canal + 12 1 3 7 4500. 10 5.45E-04 1068.318 2.0 0.00003 0.03 2 1.0 0 canal + 13 1 4 8 6000. 10 5.45E-04 1065.455 2.0 0.00003 0.03 2 1.0 0 canal + 14 1 5 8 5000. 10 5.45E-04 1062.455 2.0 0.00003 0.03 2 1.0 0 canal + 15 1 6 8 2000. 10 5.45E-04 1060.545 2.0 0.00003 0.03 2 1.0 0 canal + 16 1 5 10 2500. 10 1.81E-03 1077.727 3.0 0.00003 0.03 1 1.0 0 + 17 1 5 9 5000. 10 1.81E-03 1070.909 3.0 0.00003 0.03 2 1.0 0 + 18 1 6 8 3500. 10 1.81E-03 1063.182 3.0 0.00003 0.03 2 1.0 0 + 19 1 6 8 4000. 15 1.00E-03 1058.000 3.0 0.00003 0.03 3 1.0 0 + 20 1 7 7 5000. 15 1.00E-03 1053.500 3.0 0.00003 0.03 2 1.0 0 + 21 1 8 7 3500. 15 1.00E-03 1049.250 3.0 0.00003 0.03 2 1.0 0 + 22 1 8 6 2500. 15 1.00E-03 1046.250 3.0 0.00003 0.03 2 1.0 0 + 23 1 9 6 5000. 12 9.09E-04 1042.727 3.0 0.00003 0.03 3 1.0 0 + 24 1 10 7 5000. 12 9.09E-04 1038.182 3.0 0.00003 0.03 2 1.0 0 + 25 1 11 7 5000. 12 9.09E-04 1033.636 3.0 0.00003 0.03 2 1.0 0 + 26 1 12 7 5000. 12 9.09E-04 1029.091 3.0 0.00003 0.03 2 1.0 0 + 27 1 13 7 2000. 12 9.09E-04 1025.909 3.0 0.00003 0.03 2 1.0 0 + 28 1 14 9 5000. 55 9.67E-04 1037.581 3.0 0.00006 0.025 1 1.0 0 + 29 1 13 8 5500. 55 9.67E-04 1032.500 3.0 0.00006 0.025 2 1.0 0 + 30 1 13 7 5000. 55 9.67E-04 1027.419 3.0 0.00006 0.025 2 1.0 0 + 31 1 13 6 5000. 40 1.25E-03 1021.875 3.0 0.00006 0.025 3 1.0 0 + 32 1 13 5 5000. 40 1.25E-03 1015.625 3.0 0.00006 0.025 2 1.0 0 + 33 1 13 4 5000. 40 1.25E-03 1009.375 3.0 0.00006 0.025 2 1.0 0 + 34 1 13 3 5000. 40 1.25E-03 1003.125 3.0 0.00006 0.025 2 1.0 0 + 35 1 13 2 5000. 40 1.25E-03 996.8750 3.0 0.00006 0.025 2 1.0 0 + 36 1 13 1 3000. 40 1.25E-03 991.8750 3.0 0.00006 0.025 2 1.0 0 + 37 0 0 0 5000. 40 1.25E-03 985.6250 3.0 0.00006 0.025 1 1.0 0 END PACKAGEDATA BEGIN CONNECTIONDATA diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 796bc5bcfeb..0bce538bf7a 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -576,13 +576,13 @@ | GWF | SFR | OPTIONS | TIME_CONVERSION | DOUBLE PRECISION | real value that is used to convert user-specified Manning's roughness coefficients from seconds to model time units. TIME\_CONVERSION should be set to 1.0, 60.0, 3,600.0, 86,400.0, and 31,557,600.0 when using time units (TIME\_UNITS) of seconds, minutes, hours, days, or years in the simulation, respectively. TIME\_CONVERSION does not need to be specified if TIME\_UNITS are seconds. | | GWF | SFR | DIMENSIONS | NREACHES | INTEGER | integer value specifying the number of stream reaches. There must be NREACHES entries in the PACKAGEDATA block. | | GWF | SFR | PACKAGEDATA | RNO | INTEGER | integer value that defines the reach number associated with the specified PACKAGEDATA data on the line. RNO must be greater than zero and less than or equal to NREACHES. Reach information must be specified for every reach or the program will terminate with an error. The program will also terminate with an error if information for a reach is specified more than once. | -| GWF | SFR | PACKAGEDATA | CELLID | INTEGER (NCELLDIM) | The keyword `NONE' must be specified for reaches that are not connected to an underlying GWF cell. The keyword `NONE' is used for reaches that are in cells that have IDOMAIN values less than one or are in areas not covered by the GWF model grid. Reach-aquifer flow is not calculated if the keyword `NONE' is specified. | +| GWF | SFR | PACKAGEDATA | CELLID | INTEGER (NCELLDIM) | is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. For reaches that are not connected to an underlying GWF cell, a zero should be specified for each grid dimension. For example, for a DIS grid a CELLID of 0 0 0 should be specified. Reach-aquifer flow is not calculated for unconnected reaches. The keyword NONE can be still be specified to identify unconnected reaches for backward compatibility with previous versions of MODFLOW 6 but eventually NONE will be deprecated and will cause MODFLOW 6 to terminate with an error. | | GWF | SFR | PACKAGEDATA | RLEN | DOUBLE PRECISION | real value that defines the reach length. RLEN must be greater than zero. | | GWF | SFR | PACKAGEDATA | RWID | DOUBLE PRECISION | real value that defines the reach width. RWID must be greater than zero. | | GWF | SFR | PACKAGEDATA | RGRD | DOUBLE PRECISION | real value that defines the stream gradient (slope) across the reach. RGRD must be greater than zero. | | GWF | SFR | PACKAGEDATA | RTP | DOUBLE PRECISION | real value that defines the bottom elevation of the reach. | -| GWF | SFR | PACKAGEDATA | RBTH | DOUBLE PRECISION | real value that defines the thickness of the reach streambed. RBTH can be any value if CELLID is `NONE'. Otherwise, RBTH must be greater than zero. | -| GWF | SFR | PACKAGEDATA | RHK | DOUBLE PRECISION | real value that defines the hydraulic conductivity of the reach streambed. RHK can be any positive value if CELLID is `NONE'. Otherwise, RHK must be greater than zero. | +| GWF | SFR | PACKAGEDATA | RBTH | DOUBLE PRECISION | real value that defines the thickness of the reach streambed. RBTH can be any value if the reach is not connected to an underlying GWF cell. Otherwise, RBTH must be greater than zero. | +| GWF | SFR | PACKAGEDATA | RHK | DOUBLE PRECISION | real value that defines the hydraulic conductivity of the reach streambed. RHK can be any positive value if the reach is not connected to an underlying GWF cell. Otherwise, RHK must be greater than zero. | | GWF | SFR | PACKAGEDATA | MAN | STRING | real or character value that defines the Manning's roughness coefficient for the reach. MAN must be greater than zero. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | | GWF | SFR | PACKAGEDATA | NCON | INTEGER | integer value that defines the number of reaches connected to the reach. If a value of zero is specified for NCON an entry for RNO is still required in the subsequent CONNECTIONDATA block. | | GWF | SFR | PACKAGEDATA | USTRF | DOUBLE PRECISION | real value that defines the fraction of upstream flow from each upstream reach that is applied as upstream inflow to the reach. The sum of all USTRF values for all reaches connected to the same upstream reach must be equal to one and USTRF must be greater than or equal to zero. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | @@ -938,7 +938,7 @@ | GWT | SSM | OPTIONS | PRINT_FLOWS | KEYWORD | keyword to indicate that the list of SSM flow rates will be printed to the listing file for every stress period time step in which ``BUDGET PRINT'' is specified in Output Control. If there is no Output Control option and ``PRINT\_FLOWS'' is specified, then flow rates are printed for the last time step of each stress period. | | GWT | SSM | OPTIONS | SAVE_FLOWS | KEYWORD | keyword to indicate that SSM flow terms will be written to the file specified with ``BUDGET FILEOUT'' in Output Control. | | GWT | SSM | SOURCES | PNAME | STRING | name of the flow package for which an auxiliary variable contains a source concentration. If this flow package is represented using an advanced transport package (SFT, LKT, MWT, or UZT), then the advanced transport package will override SSM terms specified here. | -| GWT | SSM | SOURCES | SRCTYPE | STRING | keyword indicating how concentration will be assigned for sources and sinks. Keyword must be specified as either AUX or AUXMIXED. For both options the user must provide an auxiliary variable in the corresponding flow package. The auxiliary variable must have the same name as the AUXNAME value that follows. If the AUX keyword is specified, then the auxiliary variable specified by the user will be assigned as the concenration value for groundwater sources (flows with a positive sign). For negative flow rates (sinks), groundwater will be withdrawn from the cell at the simulated concentration of the cell. The AUXMIXED option provides an alternative method for how to determine the concentration of sinks. If the cell concentration is larger than the user-specified auxiliary concentration, then the concentration of groundwater withdrawn from the cell will be assigned as the user-specified concentration. Alternatively, if the user-specified auxiliary concentration is larger than the cell concentration, then groundwater will be withdrawn at the cell concentration. Thus, the AUXMIXED option is designed to work with the Evapotranspiration (EVT) and Recharge (RCH) Packages where water may be withdrawn at a concentration that is less than the cell concentration. | +| GWT | SSM | SOURCES | SRCTYPE | STRING | keyword indicating how concentration will be assigned for sources and sinks. Keyword must be specified as either AUX or AUXMIXED. For both options the user must provide an auxiliary variable in the corresponding flow package. The auxiliary variable must have the same name as the AUXNAME value that follows. If the AUX keyword is specified, then the auxiliary variable specified by the user will be assigned as the concentration value for groundwater sources (flows with a positive sign). For negative flow rates (sinks), groundwater will be withdrawn from the cell at the simulated concentration of the cell. The AUXMIXED option provides an alternative method for how to determine the concentration of sinks. If the cell concentration is larger than the user-specified auxiliary concentration, then the concentration of groundwater withdrawn from the cell will be assigned as the user-specified concentration. Alternatively, if the user-specified auxiliary concentration is larger than the cell concentration, then groundwater will be withdrawn at the cell concentration. Thus, the AUXMIXED option is designed to work with the Evapotranspiration (EVT) and Recharge (RCH) Packages where water may be withdrawn at a concentration that is less than the cell concentration. | | GWT | SSM | SOURCES | AUXNAME | STRING | name of the auxiliary variable in the package PNAME. This auxiliary variable must exist and be specified by the user in that package. The values in this auxiliary variable will be used to set the concentration associated with the flows for that boundary package. | | GWT | SSM | FILEINPUT | PNAME | STRING | name of the flow package for which an SPC6 input file contains a source concentration. If this flow package is represented using an advanced transport package (SFT, LKT, MWT, or UZT), then the advanced transport package will override SSM terms specified here. | | GWT | SSM | FILEINPUT | SPC6 | KEYWORD | keyword to specify that record corresponds to a source sink mixing input file. | @@ -1094,7 +1094,7 @@ | GWT | MWT | PACKAGEDATA | BOUNDNAME | STRING | name of the well 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. | | GWT | MWT | PERIOD | IPER | INTEGER | integer value specifying the starting stress period number for which the data specified in the PERIOD block apply. IPER must be less than or equal to NPER in the TDIS Package and greater than zero. The IPER value assigned to a stress period block must be greater than the IPER value assigned for the previous PERIOD block. The information specified in the PERIOD block will continue to apply for all subsequent stress periods, unless the program encounters another PERIOD block. | | GWT | MWT | PERIOD | MAWNO | INTEGER | integer value that defines the well number associated with the specified PERIOD data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. | -| GWT | MWT | PERIOD | MWTSETTING | KEYSTRING | line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWTSETTING string include: STATUS, CONCENTRATION, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the concentration of associated with the corresponding flow terms. Concentrations cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated concentration of the well. | +| GWT | MWT | PERIOD | MWTSETTING | KEYSTRING | line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWTSETTING string include: STATUS, CONCENTRATION, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the concentration associated with the corresponding flow terms. Concentrations cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated concentration of the well. | | GWT | MWT | PERIOD | STATUS | STRING | keyword option to define well status. STATUS can be ACTIVE, INACTIVE, or CONSTANT. By default, STATUS is ACTIVE, which means that concentration will be calculated for the well. If a well is inactive, then there will be no solute mass fluxes into or out of the well and the inactive value will be written for the well concentration. If a well is constant, then the concentration for the well will be fixed at the user specified value. | | GWT | MWT | PERIOD | CONCENTRATION | STRING | real or character value that defines the concentration for the well. The specified CONCENTRATION is only applied if the well is a constant concentration well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | | GWT | MWT | PERIOD | RATE | STRING | real or character value that defines the injection solute concentration $(ML^{-3})$ for the well. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value. | diff --git a/doc/mf6io/mf6ivar/tex/gwf-sfr-desc.tex b/doc/mf6io/mf6ivar/tex/gwf-sfr-desc.tex index b0664a7f91f..99738579c0a 100644 --- a/doc/mf6io/mf6ivar/tex/gwf-sfr-desc.tex +++ b/doc/mf6io/mf6ivar/tex/gwf-sfr-desc.tex @@ -67,7 +67,7 @@ \begin{description} \item \texttt{rno}---integer value that defines the reach number associated with the specified PACKAGEDATA data on the line. RNO must be greater than zero and less than or equal to NREACHES. Reach information must be specified for every reach or the program will terminate with an error. The program will also terminate with an error if information for a reach is specified more than once. -\item \texttt{cellid}---The keyword `NONE' must be specified for reaches that are not connected to an underlying GWF cell. The keyword `NONE' is used for reaches that are in cells that have IDOMAIN values less than one or are in areas not covered by the GWF model grid. Reach-aquifer flow is not calculated if the keyword `NONE' is specified. +\item \texttt{cellid}---is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. For reaches that are not connected to an underlying GWF cell, a zero should be specified for each grid dimension. For example, for a DIS grid a CELLID of 0 0 0 should be specified. Reach-aquifer flow is not calculated for unconnected reaches. The keyword NONE can be still be specified to identify unconnected reaches for backward compatibility with previous versions of MODFLOW 6 but eventually NONE will be deprecated and will cause MODFLOW 6 to terminate with an error. \item \texttt{rlen}---real value that defines the reach length. RLEN must be greater than zero. @@ -77,9 +77,9 @@ \item \texttt{rtp}---real value that defines the bottom elevation of the reach. -\item \texttt{rbth}---real value that defines the thickness of the reach streambed. RBTH can be any value if CELLID is `NONE'. Otherwise, RBTH must be greater than zero. +\item \texttt{rbth}---real value that defines the thickness of the reach streambed. RBTH can be any value if the reach is not connected to an underlying GWF cell. Otherwise, RBTH must be greater than zero. -\item \texttt{rhk}---real value that defines the hydraulic conductivity of the reach streambed. RHK can be any positive value if CELLID is `NONE'. Otherwise, RHK must be greater than zero. +\item \texttt{rhk}---real value that defines the hydraulic conductivity of the reach streambed. RHK can be any positive value if the reach is not connected to an underlying GWF cell. Otherwise, RHK must be greater than zero. \item \textcolor{blue}{\texttt{man}---real or character value that defines the Manning's roughness coefficient for the reach. MAN must be greater than zero. If the Options block includes a TIMESERIESFILE entry (see the ``Time-Variable Input'' section), values can be obtained from a time series by entering the time-series name in place of a numeric value.} diff --git a/src/Model/GroundWaterFlow/gwf3dis8.f90 b/src/Model/GroundWaterFlow/gwf3dis8.f90 index f13f34ee9f9..0d3e1d334c7 100644 --- a/src/Model/GroundWaterFlow/gwf3dis8.f90 +++ b/src/Model/GroundWaterFlow/gwf3dis8.f90 @@ -8,6 +8,7 @@ module GwfDisModule ubdsv06 use SimModule, only: count_errors, store_error, store_error_unit, & store_error_filename + use SimVariablesModule, only: errmsg use MemoryManagerModule, only: mem_allocate use TdisModule, only: kstp, kper, pertim, totim, delt @@ -420,7 +421,6 @@ subroutine grid_finalize(this) ! -- dummy class(GwfDisType) :: this ! -- locals - character(len=300) :: ermsg integer(I4B) :: n, i, j, k integer(I4B) :: node integer(I4B) :: noder @@ -469,8 +469,8 @@ subroutine grid_finalize(this) dz = top - this%bot3d(j, i, k) if (dz <= DZERO) then n = n + 1 - write (ermsg, fmt=fmtdz) k, i, j, top, this%bot3d(j, i, k) - call store_error(ermsg) + write (errmsg, fmt=fmtdz) k, i, j, top, this%bot3d(j, i, k) + call store_error(errmsg) end if end do end do @@ -749,7 +749,6 @@ subroutine nodeu_to_array(this, nodeu, arr) integer(I4B), intent(in) :: nodeu integer(I4B), dimension(:), intent(inout) :: arr ! -- local - character(len=LINELENGTH) :: errmsg integer(I4B) :: isize integer(I4B) :: i, j, k ! ------------------------------------------------------------------------------ @@ -792,7 +791,6 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck ! -- local - character(len=LINELENGTH) :: errmsg ! ------------------------------------------------------------------------------ ! ! -- check the node number if requested @@ -800,8 +798,9 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) ! ! -- If within valid range, convert to reduced nodenumber if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (errmsg, '(a,i10)') & - 'Nodenumber less than 1 or greater than nodes:', nodeu + write (errmsg, '(a,i0,a)') & + 'Node number (', nodeu, & + ') less than 1 or greater than the number of nodes.' call store_error(errmsg) nodenumber = 0 else @@ -836,7 +835,6 @@ function get_nodenumber_idx3(this, k, i, j, icheck) & integer(I4B), intent(in) :: k, i, j integer(I4B), intent(in) :: icheck ! -- local - character(len=LINELENGTH) :: errmsg integer(I4B) :: nodeu ! formats character(len=*), parameter :: fmterr = & @@ -864,8 +862,8 @@ function get_nodenumber_idx3(this, k, i, j, icheck) & ! ! -- Error if outside of range if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (errmsg, '(a,i10)') & - 'Nodenumber less than 1 or greater than nodes:', nodeu + write (errmsg, '(a,i0,a)') & + 'Node number (', nodeu, ')less than 1 or greater than nodes.' call store_error(errmsg) end if end if @@ -967,7 +965,6 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & integer(I4B) :: k, i, j, nlay, nrow, ncol integer(I4B) :: lloclocal, ndum, istat, n real(DP) :: r - character(len=LINELENGTH) :: ermsg, fname ! ------------------------------------------------------------------------------ ! if (present(flag_string)) then @@ -1001,28 +998,35 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & end if end if ! + errmsg = "" + ! if (k < 1 .or. k > nlay) then - write (ermsg, *) ' Layer number in list is outside of the grid', k - call store_error(ermsg) + write (errmsg, '(a,i0,a)') & + 'Layer number in list (', k, ') is outside of the grid.' end if if (i < 1 .or. i > nrow) then - write (ermsg, *) ' Row number in list is outside of the grid', i - call store_error(ermsg) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Row number in list (', i, & + ') is outside of the grid.' end if if (j < 1 .or. j > ncol) then - write (ermsg, *) ' Column number in list is outside of the grid', j - call store_error(ermsg) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Column number in list (', j, & + ') is outside of the grid.' end if + ! nodeu = get_node(k, i, j, nlay, nrow, ncol) ! if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (ermsg, *) ' Node number in list is outside of the grid', nodeu - call store_error(ermsg) - inquire (unit=in, name=fname) - call store_error('Error converting in file: ') - call store_error(trim(adjustl(fname))) - call store_error('Cell number cannot be determined in line: ') - call store_error(trim(adjustl(line))) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), & + "Node number in list (", nodeu, ") is outside of the grid. "// & + "Cell number cannot be determined in line '"// & + trim(adjustl(line))//"'." + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then + call store_error(errmsg) call store_error_unit(in) end if ! @@ -1060,7 +1064,6 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & integer(I4B) :: k, i, j, nlay, nrow, ncol integer(I4B) :: istat real(DP) :: r - character(len=LINELENGTH) :: ermsg, fname ! ------------------------------------------------------------------------------ ! if (present(flag_string)) then @@ -1095,28 +1098,35 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & end if end if ! + errmsg = "" + ! if (k < 1 .or. k > nlay) then - write (ermsg, *) ' Layer number in list is outside of the grid', k - call store_error(ermsg) + write (errmsg, '(a,i0,a)') & + 'Layer number in list (', k, ') is outside of the grid.' end if if (i < 1 .or. i > nrow) then - write (ermsg, *) ' Row number in list is outside of the grid', i - call store_error(ermsg) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Row number in list (', i, & + ') is outside of the grid.' end if if (j < 1 .or. j > ncol) then - write (ermsg, *) ' Column number in list is outside of the grid', j - call store_error(ermsg) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Column number in list (', j, & + ') is outside of the grid.' end if + ! nodeu = get_node(k, i, j, nlay, nrow, ncol) ! if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (ermsg, *) ' Node number in list is outside of the grid', nodeu - call store_error(ermsg) - inquire (unit=inunit, name=fname) - call store_error('Error converting in file: ') - call store_error(trim(adjustl(fname))) - call store_error('Cell number cannot be determined in cellid: ') - call store_error(trim(adjustl(cellid))) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), & + "Cell number cannot be determined for cellid ("// & + trim(adjustl(cellid))//") and results in a user "// & + "node number (", nodeu, ") that is outside of the grid." + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then + call store_error(errmsg) call store_error_unit(inunit) end if ! @@ -1674,7 +1684,6 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr - character(len=LINELENGTH) :: errmsg ! ------------------------------------------------------------------------------ ! ! -- set variables @@ -1696,7 +1705,7 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & nodeu = get_node(1, ir, ic, nlay, nrow, ncol) il = this%ibuff(nodeu) if (il < 1 .or. il > nlay) then - write (errmsg, *) 'INVALID LAYER NUMBER: ', il + write (errmsg, '(a,1x,i0)') 'Invalid layer number:', il call store_error(errmsg, terminate=.TRUE.) end if nodeu = get_node(il, ir, ic, nlay, nrow, ncol) @@ -1713,9 +1722,9 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & ! -- Check for errors nbound = ipos - 1 if (ierr > 0) then - write (errmsg, *) 'MAXBOUND DIMENSION IS TOO SMALL.' - call store_error(errmsg) - write (errmsg, *) 'INCREASE MAXBOUND TO: ', ierr + write (errmsg, '(a,1x,i0)') & + 'MAXBOUND dimension is too small.'// & + 'INCREASE MAXBOUND TO:', ierr call store_error(errmsg, terminate=.TRUE.) end if ! @@ -1732,7 +1741,7 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & call ReadArray(inunit, nodelist, aname, this%ndim, maxbnd, iout, 0) do noder = 1, maxbnd if (noder < 1 .or. noder > this%nodes) then - write (errmsg, *) 'INVALID NODE NUMBER: ', noder + write (errmsg, '(a,1x,i0)') 'Invalid node number:', noder call store_error(errmsg, terminate=.TRUE.) end if end do diff --git a/src/Model/GroundWaterFlow/gwf3disu8.f90 b/src/Model/GroundWaterFlow/gwf3disu8.f90 index 003ba6091c8..404ea357340 100644 --- a/src/Model/GroundWaterFlow/gwf3disu8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disu8.f90 @@ -1156,8 +1156,9 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) ! if (icheck /= 0) then if (nodeu < 1 .or. nodeu > this%nodes) then - write (errmsg, '(a,i10)') & - 'Nodenumber less than 1 or greater than nodes:', nodeu + write (errmsg, '(a,i0,a,i0,a)') & + 'Node number (', nodeu, ') is less than 1 or greater than nodes (', & + this%nodes, ').' call store_error(errmsg) end if end if @@ -1416,7 +1417,6 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & ! -- local integer(I4B) :: lloclocal, ndum, istat, n real(DP) :: r - character(len=LINELENGTH) :: fname ! ------------------------------------------------------------------------------ ! if (present(flag_string)) then @@ -1444,13 +1444,11 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & end if ! if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (errmsg, *) ' Node number in list is outside of the grid', nodeu + write (errmsg, '(a,i0,a)') & + "Node number in list (", nodeu, ") is outside of the grid. "// & + "Cell number cannot be determined in line '"// & + trim(adjustl(line))//"'." call store_error(errmsg) - inquire (unit=in, name=fname) - call store_error('Error converting in file: ') - call store_error(trim(adjustl(fname))) - call store_error('Cell number cannot be determined in line: ') - call store_error(trim(adjustl(line))) call store_error_unit(in) end if ! @@ -1487,7 +1485,6 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & integer(I4B) :: lloclocal, istart, istop, ndum, n integer(I4B) :: istat real(DP) :: r - character(len=LINELENGTH) :: fname ! ------------------------------------------------------------------------------ ! if (present(flag_string)) then @@ -1516,13 +1513,11 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & end if ! if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (errmsg, *) ' Node number in list is outside of the grid', nodeu + write (errmsg, '(a,i0,a)') & + "Cell number cannot be determined for cellid ("// & + trim(adjustl(cellid))//") and results in a user "// & + "node number (", nodeu, ") that is outside of the grid." call store_error(errmsg) - inquire (unit=inunit, name=fname) - call store_error('Error converting in file: ') - call store_error(trim(adjustl(fname))) - call store_error('Cell number cannot be determined in cellid: ') - call store_error(trim(adjustl(cellid))) call store_error_unit(inunit) end if ! diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90 index 63630cc0190..17895b55f38 100644 --- a/src/Model/GroundWaterFlow/gwf3disv8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disv8.f90 @@ -9,6 +9,7 @@ module GwfDisvModule ubdsv06 use SimModule, only: count_errors, store_error, store_error_unit, & store_error_filename + use SimVariablesModule, only: errmsg use DisvGeom, only: DisvGeomType use MemoryManagerModule, only: mem_allocate use TdisModule, only: kstp, kper, pertim, totim, delt @@ -444,7 +445,6 @@ subroutine grid_finalize(this) integer(I4B) :: node, noder, j, k real(DP) :: top real(DP) :: dz - character(len=300) :: ermsg ! -- formats character(len=*), parameter :: fmtdz = & "('CELL (',i0,',',i0,') THICKNESS <= 0. ', & @@ -484,8 +484,8 @@ subroutine grid_finalize(this) end if dz = top - this%bot2d(j, k) if (dz <= DZERO) then - write (ermsg, fmt=fmtdz) k, j, top, this%bot2d(j, k) - call store_error(ermsg) + write (errmsg, fmt=fmtdz) k, j, top, this%bot2d(j, k) + call store_error(errmsg) end if end if end do @@ -724,7 +724,6 @@ subroutine connect(this) integer(I4B) :: narea_eq_zero integer(I4B) :: narea_lt_zero real(DP) :: area - character(len=LINELENGTH) :: errmsg ! ------------------------------------------------------------------------------ ! ! -- Initialize @@ -740,14 +739,14 @@ subroutine connect(this) end do if (area < DZERO) then narea_lt_zero = narea_lt_zero + 1 - write (errmsg, '(a,i0)') & - &'Calculated CELL2D area less than zero for cell ', j + write (errmsg, '(a,i0,a)') & + &'Calculated CELL2D area less than zero for cell ', j, '.' call store_error(errmsg) end if if (area == DZERO) then narea_eq_zero = narea_eq_zero + 1 - write (errmsg, '(a,i0)') & - 'Calculated CELL2D area is zero for cell ', j + write (errmsg, '(a,i0,a)') & + 'Calculated CELL2D area is zero for cell ', j, '.' call store_error(errmsg) end if end do @@ -755,16 +754,16 @@ subroutine connect(this) ! -- check for errors if (count_errors() > 0) then if (narea_lt_zero > 0) then - write (errmsg, '(i0, a)') narea_lt_zero, & - ' cell(s) have an area less than zero. Calculated cell & - &areas must be greater than zero. Negative areas often & + write (errmsg, '(i0,a)') narea_lt_zero, & + ' cell(s) have an area less than zero. Calculated cell & + &areas must be greater than zero. Negative areas often & &mean vertices are not listed in clockwise order.' call store_error(errmsg) end if if (narea_eq_zero > 0) then - write (errmsg, '(i0, a)') narea_eq_zero, & - ' cell(s) have an area equal to zero. Calculated cell & - &areas must be greater than zero. Calculated cell & + write (errmsg, '(i0,a)') narea_eq_zero, & + ' cell(s) have an area equal to zero. Calculated cell & + &areas must be greater than zero. Calculated cell & &areas equal to zero indicate that the cell is not defined & &by a valid polygon.' call store_error(errmsg) @@ -974,7 +973,6 @@ subroutine nodeu_to_array(this, nodeu, arr) integer(I4B), intent(in) :: nodeu integer(I4B), dimension(:), intent(inout) :: arr ! -- local - character(len=LINELENGTH) :: errmsg integer(I4B) :: isize integer(I4B) :: i, j, k ! ------------------------------------------------------------------------------ @@ -984,7 +982,7 @@ subroutine nodeu_to_array(this, nodeu, arr) if (isize /= this%ndim) then write (errmsg, '(a,i0,a,i0,a)') & 'Program error: nodeu_to_array size of array (', isize, & - ') is not equal to the discretization dimension (', this%ndim, ')' + ') is not equal to the discretization dimension (', this%ndim, ').' call store_error(errmsg, terminate=.TRUE.) end if ! @@ -1015,7 +1013,6 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) integer(I4B), intent(in) :: nodeu integer(I4B), intent(in) :: icheck ! -- local - character(len=LINELENGTH) :: errmsg ! ------------------------------------------------------------------------------ ! ! -- check the node number if requested @@ -1023,10 +1020,11 @@ function get_nodenumber_idx1(this, nodeu, icheck) result(nodenumber) ! ! -- If within valid range, convert to reduced nodenumber if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (errmsg, '(a,i10)') & - 'Nodenumber less than 1 or greater than nodes:', nodeu - call store_error(errmsg) nodenumber = 0 + write (errmsg, '(a,i0,a,i0,a)') & + 'Node number (', nodeu, ') is less than 1 or greater than nodes (', & + this%nodesuser, ').' + call store_error(errmsg) else nodenumber = nodeu if (this%nodes < this%nodesuser) nodenumber = this%nodereduced(nodeu) @@ -1058,7 +1056,6 @@ function get_nodenumber_idx2(this, k, j, icheck) & integer(I4B), intent(in) :: k, j integer(I4B), intent(in) :: icheck ! -- local - character(len=LINELENGTH) :: errmsg integer(I4B) :: nodeu ! formats character(len=*), parameter :: fmterr = & @@ -1076,17 +1073,30 @@ function get_nodenumber_idx2(this, k, j, icheck) & ! -- check the node number if requested if (icheck /= 0) then ! - if (k < 1 .or. k > this%nlay) & - call store_error('Layer less than one or greater than nlay') - if (j < 1 .or. j > this%ncpl) & - call store_error('Node number less than one or greater than ncpl') + errmsg = "" + ! + if (k < 1 .or. k > this%nlay) then + write (errmsg, '(a,i0,a)') & + 'Layer number in list (', k, ') is outside of the grid.' + end if + if (j < 1 .or. j > this%ncpl) then + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Node number in list (', j, & + ') is outside of the grid.' + end if ! ! -- Error if outside of range if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (errmsg, '(a,i10)') & - 'Nodenumber less than 1 or greater than nodes:', nodeu + write (errmsg, '(a,1x,a,i0,a,i0,a)') & + trim(adjustl(errmsg)), & + 'Node number (', nodeu, ') is less than 1 or greater '// & + 'than nodes (', this%nodesuser, ').' + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then call store_error(errmsg) end if + ! end if ! ! -- return @@ -1378,7 +1388,6 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & integer(I4B) :: j, k, nlay, nrow, ncpl integer(I4B) :: lloclocal, ndum, istat, n real(DP) :: r - character(len=LINELENGTH) :: ermsg, fname ! ------------------------------------------------------------------------------ ! if (present(flag_string)) then @@ -1411,24 +1420,30 @@ function nodeu_from_string(this, lloc, istart, istop, in, iout, line, & end if end if ! + errmsg = '' + ! if (k < 1 .or. k > nlay) then - write (ermsg, *) ' Layer number in list is outside of the grid', k - call store_error(ermsg) + write (errmsg, '(a,i0,a)') & + 'Layer number in list (', k, ') is outside of the grid.' end if if (j < 1 .or. j > ncpl) then - write (ermsg, *) ' Cell2d number in list is outside of the grid', j - call store_error(ermsg) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Cell2d number in list (', j, & + ') is outside of the grid.' end if + ! nodeu = get_node(k, 1, j, nlay, nrow, ncpl) ! if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (ermsg, *) ' Node number in list is outside of the grid', nodeu - call store_error(ermsg) - inquire (unit=in, name=fname) - call store_error('Error converting in file: ') - call store_error(trim(adjustl(fname))) - call store_error('Cell number cannot be determined in line: ') - call store_error(trim(adjustl(line))) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), & + "Node number in list (", nodeu, ") is outside of the grid. "// & + "Cell number cannot be determined in line '"// & + trim(adjustl(line))//"'." + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then + call store_error(errmsg) call store_error_unit(in) end if ! @@ -1465,7 +1480,6 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & integer(I4B) :: lloclocal, ndum, istat, n integer(I4B) :: istart, istop real(DP) :: r - character(len=LINELENGTH) :: ermsg, fname ! ------------------------------------------------------------------------------ ! if (present(flag_string)) then @@ -1499,24 +1513,30 @@ function nodeu_from_cellid(this, cellid, inunit, iout, flag_string, & end if end if ! + errmsg = '' + ! if (k < 1 .or. k > nlay) then - write (ermsg, *) ' Layer number in list is outside of the grid', k - call store_error(ermsg) + write (errmsg, '(a,i0,a)') & + 'Layer number in list (', k, ') is outside of the grid.' end if if (j < 1 .or. j > ncpl) then - write (ermsg, *) ' Cell2d number in list is outside of the grid', j - call store_error(ermsg) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), 'Cell2d number in list (', j, & + ') is outside of the grid.' end if + ! nodeu = get_node(k, 1, j, nlay, nrow, ncpl) ! if (nodeu < 1 .or. nodeu > this%nodesuser) then - write (ermsg, *) ' Node number in list is outside of the grid', nodeu - call store_error(ermsg) - inquire (unit=inunit, name=fname) - call store_error('Error converting in file: ') - call store_error(trim(adjustl(fname))) - call store_error('Cell number cannot be determined in cellid: ') - call store_error(trim(adjustl(cellid))) + write (errmsg, '(a,1x,a,i0,a)') & + trim(adjustl(errmsg)), & + "Cell number cannot be determined for cellid ("// & + trim(adjustl(cellid))//") and results in a user "// & + "node number (", nodeu, ") that is outside of the grid." + end if + ! + if (len_trim(adjustl(errmsg)) > 0) then + call store_error(errmsg) call store_error_unit(inunit) end if ! @@ -1909,7 +1929,6 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & integer(I4B), intent(in) :: iout ! -- local integer(I4B) :: il, ir, ic, ncol, nrow, nlay, nval, nodeu, noder, ipos, ierr - character(len=LINELENGTH) :: errmsg ! ------------------------------------------------------------------------------ ! ! -- set variables @@ -1928,7 +1947,8 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & nodeu = get_node(1, ir, ic, nlay, nrow, ncol) il = this%ibuff(nodeu) if (il < 1 .or. il > nlay) then - write (errmsg, *) 'Invalid layer number: ', il + write (errmsg, '(a,i0,a)') & + 'Invalid layer number (', il, ').' call store_error(errmsg, terminate=.TRUE.) end if nodeu = get_node(il, ir, ic, nlay, nrow, ncol) @@ -1945,8 +1965,8 @@ subroutine nlarray_to_nodelist(this, nodelist, maxbnd, nbound, aname, & ! -- Check for errors nbound = ipos - 1 if (ierr > 0) then - write (errmsg, '(a, i0)') & - 'MAXBOUND dimension is too small. Increase MAXBOUND to ', ierr + write (errmsg, '(a,i0,a)') & + 'MAXBOUND dimension is too small. Increase MAXBOUND to ', ierr, '.' call store_error(errmsg, terminate=.TRUE.) end if ! diff --git a/src/Model/GroundWaterFlow/gwf3sfr8.f90 b/src/Model/GroundWaterFlow/gwf3sfr8.f90 index 2a043350ec7..a7fcee06672 100644 --- a/src/Model/GroundWaterFlow/gwf3sfr8.f90 +++ b/src/Model/GroundWaterFlow/gwf3sfr8.f90 @@ -28,7 +28,7 @@ module SfrModule use BudgetObjectModule, only: BudgetObjectType, budgetobject_cr use TableModule, only: TableType, table_cr use ObserveModule, only: ObserveType - use InputOutputModule, only: extract_idnum_or_bndname + use InputOutputModule, only: extract_idnum_or_bndname, upcase use BaseDisModule, only: DisBaseType use SimModule, only: count_errors, store_error, store_error_unit, & store_warning, deprecation_warning @@ -844,7 +844,6 @@ subroutine sfr_read_packagedata(this) ! -- local variables character(len=LINELENGTH) :: text character(len=LINELENGTH) :: cellid - character(len=LINELENGTH) :: keyword character(len=10) :: cnum character(len=LENBOUNDNAME) :: bndName character(len=LENBOUNDNAME) :: bndNameTemp @@ -904,19 +903,30 @@ subroutine sfr_read_packagedata(this) call this%parser%GetCellid(this%dis%ndim, cellid, flag_string=.true.) this%igwfnode(n) = this%dis%noder_from_cellid(cellid, this%inunit, & this%iout, & - flag_string=.true.) + flag_string=.true., & + allow_zero=.true.) this%igwftopnode(n) = this%igwfnode(n) ! ! -- read the cellid string and determine if 'none' is specified if (this%igwfnode(n) < 1) then - call this%parser%GetStringCaps(keyword) this%ianynone = this%ianynone + 1 - if (keyword /= 'NONE') then + call upcase(cellid) + if (cellid == 'NONE') then + call this%parser%GetStringCaps(cellid) + ! + ! -- create warning message write (cnum, '(i0)') n - errmsg = 'Cell ID ('//trim(cellid)// & - ') for unconnected reach '//trim(cnum)// & - ' must be NONE' - call store_error(errmsg) + warnmsg = 'CELLID for unconnected reach '//trim(cnum)// & + ' specified to be NONE. Unconnected reaches '// & + 'should be specified with a zero for each grid '// & + 'dimension. For example, for a DIS grid a CELLID '// & + 'of 0 0 0 should be specified for unconnected reaches' + ! + ! -- create deprecation warning + call deprecation_warning('PACKAGEDATA', 'CELLID=NONE', '6.5.0', & + warnmsg, this%parser%GetUnit()) + else + end if end if ! -- get reach length From 8b9423855a2ae03e484fee353643c54fe4685e58 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Fri, 28 Jul 2023 15:18:26 -0400 Subject: [PATCH 06/13] ci: test with gfortran 13 (#1298) * consolidate gfortran test matrix * introduce gfortran 13 on linux and mac * limit autotests to latest on each platform --- .github/workflows/ci.yml | 130 +++++++++++++-------------------------- 1 file changed, 42 insertions(+), 88 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 4a0eef2e97b..041ba51cb2b 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -17,7 +17,7 @@ on: - 'doc/**' jobs: lint: - name: Lint (fprettify) + name: Check format runs-on: ubuntu-latest defaults: run: @@ -38,14 +38,14 @@ jobs: run: python .github/common/fortran_format_check.py build: - name: Build (gfortran 12) + name: Build runs-on: ubuntu-22.04 defaults: run: shell: bash -l {0} env: FC: gfortran - GCC_V: 12 + GCC_V: 13 steps: - name: Checkout modflow6 @@ -74,14 +74,14 @@ jobs: run: meson test --verbose --no-rebuild -C builddir smoke_test: - name: Smoke test (gfortran 12) + name: Smoke test runs-on: ubuntu-22.04 defaults: run: shell: bash -l {0} env: FC: gfortran - GCC_V: 12 + GCC_V: 13 steps: - name: Checkout modflow6 uses: actions/checkout@v3 @@ -127,8 +127,8 @@ jobs: pytest -v -n auto --durations 0 -S fi - test_gfortran_latest: - name: Test (gfortran 12) + test_gfortran: + name: Test (gfortran) needs: - lint - build @@ -137,13 +137,33 @@ jobs: strategy: fail-fast: false matrix: - os: [ ubuntu-22.04, macos-12, windows-2022 ] + # compatible combinations from https://github.com/awvwgk/setup-fortran#runner-compatibility + include: + - {os: ubuntu-20.04, gcc_v: 7, test: false} + - {os: ubuntu-20.04, gcc_v: 8, test: false} + - {os: ubuntu-20.04, gcc_v: 9, test: false} + - {os: ubuntu-20.04, gcc_v: 10, test: false} + - {os: ubuntu-20.04, gcc_v: 11, test: false} + - {os: ubuntu-22.04, gcc_v: 12, test: false} + # only run autotests on latest version on each platform + - {os: ubuntu-22.04, gcc_v: 13, test: true} + - {os: macos-12, gcc_v: 7, test: false} + - {os: macos-12, gcc_v: 8, test: false} + - {os: macos-12, gcc_v: 9, test: false} + - {os: macos-12, gcc_v: 10, test: false} + - {os: macos-12, gcc_v: 11, test: false} + - {os: macos-12, gcc_v: 12, test: false} + - {os: macos-12, gcc_v: 13, test: true} + - {os: windows-2022, gcc_v: 9, test: false} + - {os: windows-2022, gcc_v: 10, test: false} + - {os: windows-2022, gcc_v: 11, test: false} + - {os: windows-2022, gcc_v: 12, test: true} + defaults: run: shell: bash -l {0} env: FC: gfortran - GCC_V: 12 steps: - name: Checkout modflow6 uses: actions/checkout@v3 @@ -162,11 +182,11 @@ jobs: repository: MODFLOW-USGS/modflow6-examples path: modflow6-examples - - name: Setup GNU Fortran ${{ env.GCC_V }} + - name: Setup GNU Fortran ${{ matrix.gcc_v }} uses: awvwgk/setup-fortran@main with: compiler: gcc - version: ${{ env.GCC_V }} + version: ${{ matrix.gcc_v }} - name: Setup Micromamba uses: mamba-org/setup-micromamba@v1 @@ -183,13 +203,18 @@ jobs: run: | meson setup builddir -Ddebug=false --prefix=$(pwd) --libdir=bin meson install -C builddir - meson test --verbose --no-rebuild -C builddir + + - name: Unit test modflow6 + working-directory: modflow6 + run: meson test --verbose --no-rebuild -C builddir - name: Update flopy + if: matrix.test working-directory: modflow6/autotest run: python update_flopy.py - name: Get executables + if: matrix.test working-directory: modflow6/autotest env: GITHUB_TOKEN: ${{ github.token }} @@ -197,6 +222,7 @@ jobs: pytest -v --durations 0 get_exes.py - name: Test modflow6 + if: matrix.test working-directory: modflow6/autotest env: REPOS_PATH: ${{ github.workspace }} @@ -207,17 +233,15 @@ jobs: pytest -v -n auto --durations 0 -m "not large" fi - # steps below run only on Linux to test distribution procedures, e.g. - # compiling binaries, building documentation - name: Checkout usgslatex - if: runner.os == 'Linux' + if: matrix.test && runner.os == 'Linux' uses: actions/checkout@v3 with: repository: MODFLOW-USGS/usgslatex path: usgslatex - name: Install TeX Live - if: runner.os == 'Linux' + if: matrix.test && runner.os == 'Linux' run: | sudo apt-get update sudo apt install texlive-science \ @@ -227,87 +251,17 @@ jobs: texlive-fonts-extra - name: Install USGS LaTeX style files and Univers font - if: runner.os == 'Linux' + if: matrix.test && runner.os == 'Linux' working-directory: usgslatex/usgsLaTeX run: sudo ./install.sh --all-users - name: Test distribution scripts + if: matrix.test working-directory: modflow6/distribution env: GITHUB_TOKEN: ${{ github.token }} run: pytest -v --durations 0 - test_gfortran_previous: - name: Test gfortran (${{ matrix.GCC_V }}, ${{ matrix.os }}) - needs: - - lint - - build - - smoke_test - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - os: [ ubuntu-20.04 ] - GCC_V: [ 7, 8, 9, 10, 11 ] - defaults: - run: - shell: bash -l {0} - env: - FC: gfortran - steps: - - - name: Checkout modflow6 - uses: actions/checkout@v3 - with: - path: modflow6 - - - name: Checkout modflow6-testmodels - uses: actions/checkout@v3 - with: - repository: MODFLOW-USGS/modflow6-testmodels - path: modflow6-testmodels - - - name: Setup GNU Fortran ${{ matrix.GCC_V }} - uses: awvwgk/setup-fortran@main - with: - compiler: gcc - version: ${{ matrix.GCC_V }} - - - name: Setup Micromamba - uses: mamba-org/setup-micromamba@v1 - with: - environment-file: modflow6/environment.yml - cache-downloads: true - cache-environment: true - - - name: Update flopy - working-directory: modflow6/autotest - run: python update_flopy.py - - - name: Build modflow6 - working-directory: modflow6 - run: | - meson setup builddir -Ddebug=false --prefix=$(pwd) --libdir=bin - meson install -C builddir - meson test --verbose --no-rebuild -C builddir - - - name: Get executables - working-directory: modflow6/autotest - env: - GITHUB_TOKEN: ${{ github.token }} - run: pytest -v --durations 0 get_exes.py - - - name: Test modflow6 - working-directory: modflow6/autotest - env: - REPOS_PATH: ${{ github.workspace }} - run: | - if [ "${{ github.ref_name }}" == "master" ]; then - pytest -v -n auto --durations 0 -m "not large and not developmode" - else - pytest -v -n auto --durations 0 -m "not large" - fi - test_ifort: name: Test (ifort) needs: From 13d9b59ee0bd864ba3e05e17def1181bf68c4d64 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Fri, 28 Jul 2023 21:04:21 -0400 Subject: [PATCH 07/13] ci: expand Intel testing, use awvwgk/setup-fortran (#1303) --- .github/workflows/ci.yml | 124 ++++++++++++++++++++++++--------------- autotest/get_exes.py | 24 ++++++-- 2 files changed, 95 insertions(+), 53 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 041ba51cb2b..fa64ce99364 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -70,6 +70,10 @@ jobs: - name: Meson compile run: meson compile -C builddir + - name: Show build log + if: failure() + run: cat builddir/meson-logs/meson-log.txt + - name: Meson test run: meson test --verbose --no-rebuild -C builddir @@ -106,7 +110,15 @@ jobs: run: | meson setup builddir -Ddebug=false --prefix=$(pwd) --libdir=bin meson install -C builddir - meson test --verbose --no-rebuild -C builddir + + - name: Show build log + if: failure() + working-directory: modflow6 + run: cat builddir/meson-logs/meson-log.txt + + - name: Unit test programs + working-directory: modflow6 + run: meson test --verbose --no-rebuild -C builddir - name: Update flopy working-directory: modflow6/autotest @@ -204,7 +216,12 @@ jobs: meson setup builddir -Ddebug=false --prefix=$(pwd) --libdir=bin meson install -C builddir - - name: Unit test modflow6 + - name: Show build log + if: failure() + working-directory: modflow6 + run: cat builddir/meson-logs/meson-log.txt + + - name: Unit test programs working-directory: modflow6 run: meson test --verbose --no-rebuild -C builddir @@ -262,8 +279,8 @@ jobs: GITHUB_TOKEN: ${{ github.token }} run: pytest -v --durations 0 - test_ifort: - name: Test (ifort) + test_intel_fortran: + name: Test (Intel) needs: - lint - build @@ -272,7 +289,47 @@ jobs: strategy: fail-fast: false matrix: - os: [ ubuntu-latest, macos-latest, windows-latest ] + include: + ### ifx + ## 2022.2.x autotests disabled + # - mf5to6 test_evt: failure to converge + # - mf6 Keating_[disu_]dev: bad head comparison + - {os: ubuntu-22.04, compiler: intel, version: 2022.2.1, test: false} + - {os: ubuntu-22.04, compiler: intel, version: 2022.2, test: false} + ## 2021.1 segfault in meson serial sim test + # - {os: ubuntu-22.04, compiler: intel, version: 2022.1, test: false} + ## 2022.0 & 2021.[1,2,4] segfault at compile time + # - {os: ubuntu-22.04, compiler: intel, version: "2022.0", test: false} + # - {os: ubuntu-22.04, compiler: intel, version: 2021.4, test: false} + # - {os: ubuntu-22.04, compiler: intel, version: 2021.2, test: false} + # - {os: ubuntu-22.04, compiler: intel, version: 2021.1, test: false} + ## ifx not yet supported on macOS + # - {os: macos-12, compiler: intel, version: 2023.2, test: true} + ## 2023.[0,1] fail to compile + # - {os: windows-2022, compiler: intel, version: 2023.1, test: false} + # - {os: windows-2022, compiler: intel, version: "2023.0", test: false} + - {os: windows-2022, compiler: intel, version: 2022.2, test: false} + ## 2022.1 fail to link + # - {os: windows-2022, compiler: intel, version: 2022.1, test: false} + + ### ifort + ## only autotest latest on each platform + - {os: ubuntu-22.04, compiler: intel-classic, version: "2021.10", test: false} + - {os: ubuntu-22.04, compiler: intel-classic, version: 2021.9, test: false} + - {os: ubuntu-22.04, compiler: intel-classic, version: 2021.8, test: false} + - {os: ubuntu-22.04, compiler: intel-classic, version: 2021.7, test: true} + - {os: ubuntu-22.04, compiler: intel-classic, version: 2021.6, test: false} + - {os: macos-12, compiler: intel-classic, version: "2021.10", test: false} + - {os: macos-12, compiler: intel-classic, version: 2021.9, test: false} + - {os: macos-12, compiler: intel-classic, version: 2021.8, test: false} + - {os: macos-12, compiler: intel-classic, version: 2021.7, test: true} + - {os: macos-12, compiler: intel-classic, version: 2021.6, test: false} + - {os: windows-2022, compiler: intel-classic, version: "2021.10", test: false} + - {os: windows-2022, compiler: intel-classic, version: 2021.9, test: false} + - {os: windows-2022, compiler: intel-classic, version: 2021.8, test: false} + - {os: windows-2022, compiler: intel-classic, version: 2021.7, test: true} + - {os: windows-2022, compiler: intel-classic, version: 2021.6, test: false} + defaults: run: shell: bash -l {0} @@ -300,50 +357,44 @@ jobs: cache-downloads: true - name: Setup Intel Fortran - uses: modflowpy/install-intelfortran-action@v1 + uses: awvwgk/setup-fortran@main + with: + compiler: ${{ matrix.compiler }} + version: ${{ matrix.version }} - name: Update version files working-directory: modflow6/distribution run: python update_version.py - name: Build modflow6 - if: runner.os != 'Windows' working-directory: modflow6 run: | meson setup builddir -Ddebug=false --prefix=$(pwd) --libdir=bin meson install -C builddir - meson test --verbose --no-rebuild -C builddir - - name: Build modflow6 (Windows) - if: runner.os == 'Windows' + - name: Show build log + if: failure() working-directory: modflow6 - shell: pwsh - run: | - meson setup builddir -Ddebug=false --prefix=$(pwd) --libdir=bin - meson install -C builddir - meson test --verbose --no-rebuild -C builddir + run: cat builddir/meson-logs/meson-log.txt + + - name: Unit test programs + working-directory: modflow6 + run: meson test --verbose --no-rebuild -C builddir - name: Update flopy + if: matrix.test working-directory: modflow6/autotest run: python update_flopy.py - name: Get executables - if: runner.os != 'Windows' - working-directory: modflow6/autotest - env: - GITHUB_TOKEN: ${{ github.token }} - run: pytest -v --durations 0 get_exes.py - - - name: Get executables (Windows) - if: runner.os == 'Windows' + if: matrix.test working-directory: modflow6/autotest - shell: pwsh env: GITHUB_TOKEN: ${{ github.token }} run: pytest -v --durations 0 get_exes.py - name: Test programs - if: runner.os != 'Windows' + if: matrix.test working-directory: modflow6/autotest env: REPOS_PATH: ${{ github.workspace }} @@ -354,30 +405,9 @@ jobs: pytest -v -n auto --durations 0 -m "not large" fi - - name: Test programs (Windows) - if: runner.os == 'Windows' - working-directory: modflow6/autotest - shell: pwsh - env: - REPOS_PATH: ${{ github.workspace }} - run: | - if ( "${{ github.ref_name }}" -eq "master" ) { - pytest -v -n auto --durations 0 -m "not large and not developmode" - } else { - pytest -v -n auto --durations 0 -m "not large" - } - - name: Test scripts - if: runner.os != 'Windows' - working-directory: modflow6/distribution - env: - GITHUB_TOKEN: ${{ github.token }} - run: pytest -v --durations 0 - - - name: Test scripts (Windows) - if: runner.os == 'Windows' + if: matrix.test working-directory: modflow6/distribution - shell: pwsh env: GITHUB_TOKEN: ${{ github.token }} run: pytest -v --durations 0 diff --git a/autotest/get_exes.py b/autotest/get_exes.py index 0a071d50fea..f8c92ccfc85 100644 --- a/autotest/get_exes.py +++ b/autotest/get_exes.py @@ -1,7 +1,9 @@ import argparse +from os import environ from pathlib import Path from tempfile import TemporaryDirectory from warnings import warn +from platform import system import flopy import pytest @@ -9,7 +11,7 @@ from flaky import flaky from modflow_devtools.build import meson_build from modflow_devtools.download import download_and_unzip, get_release -from modflow_devtools.misc import get_ostag +from modflow_devtools.misc import get_ostag, is_in_ci, set_env repository = "MODFLOW-USGS/modflow6" top_bin_path = project_root_path / "bin" @@ -73,11 +75,21 @@ def test_rebuild_release(rebuilt_bin_path: Path): f.write(f"{line}\n") # rebuild with Meson - meson_build( - project_path=source_files_path.parent, - build_path=download_path / "builddir", - bin_path=rebuilt_bin_path, - ) + def rebuild(): + meson_build( + project_path=source_files_path.parent, + build_path=download_path / "builddir", + bin_path=rebuilt_bin_path, + ) + + # temp workaround until next release, + # ifx fails to build 6.4.2 on Windows + # most likely due to backspace issues + if system() == "Windows" and environ.get("FC") == "ifx": + with set_env(FC="ifort", CC="icl"): + rebuild() + else: + rebuild() @flaky(max_runs=3) From 5ad8543fc539810ea476df7a88263b223b54a511 Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Wed, 2 Aug 2023 11:11:53 -0400 Subject: [PATCH 08/13] docs(DEVELOPER.md): add intel compat table, describe branching model, update TOC (#1320) --- DEVELOPER.md | 50 +++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 11 deletions(-) diff --git a/DEVELOPER.md b/DEVELOPER.md index 1cfbba5524a..99e7b4d086d 100644 --- a/DEVELOPER.md +++ b/DEVELOPER.md @@ -7,6 +7,7 @@ To build and test a parallel version of the program, first read the instructions + - [Prerequisites](#prerequisites) - [Git](#git) - [Fortran compiler](#fortran-compiler) @@ -47,11 +48,13 @@ To build and test a parallel version of the program, first read the instructions - [Testing makefiles](#testing-makefiles) - [Installing `make` on Windows](#installing-make-on-windows) - [Using Conda from Git Bash](#using-conda-from-git-bash) -- [Git Strategy for Managing Long-Lived Branches](#git-strategy-for-managing-long-lived-branches) - - [Create a Backup](#create-a-backup) - - [Squash Feature Branch Commits](#squash-feature-branch-commits) - - [Rebase Feature Branch with Develop](#rebase-feature-branch-with-develop) - - [Cleanup](#cleanup) +- [Branching model](#branching-model) + - [Overview](#overview) + - [Managing long-lived branches](#managing-long-lived-branches) + - [Backup](#backup) + - [Squash](#squash) + - [Rebase](#rebase) + - [Cleanup](#cleanup) @@ -105,10 +108,27 @@ GNU Fortran can be installed on all three major platforms. #### Intel Fortran -Intel Fortran can also be used to compile MODFLOW 6 and associated utilities. The `ifort` compiler is available in the [Intel oneAPI HPC Toolkit](https://software.intel.com/content/www/us/en/develop/tools/oneapi/hpc-toolkit/download.html). An installer is bundled with the download. +Intel Fortran can also be used to compile MODFLOW 6 and associated utilities. The `ifort` and `ifx` compilers are available in the [Intel oneAPI HPC Toolkit](https://software.intel.com/content/www/us/en/develop/tools/oneapi/hpc-toolkit/download.html). A number of environment variables must be set before using Intel Fortran. General information can be found [here](https://www.intel.com/content/www/us/en/develop/documentation/oneapi-programming-guide/top/oneapi-development-environment-setup.html), with specific instructions to configure a shell session for `ifort` [here](https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/compiler-setup/use-the-command-line/specifying-the-location-of-compiler-components.html). +While the current development version of MODFLOW 6 is broadly compatible with `ifort`, `ifx` compatibility is still limited. The following table documents whether MODFLOW 6 can be succesfully built with particular platform/compiler combinations, and info on relevant errors if not. + +**Note**: this table is not exhaustive and only details the currently tested subset of combinations. + +| Platform | Compiler | Version | Compatible? | Notes | +|:---------|:---------|:--------|:------------|:------| +| Ubuntu 22.04 | ifort | 2021.[6-10] | ✓ | | +| Ubuntu 22.04 | ifx | 2022.2.[0-1] | ✓ | some autotests fail (convergence failure, bad head comparison) | +| Ubuntu 22.04 | ifx | 2022.1 | ✗ | segfault in meson serial simulation test | +| Ubuntu 22.04 | ifx | 2022.0, 2021.[1,2,4] | ✗ | compilation failure (segfault) | +| macOS 12 (Monterey) | ifort | 2021.[6-10] | ✓ | | +| macOS N | ifx | all | ✗ | ifx support for macOS [is not planned](https://community.intel.com/t5/Blogs/Tech-Innovation/Tools/Deprecation-of-Intel-Fortran-Compiler-Classic-for-macOS/post/1472697) | +| Windows 10 (Server 2022) | ifort | 2021.[6-10] | ✓ | | +| Windows 10 (Server 2022) | ifx | 2023.[0-1] | ✗ | compilation failure | +| Windows 10 (Server 2022) | ifx | 2022.2 | ✓ | | +| Windows 10 (Server 2022) | ifx | 2022.1 | ✗ | linking failure | + ##### Windows On Windows, [Visual Studio](https://visualstudio.microsoft.com) and a number of libraries must be installed for `ifort` to work. The required libraries can be installed by ticking the "Desktop Development with C++" checkbox in the Visual Studio Installer's Workloads tab. @@ -484,13 +504,21 @@ After this, `conda` commands should be available. This command may be added to a `.bashrc` or `.bash_profile` file in your home directory to permanently configure Git Bash for Conda. -## Git Strategy for Managing Long-Lived Branches +## Branching model + +This section documents MODFLOW 6 branching strategy and other VCS-related procedures. + +### Overview + +This project follows the [git flow](https://nvie.com/posts/a-successful-git-branching-model/): development occurs on the `develop` branch, while `master` is reserved for the state of the latest release. Development PRs are typically squashed to `develop` to avoid merge commits. At release time, release branches are merged to `master`, and then `master` is merged back into `develop`. + +### Managing long-lived branches When a feature branch takes a long time to develop, it is easy to become out of sync with the develop branch. Depending on the situation, it may be advisable to periodically squash the commits on the feature branch and rebase the change set with develop. The following approach for updating a long-lived feature branch has proven robust. In the example below, the feature branch is assumed to be called `feat-xyz`. -### Create a Backup +#### Backup Begin by creating a backup copy of the feature branch in case anything goes terribly wrong. @@ -500,7 +528,7 @@ git checkout -b feat-xyz-backup git checkout feat-xyz ``` -### Squash Feature Branch Commits +#### Squash Next, consider squashing commits on the feature branch. If there are many commits, it is beneficial to squash them before trying to rebase with develop. There is a nice article on [squashing commits into one using git](https://www.internalpointers.com/post/squash-commits-into-one-git), which has been very useful for consolidating commits on a long-lived modflow6 feature branch. @@ -519,7 +547,7 @@ git push origin feat-xyz --force The `--force` flag's short form is `-f`. -### Rebase Feature Branch with Develop +#### Rebase Now that the commits on `feat-xyz` have been consolidated, it is time to rebase with develop. If there are multiple commits in `feat-xyz` that make changes, undo them, rename files, and/or move things around in subsequent commits, then there may be multiple sets of merge conflicts that will need to be resolved as the rebase works its way through the commit change sets. This is why it is beneficial to squash the feature commits before rebasing with develop. @@ -539,7 +567,7 @@ At this point, you will want to force push the updated feature branch to origin git push origin feat-xyz --force ``` -### Cleanup +#### Cleanup Lastly, if you are satisfied with the results and confident the procedure went well, then you can delete the backup that you created at the start. From b2ee6cfb3a033ebf99293b07590d9dae5317498f Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Thu, 3 Aug 2023 12:25:43 -0500 Subject: [PATCH 09/13] refactor(LAK): specify no lake bed using DNODATA (#1321) - [x] deprecate use of NONE keyword to define lake-gwf conductance only a function of aquifer properties - [x] update LAKE tests that used NONE keyword to use DNODATA - [x] add new LAK test to evaluate use of DNODATA, NONE, and negative number (failure) for bedleak - [x] update mf6io - [x] update release notes - [x] update mwt and ssm tex files - a result of running `python mf6ivar.py` --- autotest/simulation.py | 1 + autotest/test_gwf_auxvars.py | 6 +- autotest/test_gwf_lak_bedleak.py | 193 +++++++++++++++++++++++++ autotest/test_gwt_lkt01.py | 8 +- autotest/test_gwt_lkt02.py | 12 +- autotest/test_gwt_lkt03.py | 8 +- autotest/test_gwt_lkt04.py | 8 +- doc/ReleaseNotes/develop.tex | 4 +- doc/mf6io/mf6ivar/dfn/gwf-lak.dfn | 2 +- doc/mf6io/mf6ivar/md/mf6ivar.md | 2 +- doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex | 2 +- doc/mf6io/mf6ivar/tex/gwt-mwt-desc.tex | 2 +- doc/mf6io/mf6ivar/tex/gwt-ssm-desc.tex | 2 +- src/Model/GroundWaterFlow/gwf3lak8.f90 | 80 ++++++---- 14 files changed, 273 insertions(+), 57 deletions(-) create mode 100644 autotest/test_gwf_lak_bedleak.py diff --git a/autotest/simulation.py b/autotest/simulation.py index 73516fbfdee..9ee16586398 100644 --- a/autotest/simulation.py +++ b/autotest/simulation.py @@ -17,6 +17,7 @@ from flopy.utils.compare import compare_heads from modflow_devtools.misc import is_in_ci +DNODATA = 3.e+30 sfmt = "{:25s} - {}" extdict = { "hds": "head", diff --git a/autotest/test_gwf_auxvars.py b/autotest/test_gwf_auxvars.py index 412b53a8385..95cd656527b 100644 --- a/autotest/test_gwf_auxvars.py +++ b/autotest/test_gwf_auxvars.py @@ -5,7 +5,7 @@ import numpy as np import pytest from framework import TestFramework -from simulation import TestSimulation +from simulation import TestSimulation, DNODATA ex = ["aux01"] auxvar1 = 101.0 @@ -183,8 +183,8 @@ def build_model(idx, dir): ] # connectiondata = [ - [0, 0, (0, 1, 1), "vertical", "none", 0.0, 0.0, 0.0, 0.0], - [1, 0, (0, 2, 2), "vertical", "none", 0.0, 0.0, 0.0, 0.0], + [0, 0, (0, 1, 1), "vertical", DNODATA, 0.0, 0.0, 0.0, 0.0], + [1, 0, (0, 2, 2), "vertical", DNODATA, 0.0, 0.0, 0.0, 0.0], ] lak = flopy.mf6.ModflowGwflak( gwf, diff --git a/autotest/test_gwf_lak_bedleak.py b/autotest/test_gwf_lak_bedleak.py new file mode 100644 index 00000000000..0787ff3cac1 --- /dev/null +++ b/autotest/test_gwf_lak_bedleak.py @@ -0,0 +1,193 @@ +import os +import sys + +import flopy +import numpy as np +import pytest +from framework import TestFramework +from simulation import TestSimulation, DNODATA + +ex = ["bedleak", "bedleak_fail", "bedleak_none"] + + +def build_model(idx, dir): + nlay, nrow, ncol = 1, 10, 10 + nper = 1 + perlen = [ + 1.0, + ] + nstp = [ + 1, + ] + tsmult = [ + 1.0, + ] + + lenx = 300.0 + delr = delc = lenx / float(nrow) + strt = 100.0 + + nouter, ninner = 100, 300 + hclose, rclose, relax = 1e-9, 1e-3, 0.97 + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = ex[idx] + + # build MODFLOW 6 files + ws = dir + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name="mf6", sim_ws=ws + ) + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, + time_units="DAYS", + nper=nper, + perioddata=tdis_rc, + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="DBD", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration="BICGSTAB", + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relax, + ) + + # create gwf model + gwf = flopy.mf6.ModflowGwf(sim, modelname=name) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=90.0, + botm=0.0, + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, save_flows=True, icelltype=1, k=1.0, k33=0.01 + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.0, + sy=0.1, + steady_state={0: True}, + ) + + # chd files + chdlist0 = [] + chdlist0.append([(0, 0, 0), 100.0]) + chdlist0.append([(0, nrow - 1, ncol - 1), 95.0]) + + chdspdict = {0: chdlist0} + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspdict, + save_flows=False, + ) + + # lak package + if "fail" in name: + bedleak = -100.0 + elif "none" in name: + bedleak = "none" + else: + bedleak = DNODATA + + # [] [] + packagedata = [ + [0, 100.0, 1, "lake1"], + [1, 100.0, 1, "lake2"], + ] + # + connectiondata = [ + [0, 0, (0, 1, 1), "vertical", bedleak, 0.0, 0.0, 0.0, 0.0], + [1, 0, (0, 2, 2), "vertical", bedleak, 0.0, 0.0, 0.0, 0.0], + ] + lak = flopy.mf6.ModflowGwflak( + gwf, + boundnames=True, + surfdep=1.0, + print_input=True, + print_stage=True, + print_flows=True, + save_flows=True, + budget_filerecord=f"{name}.lak.bud", + nlakes=len(packagedata), + packagedata=packagedata, + connectiondata=connectiondata, + ) + # lak.remove() + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + + return sim, None + + +def eval_model(sim): + print("evaluating model...") + + name = ex[sim.idxsim] + + # lak budget + if "fail" not in name: + fpth = os.path.join(sim.simpath, f"{name}.lak.bud") + bobj = flopy.utils.CellBudgetFile(fpth, precision="double") + bobj.list_unique_records() + records = bobj.get_data(text="GWF") + for r in records: + assert np.allclose(r["q"][0], -4.79616347e-12) + assert np.allclose(r["q"][1], -6.19237994e-12) + + +@pytest.mark.parametrize( + "idx, name", + list(enumerate(ex)), +) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + if "fail" in ws: + require_fail = True + else: + require_fail = False + test = TestFramework() + test.build(build_model, idx, ws) + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_model, + idxsim=idx, + require_failure=require_fail, + ), + str(function_tmpdir), + ) diff --git a/autotest/test_gwt_lkt01.py b/autotest/test_gwt_lkt01.py index 4971e5c9412..44ef5eb63c5 100644 --- a/autotest/test_gwt_lkt01.py +++ b/autotest/test_gwt_lkt01.py @@ -9,7 +9,7 @@ import numpy as np import pytest from framework import TestFramework -from simulation import TestSimulation +from simulation import TestSimulation, DNODATA ex = ["lkt_01"] @@ -134,13 +134,13 @@ def build_model(idx, dir): con_data = [] # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth ) con_data.append( - (0, 0, (0, 0, 1), "HORIZONTAL", "None", 10, 10, connlen, connwidth) + (0, 0, (0, 0, 1), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth) ) con_data.append( - (0, 1, (0, 0, 3), "HORIZONTAL", "None", 10, 10, connlen, connwidth) + (0, 1, (0, 0, 3), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth) ) con_data.append( - (0, 2, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth) + (0, 2, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) p_data = [ (0, "STATUS", "CONSTANT"), diff --git a/autotest/test_gwt_lkt02.py b/autotest/test_gwt_lkt02.py index 616fb514397..73b4daf6796 100644 --- a/autotest/test_gwt_lkt02.py +++ b/autotest/test_gwt_lkt02.py @@ -7,7 +7,7 @@ import numpy as np import pytest from framework import TestFramework -from simulation import TestSimulation +from simulation import TestSimulation, DNODATA ex = ["lkt_02"] @@ -136,21 +136,21 @@ def build_model(idx, dir): # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth ) # lake 1 con_data.append( - (0, 0, (0, 0, 1), "HORIZONTAL", "None", 10, 10, connlen, connwidth) + (0, 0, (0, 0, 1), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth) ) con_data.append( - (0, 1, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth) + (0, 1, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) # lake 2 con_data.append( - (1, 0, (0, 0, 3), "VERTICAL", "None", 10, 10, connlen, connwidth) + (1, 0, (0, 0, 3), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) # lake 3 con_data.append( - (2, 0, (0, 0, 4), "VERTICAL", "None", 10, 10, connlen, connwidth) + (2, 0, (0, 0, 4), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) con_data.append( - (2, 1, (0, 0, 5), "HORIZONTAL", "None", 10, 10, connlen, connwidth) + (2, 1, (0, 0, 5), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth) ) p_data = [ diff --git a/autotest/test_gwt_lkt03.py b/autotest/test_gwt_lkt03.py index efd98105085..0d3c55e99ae 100644 --- a/autotest/test_gwt_lkt03.py +++ b/autotest/test_gwt_lkt03.py @@ -7,7 +7,7 @@ import numpy as np import pytest from framework import TestFramework -from simulation import TestSimulation +from simulation import TestSimulation, DNODATA ex = ["lkt_03"] @@ -134,15 +134,15 @@ def build_model(idx, dir): # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth ) # lake 1 con_data.append( - (0, 0, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth) + (0, 0, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) # lake 2 con_data.append( - (1, 0, (0, 0, 3), "VERTICAL", "None", 10, 10, connlen, connwidth) + (1, 0, (0, 0, 3), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) # lake 3 con_data.append( - (2, 0, (0, 0, 4), "VERTICAL", "None", 10, 10, connlen, connwidth) + (2, 0, (0, 0, 4), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) p_data = [ diff --git a/autotest/test_gwt_lkt04.py b/autotest/test_gwt_lkt04.py index f8bb2cf2c65..888c3388042 100644 --- a/autotest/test_gwt_lkt04.py +++ b/autotest/test_gwt_lkt04.py @@ -10,7 +10,7 @@ import numpy as np import pytest from framework import TestFramework -from simulation import TestSimulation +from simulation import TestSimulation, DNODATA ex = ["lkt_04"] @@ -134,13 +134,13 @@ def build_model(idx, dir, exe): con_data = [] # con_data=(lakeno,iconn,(cellid),claktype,bedleak,belev,telev,connlen,connwidth ) con_data.append( - (0, 0, (0, 0, 1), "HORIZONTAL", "None", 10, 10, connlen, connwidth) + (0, 0, (0, 0, 1), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth) ) con_data.append( - (0, 1, (0, 0, 3), "HORIZONTAL", "None", 10, 10, connlen, connwidth) + (0, 1, (0, 0, 3), "HORIZONTAL", DNODATA, 10, 10, connlen, connwidth) ) con_data.append( - (0, 2, (0, 0, 2), "VERTICAL", "None", 10, 10, connlen, connwidth) + (0, 2, (0, 0, 2), "VERTICAL", DNODATA, 10, 10, connlen, connwidth) ) p_data = [ (0, "STATUS", "ACTIVE"), diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index b1fcec5c8c8..a7c72979d49 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -41,8 +41,8 @@ \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 unconnected reaches. Warning messages will be issued if NONE is specified for unconnected reaches. Specifying a CELLID of NONE will eventually be deprecated and will cause MODFLOW 6 to terminate with an error. - % \item xxx + \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 unconnected reaches. Warning messages will be issued if NONE is specified for unconnected reaches. 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 to identify lake-GWF connections where conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be 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 xxx \end{itemize} diff --git a/doc/mf6io/mf6ivar/dfn/gwf-lak.dfn b/doc/mf6io/mf6ivar/dfn/gwf-lak.dfn index 36a67c16dca..d229382bbc7 100644 --- a/doc/mf6io/mf6ivar/dfn/gwf-lak.dfn +++ b/doc/mf6io/mf6ivar/dfn/gwf-lak.dfn @@ -470,7 +470,7 @@ tagged false in_record true reader urword longname bed leakance -description character string or real value that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero or specified to be NONE. If BEDLEAK is specified to be NONE, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. +description real value or character string that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero, equal to the DNODATA value (3.0E+30), or specified to be NONE. If DNODATA or NONE is specified for BEDLEAK, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified. Eventually the ability to specify NONE will be deprecated and cause MODFLOW 6 to terminate with an error. block connectiondata name belev diff --git a/doc/mf6io/mf6ivar/md/mf6ivar.md b/doc/mf6io/mf6ivar/md/mf6ivar.md index 0bce538bf7a..9b45f871c9c 100644 --- a/doc/mf6io/mf6ivar/md/mf6ivar.md +++ b/doc/mf6io/mf6ivar/md/mf6ivar.md @@ -658,7 +658,7 @@ | GWF | LAK | CONNECTIONDATA | ICONN | INTEGER | integer value that defines the GWF connection number for this lake connection entry. ICONN must be greater than zero and less than or equal to NLAKECONN for lake LAKENO. | | GWF | LAK | CONNECTIONDATA | CELLID | INTEGER (NCELLDIM) | is the cell identifier, and depends on the type of grid that is used for the simulation. For a structured grid that uses the DIS input file, CELLID is the layer, row, and column. For a grid that uses the DISV input file, CELLID is the layer and CELL2D number. If the model uses the unstructured discretization (DISU) input file, CELLID is the node number for the cell. | | GWF | LAK | CONNECTIONDATA | CLAKTYPE | STRING | character string that defines the lake-GWF connection type for the lake connection. Possible lake-GWF connection type strings include: VERTICAL--character keyword to indicate the lake-GWF connection is vertical and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. HORIZONTAL--character keyword to indicate the lake-GWF connection is horizontal and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDH--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDV--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. Embedded lakes can only be connected to a single cell (NLAKECONN = 1) and there must be a lake table associated with each embedded lake. | -| GWF | LAK | CONNECTIONDATA | BEDLEAK | STRING | character string or real value that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero or specified to be NONE. If BEDLEAK is specified to be NONE, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. | +| GWF | LAK | CONNECTIONDATA | BEDLEAK | STRING | real value or character string that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero, equal to the DNODATA value (3.0E+30), or specified to be NONE. If DNODATA or NONE is specified for BEDLEAK, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified. Eventually the ability to specify NONE will be deprecated and cause MODFLOW 6 to terminate with an error. | | GWF | LAK | CONNECTIONDATA | BELEV | DOUBLE PRECISION | real value that defines the bottom elevation for a HORIZONTAL lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL, EMBEDDEDH, or EMBEDDEDV. If CLAKTYPE is HORIZONTAL and BELEV is not equal to TELEV, BELEV must be greater than or equal to the bottom of the GWF cell CELLID. If BELEV is equal to TELEV, BELEV is reset to the bottom of the GWF cell CELLID. | | GWF | LAK | CONNECTIONDATA | TELEV | DOUBLE PRECISION | real value that defines the top elevation for a HORIZONTAL lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL, EMBEDDEDH, or EMBEDDEDV. If CLAKTYPE is HORIZONTAL and TELEV is not equal to BELEV, TELEV must be less than or equal to the top of the GWF cell CELLID. If TELEV is equal to BELEV, TELEV is reset to the top of the GWF cell CELLID. | | GWF | LAK | CONNECTIONDATA | CONNLEN | DOUBLE PRECISION | real value that defines the distance between the connected GWF CELLID node and the lake for a HORIZONTAL, EMBEDDEDH, or EMBEDDEDV lake-GWF connection. CONLENN must be greater than zero for a HORIZONTAL, EMBEDDEDH, or EMBEDDEDV lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL. | diff --git a/doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex b/doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex index ad5ac4b7d71..c3886f5c8cb 100644 --- a/doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex +++ b/doc/mf6io/mf6ivar/tex/gwf-lak-desc.tex @@ -91,7 +91,7 @@ \item \texttt{claktype}---character string that defines the lake-GWF connection type for the lake connection. Possible lake-GWF connection type strings include: VERTICAL--character keyword to indicate the lake-GWF connection is vertical and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. HORIZONTAL--character keyword to indicate the lake-GWF connection is horizontal and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDH--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{11}$ tensor component defined for CELLID in the NPF package. EMBEDDEDV--character keyword to indicate the lake-GWF connection is embedded in a single cell and connection conductance calculations use the hydraulic conductivity corresponding to the $K_{33}$ tensor component defined for CELLID in the NPF package. Embedded lakes can only be connected to a single cell (NLAKECONN = 1) and there must be a lake table associated with each embedded lake. -\item \texttt{bedleak}---character string or real value that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero or specified to be NONE. If BEDLEAK is specified to be NONE, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. +\item \texttt{bedleak}---real value or character string that defines the bed leakance for the lake-GWF connection. BEDLEAK must be greater than or equal to zero, equal to the DNODATA value (3.0E+30), or specified to be NONE. If DNODATA or NONE is specified for BEDLEAK, the lake-GWF connection conductance is solely a function of aquifer properties in the connected GWF cell and lakebed sediments are assumed to be absent. Warning messages will be issued if NONE is specified. Eventually the ability to specify NONE will be deprecated and cause MODFLOW 6 to terminate with an error. \item \texttt{belev}---real value that defines the bottom elevation for a HORIZONTAL lake-GWF connection. Any value can be specified if CLAKTYPE is VERTICAL, EMBEDDEDH, or EMBEDDEDV. If CLAKTYPE is HORIZONTAL and BELEV is not equal to TELEV, BELEV must be greater than or equal to the bottom of the GWF cell CELLID. If BELEV is equal to TELEV, BELEV is reset to the bottom of the GWF cell CELLID. diff --git a/doc/mf6io/mf6ivar/tex/gwt-mwt-desc.tex b/doc/mf6io/mf6ivar/tex/gwt-mwt-desc.tex index b3bbbc87356..7c2c57c4e27 100644 --- a/doc/mf6io/mf6ivar/tex/gwt-mwt-desc.tex +++ b/doc/mf6io/mf6ivar/tex/gwt-mwt-desc.tex @@ -63,7 +63,7 @@ \item \texttt{mawno}---integer value that defines the well number associated with the specified PERIOD data on the line. MAWNO must be greater than zero and less than or equal to NMAWWELLS. -\item \texttt{mwtsetting}---line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWTSETTING string include: STATUS, CONCENTRATION, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the concentration of associated with the corresponding flow terms. Concentrations cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated concentration of the well. +\item \texttt{mwtsetting}---line of information that is parsed into a keyword and values. Keyword values that can be used to start the MWTSETTING string include: STATUS, CONCENTRATION, RAINFALL, EVAPORATION, RUNOFF, and AUXILIARY. These settings are used to assign the concentration associated with the corresponding flow terms. Concentrations cannot be specified for all flow terms. For example, the Multi-Aquifer Well Package supports a ``WITHDRAWAL'' flow term. If this withdrawal term is active, then water will be withdrawn from the well at the calculated concentration of the well. \begin{lstlisting}[style=blockdefinition] STATUS diff --git a/doc/mf6io/mf6ivar/tex/gwt-ssm-desc.tex b/doc/mf6io/mf6ivar/tex/gwt-ssm-desc.tex index 0d450f542f4..d75c8508def 100644 --- a/doc/mf6io/mf6ivar/tex/gwt-ssm-desc.tex +++ b/doc/mf6io/mf6ivar/tex/gwt-ssm-desc.tex @@ -13,7 +13,7 @@ \begin{description} \item \texttt{pname}---name of the flow package for which an auxiliary variable contains a source concentration. If this flow package is represented using an advanced transport package (SFT, LKT, MWT, or UZT), then the advanced transport package will override SSM terms specified here. -\item \texttt{srctype}---keyword indicating how concentration will be assigned for sources and sinks. Keyword must be specified as either AUX or AUXMIXED. For both options the user must provide an auxiliary variable in the corresponding flow package. The auxiliary variable must have the same name as the AUXNAME value that follows. If the AUX keyword is specified, then the auxiliary variable specified by the user will be assigned as the concenration value for groundwater sources (flows with a positive sign). For negative flow rates (sinks), groundwater will be withdrawn from the cell at the simulated concentration of the cell. The AUXMIXED option provides an alternative method for how to determine the concentration of sinks. If the cell concentration is larger than the user-specified auxiliary concentration, then the concentration of groundwater withdrawn from the cell will be assigned as the user-specified concentration. Alternatively, if the user-specified auxiliary concentration is larger than the cell concentration, then groundwater will be withdrawn at the cell concentration. Thus, the AUXMIXED option is designed to work with the Evapotranspiration (EVT) and Recharge (RCH) Packages where water may be withdrawn at a concentration that is less than the cell concentration. +\item \texttt{srctype}---keyword indicating how concentration will be assigned for sources and sinks. Keyword must be specified as either AUX or AUXMIXED. For both options the user must provide an auxiliary variable in the corresponding flow package. The auxiliary variable must have the same name as the AUXNAME value that follows. If the AUX keyword is specified, then the auxiliary variable specified by the user will be assigned as the concentration value for groundwater sources (flows with a positive sign). For negative flow rates (sinks), groundwater will be withdrawn from the cell at the simulated concentration of the cell. The AUXMIXED option provides an alternative method for how to determine the concentration of sinks. If the cell concentration is larger than the user-specified auxiliary concentration, then the concentration of groundwater withdrawn from the cell will be assigned as the user-specified concentration. Alternatively, if the user-specified auxiliary concentration is larger than the cell concentration, then groundwater will be withdrawn at the cell concentration. Thus, the AUXMIXED option is designed to work with the Evapotranspiration (EVT) and Recharge (RCH) Packages where water may be withdrawn at a concentration that is less than the cell concentration. \item \texttt{auxname}---name of the auxiliary variable in the package PNAME. This auxiliary variable must exist and be specified by the user in that package. The values in this auxiliary variable will be used to set the concentration associated with the flows for that boundary package. diff --git a/src/Model/GroundWaterFlow/gwf3lak8.f90 b/src/Model/GroundWaterFlow/gwf3lak8.f90 index 366ebc54e23..3c839ad81b0 100644 --- a/src/Model/GroundWaterFlow/gwf3lak8.f90 +++ b/src/Model/GroundWaterFlow/gwf3lak8.f90 @@ -1,6 +1,6 @@ module LakModule ! - use KindModule, only: DP, I4B + use KindModule, only: DP, I4B, LGP use ConstantsModule, only: LINELENGTH, LENBOUNDNAME, LENTIMESERIESNAME, & IWETLAKE, MAXADPIT, & DZERO, DPREC, DEM30, DEM9, DEM6, DEM5, & @@ -25,11 +25,12 @@ module LakModule use ObsModule, only: ObsType use InputOutputModule, only: get_node, URWORD, extract_idnum_or_bndname use BaseDisModule, only: DisBaseType - use SimModule, only: count_errors, store_error, store_error_unit + use SimModule, only: count_errors, store_error, store_error_unit, & + deprecation_warning use GenericUtilitiesModule, only: sim_message, is_same use BlockParserModule, only: BlockParserType use BaseDisModule, only: DisBaseType - use SimVariablesModule, only: errmsg + use SimVariablesModule, only: errmsg, warnmsg use MatrixBaseModule ! implicit none @@ -487,7 +488,7 @@ subroutine lak_read_lakes(this) character(len=9) :: cno character(len=50), dimension(:), allocatable :: caux integer(I4B) :: ierr, ival - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock integer(I4B) :: n integer(I4B) :: ii, jj integer(I4B) :: iaux @@ -723,13 +724,15 @@ subroutine lak_read_lake_connections(this) ! -- local character(len=LINELENGTH) :: keyword, cellid integer(I4B) :: ierr, ival - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock + logical(LGP) :: is_lake_bed real(DP) :: rval integer(I4B) :: j, n integer(I4B) :: nn integer(I4B) :: ipos, ipos0 integer(I4B) :: icellid, icellid0 - real(DP) :: top, bot + real(DP) :: top + real(DP) :: bot integer(I4B), dimension(:), pointer, contiguous :: nboundchk character(len=LENVARNAME) :: ctypenm @@ -838,16 +841,35 @@ subroutine lak_read_lake_connections(this) write (ctypenm, '(a16)') keyword ! -- bed leakance - !this%bedleak(ipos) = this%parser%GetDouble() + !this%bedleak(ipos) = this%parser%GetDouble() !TODO: use this when NONE keyword deprecated call this%parser%GetStringCaps(keyword) select case (keyword) case ('NONE') - this%bedleak(ipos) = -DONE + is_lake_bed = .FALSE. + this%bedleak(ipos) = DNODATA + ! + ! -- create warning message + write (warnmsg, '(2(a,1x,i0,1x),a,1pe7.1,a)') & + 'BEDLEAK for connection', j, 'in lake', n, 'is specified to '// & + 'be NONE. Lake connections where the lake-GWF connection '// & + 'conductance is solely a function of aquifer properties '// & + 'in the connected GWF cell should be specified with a '// & + 'DNODATA (', DNODATA, ') value.' + ! + ! -- create deprecation warning + call deprecation_warning('CONNECTIONDATA', 'bedleak=NONE', '6.5.0', & + warnmsg, this%parser%GetUnit()) case default - read (keyword, *) this%bedleak(ipos) + read (keyword, *) rval + if (is_same(rval, DNODATA)) then + is_lake_bed = .FALSE. + else + is_lake_bed = .TRUE. + end if + this%bedleak(ipos) = rval end select - if (keyword /= 'NONE' .and. this%bedleak(ipos) < dzero) then + if (is_lake_bed .and. this%bedleak(ipos) < DZERO) then write (errmsg, '(a,1x,i0,1x,a)') 'bedleak FOR LAKE ', n, 'MUST BE >= 0' call store_error(errmsg) end if @@ -1029,7 +1051,7 @@ subroutine lak_read_tables(this) character(len=LINELENGTH) :: line character(len=LINELENGTH) :: keyword integer(I4B) :: ierr - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock integer(I4B) :: n integer(I4B) :: iconn integer(I4B) :: ntabs @@ -1223,7 +1245,7 @@ subroutine lak_read_table(this, ilak, filename, laketable) ! -- local character(len=LINELENGTH) :: keyword integer(I4B) :: ierr - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock integer(I4B) :: iu integer(I4B) :: n integer(I4B) :: ipos @@ -1467,7 +1489,7 @@ subroutine lak_read_outlets(this) character(len=LENBOUNDNAME) :: bndName character(len=9) :: citem integer(I4B) :: ierr, ival - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock integer(I4B) :: n integer(I4B) :: jj integer(I4B), dimension(:), pointer, contiguous :: nboundchk @@ -1658,7 +1680,7 @@ subroutine lak_read_dimensions(this) ! -- local character(len=LINELENGTH) :: keyword integer(I4B) :: ierr - logical :: isfound, endOfBlock + logical(LGP) :: isfound, endOfBlock ! -- format ! ------------------------------------------------------------------------------ ! @@ -1882,10 +1904,10 @@ subroutine lak_read_initial_attr(this) end if length = this%connlength(j) end if - if (this%bedleak(j) < DZERO) then - clb(j) = -DONE + if (is_same(this%bedleak(j), DNODATA)) then + clb(j) = DNODATA else if (this%bedleak(j) > DZERO) then - clb(j) = done / this%bedleak(j) + clb(j) = DONE / this%bedleak(j) else clb(j) = DZERO end if @@ -1894,7 +1916,7 @@ subroutine lak_read_initial_attr(this) else caq(j) = DZERO end if - if (this%bedleak(j) < DZERO) then + if (is_same(this%bedleak(j), DNODATA)) then this%satcond(j) = area / caq(j) else if (clb(j) * caq(j) > DZERO) then this%satcond(j) = area / (clb(j) + caq(j)) @@ -1929,7 +1951,7 @@ subroutine lak_read_initial_attr(this) nn = this%cellid(j) area = this%warea(j) c1 = DZERO - if (clb(j) < DZERO) then + if (is_same(clb(j), DNODATA)) then cbedleak = ' NONE ' cbedcond = ' NONE ' else if (clb(j) > DZERO) then @@ -3421,7 +3443,7 @@ subroutine lak_options(this, option, found) ! -- dummy class(LakType), intent(inout) :: this character(len=*), intent(inout) :: option - logical, intent(inout) :: found + logical(LGP), intent(inout) :: found ! -- local character(len=MAXCHARLEN) :: fname, keyword real(DP) :: r @@ -3619,8 +3641,8 @@ subroutine lak_rp(this) character(len=LINELENGTH) :: title character(len=LINELENGTH) :: line character(len=LINELENGTH) :: text - logical :: isfound - logical :: endOfBlock + logical(LGP) :: isfound + logical(LGP) :: endOfBlock integer(I4B) :: ierr integer(I4B) :: node integer(I4B) :: n @@ -3830,12 +3852,12 @@ subroutine lak_cf(this, reset_mover) ! ------------------------------------------------------------------------------ ! -- dummy class(LakType) :: this - logical, intent(in), optional :: reset_mover + logical(LGP), intent(in), optional :: reset_mover ! -- local integer(I4B) :: j, n integer(I4B) :: igwfnode real(DP) :: hlak, blak - logical :: lrm + logical(LGP) :: lrm ! ------------------------------------------------------------------------------ !! !! -- Calculate lak conductance and update package RHS and HCOF @@ -5103,7 +5125,7 @@ subroutine lak_rp_obs(this) integer(I4B) :: nn2 integer(I4B) :: jj character(len=LENBOUNDNAME) :: bname - logical :: jfound + logical(LGP) :: jfound class(ObserveType), pointer :: obsrv => null() ! -------------------------------------------------------------------------- ! -- formats @@ -5420,11 +5442,11 @@ subroutine lak_solve(this, update) ! SPECIFICATIONS: ! -------------------------------------------------------------------------- use TdisModule, only: delt - logical, intent(in), optional :: update + logical(LGP), intent(in), optional :: update ! -- dummy class(LakType), intent(inout) :: this ! -- local - logical :: lupdate + logical(LGP) :: lupdate integer(I4B) :: i integer(I4B) :: j integer(I4B) :: n @@ -6639,8 +6661,8 @@ subroutine lak_calculate_density_exchange(this, iconn, stage, head, cond, & real(DP) :: elevavg real(DP) :: d1 real(DP) :: d2 - logical :: stage_below_bot - logical :: head_below_bot + logical(LGP) :: stage_below_bot + logical(LGP) :: head_below_bot ! -- formats ! ------------------------------------------------------------------------------ ! From 2f52a232a845a694ba00080f0e08445bd5cd47bd Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Thu, 3 Aug 2023 23:56:06 -0400 Subject: [PATCH 10/13] ci(large models): use setup-fortran, expand matrix, fix intel caching on Windows (#1322) --- .github/workflows/large.yml | 40 +++++++++++++++++++++++-------------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/.github/workflows/large.yml b/.github/workflows/large.yml index 291dcc16673..7d90a156a64 100644 --- a/.github/workflows/large.yml +++ b/.github/workflows/large.yml @@ -3,29 +3,44 @@ on: schedule: - cron: '0 6 * * *' # run at 6 AM UTC every day jobs: + # caching only necessary on Windows cache_ifort: name: Cache Intel OneAPI compilers runs-on: ${{ matrix.os }} strategy: fail-fast: false matrix: - os: [ ubuntu-22.04, macos-12, windows-2022 ] + include: + # ifx + - {os: windows-2022, compiler: intel, version: 2022.2} + # ifort + - {os: windows-2022, compiler: intel-classic, version: "2021.10"} + - {os: windows-2022, compiler: intel-classic, version: 2021.9} + - {os: windows-2022, compiler: intel-classic, version: 2021.8} + - {os: windows-2022, compiler: intel-classic, version: 2021.7} + - {os: windows-2022, compiler: intel-classic, version: 2021.6} steps: - - name: Setup Intel Fortran - uses: modflowpy/install-intelfortran-action@v1 + - name: Setup ${{ matrix.compiler }} ${{ matrix.version }} + uses: awvwgk/setup-fortran@main + with: + compiler: ${{ matrix.compiler }} + version: ${{ matrix.version }} test: name: Test runs-on: ubuntu-22.04 strategy: fail-fast: false matrix: - fc: [ ifort, gfortran ] - repo: [ examples, largetestmodels ] + include: + - {compiler: gcc, version: 13, repo: examples} + - {compiler: gcc, version: 13, repo: largetestmodels} + - {compiler: intel, version: 2022.2.1, repo: examples} + - {compiler: intel, version: 2022.2.1, repo: largetestmodels} + - {compiler: intel-classic, version: 2021.6, repo: examples} + - {compiler: intel-classic, version: 2021.6, repo: largetestmodels} defaults: run: shell: bash -l {0} - env: - GCC_V: 12 steps: - name: Checkout modflow6 @@ -46,16 +61,11 @@ jobs: cache-downloads: true cache-environment: true - - name: Setup gfortran ${{ env.GCC_V }} - if: matrix.FC == 'gfortran' + - name: Setup compilers (${{ matrix.compiler }} ${{ matrix.version }}) uses: awvwgk/setup-fortran@main with: - compiler: gcc - version: ${{ env.GCC_V }} - - - name: Setup ifort - if: matrix.fc == 'ifort' - uses: modflowpy/install-intelfortran-action@v1 + compiler: ${{ matrix.compiler }} + version: ${{ matrix.version }} - name: Cache modflow6 examples id: cache-examples From 5d4f56a5c785c9dd5681032f2a6c8b766533d9af Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Sat, 5 Aug 2023 14:04:57 -0400 Subject: [PATCH 11/13] test: define targets in conftest.py, refactor devtools usages (#1327) * configure exes in conftest.py instead of devtools * remove usages of Case (to be removed from devtools) * refactor test_gwf_maw[1-4] to accommodate removal --- autotest/conftest.py | 42 +- autotest/test_gwf.py | 27 -- autotest/test_gwf_maw01.py | 213 +++++++++ autotest/test_gwf_maw02.py | 312 +++++++++++++ autotest/test_gwf_maw03.py | 229 +++++++++ autotest/test_gwf_maw04.py | 366 ++++++++------- autotest/test_gwf_maw05.py | 2 +- autotest/test_gwf_maw06.py | 1 + autotest/test_gwf_maw07.py | 1 + autotest/test_gwf_maw08.py | 1 + autotest/test_gwf_maw09.py | 1 + autotest/test_gwf_maw10.py | 2 +- autotest/test_gwf_maw_cases.py | 817 --------------------------------- autotest/test_gwf_maw_obs.py | 1 - 14 files changed, 979 insertions(+), 1036 deletions(-) delete mode 100644 autotest/test_gwf.py create mode 100644 autotest/test_gwf_maw01.py create mode 100644 autotest/test_gwf_maw02.py create mode 100644 autotest/test_gwf_maw03.py delete mode 100644 autotest/test_gwf_maw_cases.py diff --git a/autotest/conftest.py b/autotest/conftest.py index 5f0073b18c2..27ff70c453d 100644 --- a/autotest/conftest.py +++ b/autotest/conftest.py @@ -1,8 +1,11 @@ import platform +import sys from pathlib import Path +from os import PathLike +from typing import Dict import pytest -from modflow_devtools.executables import Executables, build_default_exe_dict +from modflow_devtools.executables import Executables, get_suffixes pytest_plugins = ["modflow_devtools.fixtures"] project_root_path = Path(__file__).resolve().parent.parent @@ -45,7 +48,42 @@ def libmf6_path(bin_path) -> Path: @pytest.fixture(scope="session") def targets(bin_path) -> Executables: - return Executables(**build_default_exe_dict(bin_path)) + exe_ext, lib_ext = get_suffixes(sys.platform) + dl_bin = bin_path / "downloaded" + rb_bin = bin_path / "rebuilt" + tgts = dict() + + # downloaded executables + tgts["mf2000"] = dl_bin / f"mf2000{exe_ext}" + tgts["mf2005"] = dl_bin / f"mf2005dbl{exe_ext}" + tgts["mfnwt"] = dl_bin / f"mfnwtdbl{exe_ext}" + tgts["mfusg"] = dl_bin / f"mfusgdbl{exe_ext}" + tgts["mflgr"] = dl_bin / f"mflgrdbl{exe_ext}" + tgts["mf2005s"] = dl_bin / f"mf2005{exe_ext}" + tgts["mt3dms"] = dl_bin / f"mt3dms{exe_ext}" + tgts["crt"] = dl_bin / f"crt{exe_ext}" + tgts["gridgen"] = dl_bin / f"gridgen{exe_ext}" + tgts["mp6"] = dl_bin / f"mp6{exe_ext}" + tgts["mp7"] = dl_bin / f"mp7{exe_ext}" + tgts["swtv4"] = dl_bin / f"swtv4{exe_ext}" + tgts["sutra"] = dl_bin / f"sutra{exe_ext}" + tgts["triangle"] = dl_bin / f"triangle{exe_ext}" + tgts["vs2dt"] = dl_bin / f"vs2dt{exe_ext}" + tgts["zonbudusg"] = dl_bin / f"zonbudusg{exe_ext}" + + # binaries rebuilt from last release + tgts["mf6_regression"] = rb_bin / f"mf6{exe_ext}" + tgts["libmf6_regression"] = rb_bin / f"libmf6{lib_ext}" + tgts["mf5to6_regression"] = rb_bin / f"mf5to6{exe_ext}" + tgts["zbud6_regression"] = rb_bin / f"zbud6{exe_ext}" + + # local development binaries + tgts["mf6"] = bin_path / f"mf6{exe_ext}" + tgts["libmf6"] = bin_path / f"libmf6{lib_ext}" + tgts["mf5to6"] = bin_path / f"mf5to6{exe_ext}" + tgts["zbud6"] = bin_path / f"zbud6{exe_ext}" + + return Executables(**tgts) @pytest.fixture diff --git a/autotest/test_gwf.py b/autotest/test_gwf.py deleted file mode 100644 index 59e0d18f6ac..00000000000 --- a/autotest/test_gwf.py +++ /dev/null @@ -1,27 +0,0 @@ -from modflow_devtools.executables import Executables -from pytest_cases import parametrize_with_cases -from simulation import TestSimulation -from test_gwf_maw04 import GwfMaw04Cases -from test_gwf_maw_cases import GwfMawCases - - -@parametrize_with_cases("case", cases=[GwfMawCases, GwfMaw04Cases]) -def test_gwf_models(case, targets: Executables): - data, sim, cmp, exfunc = case - sim.write_simulation() - if cmp: - cmp.write_simulation() - - test = TestSimulation( - name=data.name, - exe_dict=targets, - exfunc=exfunc, - idxsim=0, # TODO: remove parameter from TestSimulation - mf6_regression=True, - require_failure=data.xfail, - make_comparison=data.compare, - ) - - test.set_model(sim.simulation_data.mfpath.get_sim_path(), testModel=False) - test.run() - test.compare() diff --git a/autotest/test_gwf_maw01.py b/autotest/test_gwf_maw01.py new file mode 100644 index 00000000000..7db88abe711 --- /dev/null +++ b/autotest/test_gwf_maw01.py @@ -0,0 +1,213 @@ +import os +from types import SimpleNamespace as Case + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +budtol = 1e-2 +bud_lst = ["GWF_IN", "GWF_OUT", "RATE_IN", "RATE_OUT"] + +well1 = Case( + observations={"maw_obs.csv": [("mh1", "head", 1)]}, + packagedata=[[0, 0.1, 50.0, 100.0, "THIEM", 1]], + connectiondata=[[0, 0, (0, 0, 1), 100.0, 50.0, 1.0, 0.1]], + perioddata=[[0, "rate", 0.0]], +) + +ex = ["maw01", "maw01nwt", "maw01nwtur"] +krylov = ["CG", "BICGSTAB", "BICGSTAB"] +newton = [None, "NEWTON", "NEWTON UNDER_RELAXATION"] +nlay = 1 +nrow = 1 +ncol = 3 +nper = 3 +delr = 300 +delc = 300 +perlen = 3 * [1] +nstp = 3 * [1] +tsmult = 3 * [1] +well = well1 +strt = 100 +hk = 1 +nouter = 100 +ninner = 300 +hclose = 1e-9 +rclose = 1e-3 +relaxation_factor = 1 +compare = False + + +def build_model(idx, ws, mf6): + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, + version="mf6", + exe_name=mf6, + sim_ws=ws, + ) + + # create tdis package + tdis_rc = [(perlen[i], nstp[i], tsmult[i]) for i in range(nper)] + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.MFModel( + sim, + model_type="gwf6", + modelname=name, + model_nam_file=f"{name}.nam", + ) + gwf.name_file.newtonoptions = newton[idx] + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=krylov[idx], + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relaxation_factor, + ) + sim.register_ims_package(ims, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=100.0, + botm=0.0, + idomain=1, + filename=f"{name}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=hk, + k33=hk, + filename=f"{name}.npf", + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.0, + sy=0.1, + steady_state={0: True}, + # transient={1: False}, + filename=f"{name}.sto", + ) + + # chd files + chdlist0 = [] + chdlist0.append([(0, 0, 0), 100.0]) + chdlist0.append([(0, 0, 2), 100.0]) + + chdlist1 = [] + chdlist1.append([(0, 0, 0), 25.0]) + chdlist1.append([(0, 0, 2), 25.0]) + + chdspdict = {0: chdlist0, 1: chdlist1, 2: chdlist0} + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspdict, + save_flows=False, + filename=f"{name}.chd", + ) + + # wel files + # wel = flopy.mf6.ModflowGwfwel(gwf, print_input=True, print_flows=True, + # maxbound=len(ws), + # periodrecarray=wd6, + # save_flows=False) + # MAW + maw = flopy.mf6.ModflowGwfmaw( + gwf, + filename=f"{name}.maw", + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + observations=well.observations, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{name}.oc", + ) + + return sim, None + + +def eval_results(sim): + print("evaluating MAW heads...") + + # MODFLOW 6 maw results + fpth = os.path.join(sim.simpath, "maw_obs.csv") + tc = np.genfromtxt(fpth, names=True, delimiter=",") + + # create known results array + tc0 = np.array([100.0, 25.0, 100.0]) + + # calculate maximum absolute error + diff = tc["MH1"] - tc0 + diffmax = np.abs(diff).max() + dtol = 1e-9 + msg = f"maximum absolute maw head difference ({diffmax}) " + + if diffmax > dtol: + sim.success = False + msg += f"exceeds {dtol}" + assert diffmax < dtol, msg + else: + sim.success = True + print(" " + msg) + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_maw02.py b/autotest/test_gwf_maw02.py new file mode 100644 index 00000000000..09eda8d61be --- /dev/null +++ b/autotest/test_gwf_maw02.py @@ -0,0 +1,312 @@ +import os +from types import SimpleNamespace as Case + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +budtol = 1e-2 +bud_lst = ["GWF_IN", "GWF_OUT", "RATE_IN", "RATE_OUT"] + +well2 = Case( + observations={"maw_obs.csv": [("mh1", "head", 1)]}, + packagedata=[ + [0, 0.1, 0.0, 100.0, "THIEM", 1], + [1, 0.1, 0.0, 100.0, "THIEM", 1], + ], + connectiondata=[ + [0, 0, (0, 0, 1), 100.0, 0.0, 1.0, 0.1], + [1, 0, (0, 0, 1), 100.0, 0.0, 1.0, 0.1], + ], + perioddata={ + 0: [ + [0, "rate", -20.0], + [0, "status", "inactive"], + [0, "rate_scaling", 1.0, 15.0], + [1, "rate", -30.0], + [1, "status", "inactive"], + [1, "rate_scaling", 5.0, 15.0], + ], + 1: [ + [0, "rate", -110.0], + [0, "status", "active"], + [1, "rate", -130.0], + [1, "status", "active"], + ], + 3: [[0, "status", "inactive"]], + 4: [[0, "status", "active"]], + }, +) + +ex = ["maw02"] +krylov = "CG" +nlay = 1 +nrow = 1 +ncol = 3 +nper = 5 +delr = 300 +delc = 300 +perlen = 5 * [1] +nstp = 5 * [1] +tsmult = 5 * [1] +well = well2 +strt = 100 +hk = 1 +nouter = 100 +ninner = 300 +hclose = 1e-9 +rclose = 1e-3 +relaxation_factor = 1 +compare = False + + +def build_model(idx, ws, mf6): + name = ex[idx] + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name=mf6, sim_ws=ws + ) + + # create tdis package + tdis_rc = [(perlen[i], nstp[i], tsmult[i]) for i in range(nper)] + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.MFModel( + sim, + model_type="gwf6", + modelname=name, + model_nam_file=f"{name}.nam", + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=krylov, + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relaxation_factor, + ) + sim.register_ims_package(ims, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=100.0, + botm=0.0, + idomain=1, + filename=f"{name}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=hk, + k33=hk, + filename=f"{name}.npf", + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=1, + ss=0.0, + sy=0.1, + steady_state={0: True}, + # transient={1: False}, + filename=f"{name}.sto", + ) + + # chd files + chdlist0 = [] + chdlist0.append([(0, 0, 0), 100.0]) + chdlist0.append([(0, 0, 2), 100.0]) + + chdlist1 = [] + chdlist1.append([(0, 0, 0), 25.0]) + chdlist1.append([(0, 0, 2), 25.0]) + + chdspdict = {0: chdlist0, 1: chdlist1, 2: chdlist0} + chd = flopy.mf6.ModflowGwfchd( + gwf, + stress_period_data=chdspdict, + save_flows=False, + filename=f"{name}.chd", + ) + + # MAW + maw = flopy.mf6.ModflowGwfmaw( + gwf, + filename=f"{name}.maw", + budget_filerecord=f"{name}.maw.cbc", + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + observations=well.observations, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + pname="MAW-1", + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")], + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{name}.oc", + ) + + return sim, None + + +def eval_results(sim): + print("evaluating MAW budgets...") + + shape3d = (nlay, nrow, ncol) + size3d = nlay * nrow * ncol + + # get results from listing file + fpth = os.path.join(sim.simpath, f"{os.path.basename(sim.name)}.lst") + budl = flopy.utils.Mf6ListBudget( + fpth, budgetkey="MAW-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP" + ) + names = list(bud_lst) + d0 = budl.get_budget(names=names)[0] + dtype = d0.dtype + nbud = d0.shape[0] + + # get results from cbc file + cbc_bud = ["GWF", "RATE"] + d = np.recarray(nbud, dtype=dtype) + for key in bud_lst: + d[key] = 0.0 + fpth = os.path.join(sim.simpath, f"{os.path.basename(sim.name)}.maw.cbc") + cobj = flopy.utils.CellBudgetFile(fpth, precision="double") + kk = cobj.get_kstpkper() + times = cobj.get_times() + cbc_vals = [] + for idx, (k, t) in enumerate(zip(kk, times)): + for text in cbc_bud: + qin = 0.0 + qout = 0.0 + v = cobj.get_data(kstpkper=k, text=text)[0] + if isinstance(v, np.recarray): + vt = np.zeros(size3d, dtype=float) + wq = [] + for jdx, node in enumerate(v["node"]): + vt[node - 1] += v["q"][jdx] + wq.append(v["q"][jdx]) + v = vt.reshape(shape3d) + if text == cbc_bud[-1]: + cbc_vals.append(wq) + for kk in range(v.shape[0]): + for ii in range(v.shape[1]): + for jj in range(v.shape[2]): + vv = v[kk, ii, jj] + if vv < 0.0: + qout -= vv + else: + qin += vv + d["totim"][idx] = t + d["time_step"][idx] = k[0] + d["stress_period"] = k[1] + key = f"{text}_IN" + d[key][idx] = qin + key = f"{text}_OUT" + d[key][idx] = qout + + maw_vals = [ + [0.000, 0.000], + [-106.11303563809453, -96.22598985147631], + [-110.000, -130.000], + [0.0, -130.000], + [-110.000, -130.000], + ] + + # evaluate if well rates in cbc file are equal to expected values + diffv = [] + for ovs, svs in zip(maw_vals, cbc_vals): + for ov, sv in zip(ovs, svs): + diffv.append(ov - sv) + diffv = np.abs(np.array(diffv)).max() + msg = f"\nmaximum absolute maw rate difference ({diffv})\n" + + # calculate difference between water budget items in the lst and cbc files + diff = np.zeros((nbud, len(bud_lst)), dtype=float) + for idx, key in enumerate(bud_lst): + diff[:, idx] = d0[key] - d[key] + diffmax = np.abs(diff).max() + msg += f"maximum absolute total-budget difference ({diffmax}) " + + # write summary + fpth = os.path.join( + sim.simpath, f"{os.path.basename(sim.name)}.bud.cmp.out" + ) + f = open(fpth, "w") + for i in range(diff.shape[0]): + if i == 0: + line = f"{'TIME':>10s}" + for idx, key in enumerate(bud_lst): + line += f"{key + '_LST':>25s}" + line += f"{key + '_CBC':>25s}" + line += f"{key + '_DIF':>25s}" + f.write(line + "\n") + line = f"{d['totim'][i]:10g}" + for idx, key in enumerate(bud_lst): + line += f"{d0[key][i]:25g}" + line += f"{d[key][i]:25g}" + line += f"{diff[i, idx]:25g}" + f.write(line + "\n") + f.close() + + if diffmax > budtol or diffv > budtol: + sim.success = False + msg += f"\n...exceeds {budtol}" + assert diffmax < budtol, msg + else: + sim.success = True + print(" " + msg) + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_maw03.py b/autotest/test_gwf_maw03.py new file mode 100644 index 00000000000..4a4c6ae71a8 --- /dev/null +++ b/autotest/test_gwf_maw03.py @@ -0,0 +1,229 @@ +import os +from types import SimpleNamespace as Case + +import flopy +import numpy as np +import pytest + +from framework import TestFramework +from simulation import TestSimulation + +budtol = 1e-2 +bud_lst = ["GWF_IN", "GWF_OUT", "RATE_IN", "RATE_OUT"] + + +def well3(name): + perioddata = { + "maw03a": [ + (0, "rate", 2000.0), + ], + "maw03b": [(0, "rate", 2000.0), (0, "head_limit", 0.4)], + "maw03c": [(0, "rate", 2000.0), (0, "rate_scaling", 0.0, 1.0)], + } + wellbottom = -1000 + return Case( + observations={ + f"{name}.maw.obs.csv": [ + ("m1head", "head", (0,)), + ("m1rate", "rate", (0,)), + ] # is this index one-based? Not if in a tuple + }, + packagedata=[[0, 0.15, wellbottom, 0.0, "THIEM", 1]], + connectiondata=[[0, 0, (0, 50, 50), 0.0, wellbottom, 0.0, 0.0]], + perioddata=perioddata[name], + ) + + +ex = ["maw03a", "maw03b", "maw03c"] +krylov = "CG" +nlay = 1 +nrow = 101 +ncol = 101 +nper = 1 +delr = 142 +delc = 142 +perlen = [1000] +nstp = [50] +tsmult = [1.2] +strt = 0 +hk = 10 +nouter = 100 +ninner = 100 +hclose = 1e-6 +rclose = 1e-6 +relaxation_factor = 1 +compare = False + + +def build_model(idx, ws, mf6): + top = 0.0 + botm = [-1000.0] + + tdis_rc = [] + for i in range(nper): + tdis_rc.append((perlen[i], nstp[i], tsmult[i])) + + name = ex[idx] + sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name=mf6) + + # create tdis package + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) + + # create gwf model + gwf = flopy.mf6.MFModel( + sim, + model_type="gwf6", + modelname=name, + model_nam_file=f"{name}.nam", + ) + + # create iterative model solution and register the gwf model with it + ims = flopy.mf6.ModflowIms( + sim, + print_option="SUMMARY", + outer_dvclose=hclose, + outer_maximum=nouter, + under_relaxation="NONE", + inner_maximum=ninner, + inner_dvclose=hclose, + rcloserecord=rclose, + linear_acceleration=krylov, + scaling_method="NONE", + reordering_method="NONE", + relaxation_factor=relaxation_factor, + ) + sim.register_ims_package(ims, [gwf.name]) + + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=top, + botm=botm, + idomain=1, + filename=f"{name}.dis", + ) + + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt, filename=f"{name}.ic") + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, + save_flows=True, + icelltype=1, + k=hk, + k33=hk, + filename=f"{name}.npf", + ) + + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=True, + iconvert=0, + ss=1.0e-5, + sy=0.1, + steady_state={0: False}, + transient={0: True}, + filename=f"{name}.sto", + ) + + # MAW + well = well3(name) + maw = flopy.mf6.ModflowGwfmaw( + gwf, + filename=f"{name}.maw", + print_input=True, + print_head=True, + print_flows=True, + save_flows=True, + observations=well.observations, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + ) + + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + headprintrecord=[ + ("COLUMNS", ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL") + ], + saverecord=[("HEAD", "ALL")], + printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + filename=f"{name}.oc", + ) + + # head observations + obs_data0 = [("head_well_cell", "HEAD", (0, 0, 0))] + obs_recarray = {f"{name}.obs.csv": obs_data0} + obs = flopy.mf6.ModflowUtlobs( + gwf, + pname="head_obs", + filename=f"{name}.obs", + digits=15, + print_input=True, + continuous=obs_recarray, + ) + + return sim, None + + +def eval_results(sim): + print("evaluating MAW heads...") + + # MODFLOW 6 maw results + test_name = sim.name + fpth = os.path.join(sim.simpath, f"{test_name}.maw.obs.csv") + tc = np.genfromtxt(fpth, names=True, delimiter=",") + + if test_name.endswith("a"): + # M1RATE should be 2000. + msg = "The injection rate should be 2000. for all times" + assert tc["M1RATE"].min() == tc["M1RATE"].max() == 2000, msg + + elif test_name.endswith("b"): + # M1RATE should have a minimum value less than 200 and + # M1HEAD should not exceed 0.400001 + msg = ( + "Injection rate should fall below 200 and the head should not" + "exceed 0.4" + ) + assert tc["M1RATE"].min() < 200.0, msg + assert tc["M1HEAD"].max() < 0.400001, msg + + elif test_name.endswith("c"): + # M1RATE should have a minimum value less than 800 + # M1HEAD should not exceed 1.0. + msg = ( + "Min injection rate should be less than 800 and well " + "head should not exceed 1.0" + ) + assert tc["M1RATE"].min() < 800.0 and tc["M1HEAD"].max() < 1.0, msg + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_maw04.py b/autotest/test_gwf_maw04.py index e77656f0c68..462a3afa735 100644 --- a/autotest/test_gwf_maw04.py +++ b/autotest/test_gwf_maw04.py @@ -1,10 +1,12 @@ import os -from types import SimpleNamespace +from types import SimpleNamespace as Case import flopy import numpy as np -from modflow_devtools.case import Case -from pytest_cases import parametrize +import pytest + +from framework import TestFramework +from simulation import TestSimulation # temporal discretization nper = 2 @@ -92,7 +94,7 @@ def well4(label): [0, 1, (1, nhalf, nhalf), botm[0], botm[1], hks, sradius[label]], ] perioddata = {1: [[0, "RATE", wellq]]} - return SimpleNamespace( + return Case( print_input=True, no_well_storage=True, packagedata=packagedata, @@ -101,198 +103,188 @@ def well4(label): ) -case = Case( - name="maw_iss305", - nlay=nlay, - nrow=nrow, - ncol=ncol, - nper=nper, - delr=delr, - perlen=perlen, - nstp=nstp, - tsmult=tsmult, - steady=steady, - strt=0, - hk=10, - nouter=100, - ninner=100, - hclose=1e-9, - rclose=1e-6, - relax=1, - top=top, - botm=botm, - confined=confined, - ss=ss, - chd_spd=chd_spd, - chd5_spd=chd5_spd, - nhalf=nhalf, - radius=radius, - wellq=wellq, - compare=False, -) -cases = [case.copy_update(name=case.name + "a", well=well4("a"),)] + [ - case.copy_update(name=case.name + label, well=well4(label), xfail=True) - for label in [ - "b", - # "c", # todo: this one passes when it should fail - "d", - "e", - "f", - ] -] +# npf data +strt = 0 +hk = 10 +# solver +nouter = 100 +ninner = 100 +hclose = 1e-9 +rclose = 1e-6 +relax = 1 -class GwfMaw04Cases: - @parametrize(data=cases, ids=[c.name for c in cases]) - def case_4(self, function_tmpdir, targets, data): - name = data.name - ws = str(function_tmpdir) +# subproblems +subprobs = ["a", "b", "c", "d", "e", "f"] +ex = [f"maw_iss305{sp}" for sp in subprobs] +wells = [well4(sp) for sp in subprobs] +xfail = [False, True, True, True, True, True] - # build MODFLOW 6 files - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name=targets["mf6"], sim_ws=ws - ) - # create tdis package - tdis_rc = [] - for idx in range(data.nper): - tdis_rc.append( - (data.perlen[idx], data.nstp[idx], data.tsmult[idx]) - ) - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", nper=data.nper, perioddata=tdis_rc - ) - # create iterative model solution - ims = flopy.mf6.ModflowIms( - sim, - inner_dvclose=data.hclose, - rcloserecord=data.rclose, - outer_dvclose=data.hclose, - ) +def build_model(idx, ws, mf6): + name = ex[idx] + well = wells[idx] + # build MODFLOW 6 files + sim = flopy.mf6.MFSimulation( + sim_name=name, version="mf6", exe_name=mf6, sim_ws=ws + ) + # create tdis package + tdis_rc = [] + for kper in range(nper): + tdis_rc.append((perlen[kper], nstp[kper], tsmult[kper])) + tdis = flopy.mf6.ModflowTdis( + sim, time_units="DAYS", nper=nper, perioddata=tdis_rc + ) - # create gwf model - gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) + # create iterative model solution + ims = flopy.mf6.ModflowIms( + sim, + inner_dvclose=hclose, + rcloserecord=rclose, + outer_dvclose=hclose, + ) - # discretization - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=data.nlay, - nrow=data.nrow, - ncol=data.ncol, - delr=data.delr, - delc=data.delr, - top=data.top, - botm=data.botm, - ) - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=data.strt) + # create gwf model + gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, save_flows=False, icelltype=data.confined, k=data.hk - ) - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=False, - iconvert=data.confined, - ss=data.ss, - steady_state={0: True}, - transient={1: True}, + # discretization + dis = flopy.mf6.ModflowGwfdis( + gwf, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delr, + top=top, + botm=botm, + ) + # initial conditions + ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) + + # node property flow + npf = flopy.mf6.ModflowGwfnpf( + gwf, save_flows=False, icelltype=confined, k=hk + ) + # storage + sto = flopy.mf6.ModflowGwfsto( + gwf, + save_flows=False, + iconvert=confined, + ss=ss, + steady_state={0: True}, + transient={1: True}, + ) + # constant head + chd = flopy.mf6.ModflowGwfchd( + gwf, stress_period_data=chd_spd, save_flows=False + ) + # multi-aquifer well + maw = flopy.mf6.ModflowGwfmaw( + gwf, + print_input=well.print_input, + no_well_storage=well.no_well_storage, + packagedata=well.packagedata, + connectiondata=well.connectiondata, + perioddata=well.perioddata, + ) + # output control + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=f"{name}.cbc", + head_filerecord=f"{name}.hds", + saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], + ) + # build MODFLOW-2005 files + if xfail[idx]: + mc = None + else: + cmppth = "mf2005" + ws = os.path.join(ws, cmppth) + mc = flopy.modflow.Modflow(name, model_ws=ws, version=cmppth) + dis = flopy.modflow.ModflowDis( + mc, + nlay=nlay, + nrow=nrow, + ncol=ncol, + nper=nper, + perlen=perlen, + nstp=nstp, + tsmult=tsmult, + steady=steady, + delr=delr, + delc=delr, + top=top, + botm=botm, ) - # constant head - chd = flopy.mf6.ModflowGwfchd( - gwf, stress_period_data=data.chd_spd, save_flows=False + bas = flopy.modflow.ModflowBas(mc, strt=strt) + lpf = flopy.modflow.ModflowLpf( + mc, + laytyp=confined, + hk=hk, + vka=hk, + ss=ss, + sy=0, ) - # multi-aquifer well - maw = flopy.mf6.ModflowGwfmaw( - gwf, - print_input=data.well.print_input, - no_well_storage=data.well.no_well_storage, - packagedata=data.well.packagedata, - connectiondata=data.well.connectiondata, - perioddata=data.well.perioddata, + chd = flopy.modflow.ModflowChd(mc, stress_period_data=chd5_spd) + # mnw2 + # empty mnw2 file to create recarrays + mnw2 = flopy.modflow.ModflowMnw2(mc) + node_data = mnw2.get_empty_node_data(2) + node_data["ztop"] = np.array([top, botm[0]]) + node_data["zbotm"] = np.array([botm[0], botm[1]]) + node_data["i"] = np.array([nhalf, nhalf]) + node_data["j"] = np.array([nhalf, nhalf]) + node_data["wellid"] = np.array(["well1", "well1"]) + node_data["losstype"] = np.array(["skin", "skin"]) + node_data["rw"] = np.array([radius, radius]) + node_data["rskin"] = np.array([sradius[name[-1]], sradius[name[-1]]]) + hks = hk * skin_mult[name[-1]] + node_data["kskin"] = np.array([hks, hks]) + dtype = [("wellid", np.unicode_, 20), ("qdes", " dtol: - config.success = False - msg += f"exceeds {dtol}" - assert diffmax < dtol, msg - else: - config.success = True - print(" " + msg) - - case2 = Case( - name="maw02", - krylov="CG", - nlay=1, - nrow=1, - ncol=3, - nper=5, - delr=300, - delc=300, - perlen=5 * [1], - nstp=5 * [1], - tsmult=5 * [1], - well=well2, - strt=100, - hk=1, - nouter=100, - ninner=300, - hclose=1e-9, - rclose=1e-3, - relaxation_factor=1, - compare=False, - ) - cases2 = [case2] - - @parametrize(data=cases2, ids=[c.name for c in cases2]) - def case_2(self, function_tmpdir, targets, data): - name = data.name - ws = str(function_tmpdir) - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name=targets["mf6"], sim_ws=ws - ) - - # create tdis package - tdis_rc = [ - (data.perlen[i], data.nstp[i], data.tsmult[i]) - for i in range(data.nper) - ] - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", nper=data.nper, perioddata=tdis_rc - ) - - # create gwf model - gwf = flopy.mf6.MFModel( - sim, - model_type="gwf6", - modelname=name, - model_nam_file=f"{name}.nam", - ) - - # create iterative model solution and register the gwf model with it - ims = flopy.mf6.ModflowIms( - sim, - print_option="SUMMARY", - outer_dvclose=data.hclose, - outer_maximum=data.nouter, - under_relaxation="NONE", - inner_maximum=data.ninner, - inner_dvclose=data.hclose, - rcloserecord=data.rclose, - linear_acceleration=data.krylov, - scaling_method="NONE", - reordering_method="NONE", - relaxation_factor=data.relaxation_factor, - ) - sim.register_ims_package(ims, [gwf.name]) - - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=data.nlay, - nrow=data.nrow, - ncol=data.ncol, - delr=data.delr, - delc=data.delc, - top=100.0, - botm=0.0, - idomain=1, - filename=f"{name}.dis", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=data.strt, filename=f"{name}.ic") - - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, - save_flows=True, - icelltype=1, - k=data.hk, - k33=data.hk, - filename=f"{name}.npf", - ) - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=True, - iconvert=1, - ss=0.0, - sy=0.1, - steady_state={0: True}, - # transient={1: False}, - filename=f"{name}.sto", - ) - - # chd files - chdlist0 = [] - chdlist0.append([(0, 0, 0), 100.0]) - chdlist0.append([(0, 0, 2), 100.0]) - - chdlist1 = [] - chdlist1.append([(0, 0, 0), 25.0]) - chdlist1.append([(0, 0, 2), 25.0]) - - chdspdict = {0: chdlist0, 1: chdlist1, 2: chdlist0} - chd = flopy.mf6.ModflowGwfchd( - gwf, - stress_period_data=chdspdict, - save_flows=False, - filename=f"{name}.chd", - ) - - # MAW - maw = flopy.mf6.ModflowGwfmaw( - gwf, - filename=f"{name}.maw", - budget_filerecord=f"{name}.maw.cbc", - print_input=True, - print_head=True, - print_flows=True, - save_flows=True, - observations=data.well.observations, - packagedata=data.well.packagedata, - connectiondata=data.well.connectiondata, - perioddata=data.well.perioddata, - pname="MAW-1", - ) - - # output control - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{name}.cbc", - head_filerecord=f"{name}.hds", - headprintrecord=[ - ("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL") - ], - saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - filename=f"{name}.oc", - ) - - return data, sim, None, self.eval_2 - - def eval_2(self, config, data): - print("evaluating MAW budgets...") - - shape3d = (data.nlay, data.nrow, data.ncol) - size3d = data.nlay * data.nrow * data.ncol - - # get results from listing file - fpth = os.path.join( - config.simpath, f"{os.path.basename(config.name)}.lst" - ) - budl = flopy.utils.Mf6ListBudget( - fpth, budgetkey="MAW-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP" - ) - names = list(self.bud_lst) - d0 = budl.get_budget(names=names)[0] - dtype = d0.dtype - nbud = d0.shape[0] - - # get results from cbc file - cbc_bud = ["GWF", "RATE"] - d = np.recarray(nbud, dtype=dtype) - for key in self.bud_lst: - d[key] = 0.0 - fpth = os.path.join( - config.simpath, f"{os.path.basename(config.name)}.maw.cbc" - ) - cobj = flopy.utils.CellBudgetFile(fpth, precision="double") - kk = cobj.get_kstpkper() - times = cobj.get_times() - cbc_vals = [] - for idx, (k, t) in enumerate(zip(kk, times)): - for text in cbc_bud: - qin = 0.0 - qout = 0.0 - v = cobj.get_data(kstpkper=k, text=text)[0] - if isinstance(v, np.recarray): - vt = np.zeros(size3d, dtype=float) - wq = [] - for jdx, node in enumerate(v["node"]): - vt[node - 1] += v["q"][jdx] - wq.append(v["q"][jdx]) - v = vt.reshape(shape3d) - if text == cbc_bud[-1]: - cbc_vals.append(wq) - for kk in range(v.shape[0]): - for ii in range(v.shape[1]): - for jj in range(v.shape[2]): - vv = v[kk, ii, jj] - if vv < 0.0: - qout -= vv - else: - qin += vv - d["totim"][idx] = t - d["time_step"][idx] = k[0] - d["stress_period"] = k[1] - key = f"{text}_IN" - d[key][idx] = qin - key = f"{text}_OUT" - d[key][idx] = qout - - maw_vals = [ - [0.000, 0.000], - [-106.11303563809453, -96.22598985147631], - [-110.000, -130.000], - [0.0, -130.000], - [-110.000, -130.000], - ] - - # evaluate if well rates in cbc file are equal to expected values - diffv = [] - for ovs, svs in zip(maw_vals, cbc_vals): - for ov, sv in zip(ovs, svs): - diffv.append(ov - sv) - diffv = np.abs(np.array(diffv)).max() - msg = f"\nmaximum absolute maw rate difference ({diffv})\n" - - # calculate difference between water budget items in the lst and cbc files - diff = np.zeros((nbud, len(self.bud_lst)), dtype=float) - for idx, key in enumerate(self.bud_lst): - diff[:, idx] = d0[key] - d[key] - diffmax = np.abs(diff).max() - msg += f"maximum absolute total-budget difference ({diffmax}) " - - # write summary - fpth = os.path.join( - config.simpath, f"{os.path.basename(config.name)}.bud.cmp.out" - ) - f = open(fpth, "w") - for i in range(diff.shape[0]): - if i == 0: - line = f"{'TIME':>10s}" - for idx, key in enumerate(self.bud_lst): - line += f"{key + '_LST':>25s}" - line += f"{key + '_CBC':>25s}" - line += f"{key + '_DIF':>25s}" - f.write(line + "\n") - line = f"{d['totim'][i]:10g}" - for idx, key in enumerate(self.bud_lst): - line += f"{d0[key][i]:25g}" - line += f"{d[key][i]:25g}" - line += f"{diff[i, idx]:25g}" - f.write(line + "\n") - f.close() - - if diffmax > self.budtol or diffv > self.budtol: - config.success = False - msg += f"\n...exceeds {self.budtol}" - assert diffmax < self.budtol, msg - else: - config.success = True - print(" " + msg) - - case3 = Case( - name="maw03", - krylov="CG", - nlay=1, - nrow=101, - ncol=101, - nper=1, - delr=142, - delc=142, - perlen=[1000], - nstp=[50], - tsmult=[1.2], - strt=0, - hk=10, - nouter=100, - ninner=100, - hclose=1e-6, - rclose=1e-6, - relaxation_factor=1, - compare=False, - ) - cases3 = [ - case3.copy_update( - name="maw03a", - well=well3("maw03a"), - ), - case3.copy_update(name="maw03b", well=well3("maw03b")), - case3.copy_update(name="maw03c", well=well3("maw03c")), - ] - - @parametrize(data=cases3, ids=[c.name for c in cases3]) - def case_3(self, function_tmpdir, targets, data): - top = 0.0 - botm = [-1000.0] - - tdis_rc = [] - for i in range(data.nper): - tdis_rc.append((data.perlen[i], data.nstp[i], data.tsmult[i])) - - name = data.name - ws = str(function_tmpdir) - sim = flopy.mf6.MFSimulation( - sim_name=name, sim_ws=ws, exe_name=targets["mf6"] - ) - - # create tdis package - tdis = flopy.mf6.ModflowTdis( - sim, time_units="DAYS", nper=data.nper, perioddata=tdis_rc - ) - - # create gwf model - gwf = flopy.mf6.MFModel( - sim, - model_type="gwf6", - modelname=name, - model_nam_file=f"{name}.nam", - ) - - # create iterative model solution and register the gwf model with it - ims = flopy.mf6.ModflowIms( - sim, - print_option="SUMMARY", - outer_dvclose=data.hclose, - outer_maximum=data.nouter, - under_relaxation="NONE", - inner_maximum=data.ninner, - inner_dvclose=data.hclose, - rcloserecord=data.rclose, - linear_acceleration=data.krylov, - scaling_method="NONE", - reordering_method="NONE", - relaxation_factor=data.relaxation_factor, - ) - sim.register_ims_package(ims, [gwf.name]) - - dis = flopy.mf6.ModflowGwfdis( - gwf, - nlay=data.nlay, - nrow=data.nrow, - ncol=data.ncol, - delr=data.delr, - delc=data.delc, - top=top, - botm=botm, - idomain=1, - filename=f"{name}.dis", - ) - - # initial conditions - ic = flopy.mf6.ModflowGwfic(gwf, strt=data.strt, filename=f"{name}.ic") - - # node property flow - npf = flopy.mf6.ModflowGwfnpf( - gwf, - save_flows=True, - icelltype=1, - k=data.hk, - k33=data.hk, - filename=f"{name}.npf", - ) - - # storage - sto = flopy.mf6.ModflowGwfsto( - gwf, - save_flows=True, - iconvert=0, - ss=1.0e-5, - sy=0.1, - steady_state={0: False}, - transient={0: True}, - filename=f"{name}.sto", - ) - - # MAW - maw = flopy.mf6.ModflowGwfmaw( - gwf, - filename=f"{name}.maw", - print_input=True, - print_head=True, - print_flows=True, - save_flows=True, - observations=data.well.observations, - packagedata=data.well.packagedata, - connectiondata=data.well.connectiondata, - perioddata=data.well.perioddata, - ) - - # output control - oc = flopy.mf6.ModflowGwfoc( - gwf, - budget_filerecord=f"{name}.cbc", - head_filerecord=f"{name}.hds", - headprintrecord=[ - ("COLUMNS", data.ncol, "WIDTH", 15, "DIGITS", 6, "GENERAL") - ], - saverecord=[("HEAD", "ALL")], - printrecord=[("HEAD", "ALL"), ("BUDGET", "ALL")], - filename=f"{name}.oc", - ) - - # head observations - obs_data0 = [("head_well_cell", "HEAD", (0, 0, 0))] - obs_recarray = {f"{name}.obs.csv": obs_data0} - obs = flopy.mf6.ModflowUtlobs( - gwf, - pname="head_obs", - filename=f"{name}.obs", - digits=15, - print_input=True, - continuous=obs_recarray, - ) - - return data, sim, None, self.eval_3 - - def eval_3(self, config, data): - print("evaluating MAW heads...") - - # MODFLOW 6 maw results - test_name = config.name - case_name = data.name - fpth = os.path.join(config.simpath, f"{test_name}.maw.obs.csv") - tc = np.genfromtxt(fpth, names=True, delimiter=",") - - if case_name == "a": - - # M1RATE should be 2000. - msg = "The injection rate should be 2000. for all times" - assert tc["M1RATE"].min() == tc["M1RATE"].max() == 2000, msg - - elif case_name == "b": - - # M1RATE should have a minimum value less than 200 and - # M1HEAD should not exceed 0.400001 - msg = ( - "Injection rate should fall below 200 and the head should not" - "exceed 0.4" - ) - assert tc["M1RATE"].min() < 200.0, msg - assert tc["M1HEAD"].max() < 0.400001, msg - - elif case_name == "c": - - # M1RATE should have a minimum value less than 800 - # M1HEAD should not exceed 1.0. - msg = ( - "Min injection rate should be less than 800 and well " - "head should not exceed 1.0" - ) - assert tc["M1RATE"].min() < 800.0 and tc["M1HEAD"].max() < 1.0, msg - - # TODO - # case4 = Case( - # name='maw_iss305', - # krylov='CG', - # nlay=2, - # nrow=101, - # ncol=101, - # nper=2, - # delr=142, - # delc=142, - # perlen=[0.0, 365.0], - # nstp=[1, 25], - # tsmult=[1.0, 1.1], - # steady=[True, False], - # well=None, - # strt=0, - # hk=10, - # nouter=100, - # ninner=100, - # hclose=1e-9, - # rclose=1e-6, - # relax=1, - # require_failure=True - # ) - # cases4 = [ - # case4.copy_update({ - # 'name': "maw_iss305a", - # 'well': well4(case4, "CUMULATIVE"), - # 'require_failure': False - # }), - # case4.copy_update({ - # 'name': "maw_iss305b", - # 'well': well4(case4, "SKIN") - # }), - # case4.copy_update({ - # 'name': "maw_iss305c", - # 'well': well4(case4, "SKIN") - # }), - # case4.copy_update({ - # 'name': "maw_iss305d", - # 'well': well4(case4, "SKIN") - # }), - # case4.copy_update({ - # 'name': "maw_iss305e", - # 'well': well4(case4, "SPECIFIED") - # }), - # case4.copy_update({ - # 'name': "maw_iss305f", - # 'well': well4(case4, "CUMULATIVE") - # }) - # ] - - # @parametrize(data=cases4, ids=[c['name'] for c in cases4]) - # def case_4(self, function_tmpdir, targets, data): - # pass - - # def eval_4(self, sim, data): - # pass diff --git a/autotest/test_gwf_maw_obs.py b/autotest/test_gwf_maw_obs.py index 7455949fc37..050de4f96b9 100644 --- a/autotest/test_gwf_maw_obs.py +++ b/autotest/test_gwf_maw_obs.py @@ -14,7 +14,6 @@ def build_model(dir, exe): - nlay, nrow, ncol = 1, 1, 3 nper = 3 perlen = [1.0, 1.0, 1.0] From 11089fe3cf9731b0650802fd33186c2fd3314a1c Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Mon, 7 Aug 2023 11:55:15 -0400 Subject: [PATCH 12/13] fix(build_docs.py): filter benchmark artifacts by branch (#1330) --- distribution/build_docs.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/distribution/build_docs.py b/distribution/build_docs.py index 37d5ad57fa6..353ba713a31 100644 --- a/distribution/build_docs.py +++ b/distribution/build_docs.py @@ -16,8 +16,12 @@ import pytest from flaky import flaky from modflow_devtools.build import meson_build -from modflow_devtools.download import (download_and_unzip, download_artifact, - get_release, list_artifacts) +from modflow_devtools.download import ( + download_and_unzip, + download_artifact, + get_release, + list_artifacts, +) from modflow_devtools.markers import requires_exe, requires_github from modflow_devtools.misc import is_in_ci, run_cmd, set_dir @@ -112,6 +116,9 @@ def download_benchmarks( key=lambda a: datetime.strptime(a["created_at"], "%Y-%m-%dT%H:%M:%SZ"), reverse=True, ) + artifacts = [ + a for a in artifacts if a["workflow_run"]["head_branch"] == "develop" # todo make configurable + ] most_recent = next(iter(artifacts), None) print(f"Found most recent benchmarks (artifact {most_recent['id']})") if most_recent: From 17a29c480a8408b773a0d01f74a9c2e94b6fdc9b Mon Sep 17 00:00:00 2001 From: w-bonelli Date: Tue, 8 Aug 2023 10:55:06 -0400 Subject: [PATCH 13/13] docs: add PR template, update bug report template (#1324) * add pull request template with suggested task checklist * add compiler version to bug report template (for installs built from src) --- .github/ISSUE_TEMPLATE/bug_report.md | 9 +++++---- .github/PULL_REQUEST_TEMPLATE/pull_request_template.md | 8 ++++++++ 2 files changed, 13 insertions(+), 4 deletions(-) create mode 100644 .github/PULL_REQUEST_TEMPLATE/pull_request_template.md diff --git a/.github/ISSUE_TEMPLATE/bug_report.md b/.github/ISSUE_TEMPLATE/bug_report.md index bbfd8d081ab..6d7a3d18f05 100644 --- a/.github/ISSUE_TEMPLATE/bug_report.md +++ b/.github/ISSUE_TEMPLATE/bug_report.md @@ -10,7 +10,7 @@ assignees: '' **Describe the bug** A clear and concise description of what the bug is. -**To Reproduce** +**To reproduce** Steps to reproduce the behavior: 1. Go to '...' 2. Click on '....' @@ -23,6 +23,7 @@ A clear and concise description of what you expected to happen. **Screenshots** If applicable, add screenshots to help explain your problem. -**Desktop (please complete the following information):** - - OS: [e.g. macOS, Linux, Windows] - - Version [e.g. 22] +**Environment** + - Operating system (e.g. macOS, Linux, Windows) and version + - MODFLOW 6 version (if installed via distribution) + - Compiler toolchain/version (if built from source) diff --git a/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md b/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md new file mode 100644 index 00000000000..1affcbe0ed5 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE/pull_request_template.md @@ -0,0 +1,8 @@ +Feel free to remove check-list items that aren't relevant to your change. + +- [ ] Closes #xxxx +- [ ] Passed autotests +- [ ] Formatted source files with `fprettify` +- [ ] Updated definition (*.dfn) files with new or modified options +- [ ] Described new options, features or behavior changes in release notes +- [ ] Updated meson files, makefiles, and Visual Studio project files if new source files added \ No newline at end of file