diff --git a/.vscode/build_vscode.py b/.vscode/build_vscode.py index d222cc7c772..208d0fdc538 100644 --- a/.vscode/build_vscode.py +++ b/.vscode/build_vscode.py @@ -1,6 +1,7 @@ import subprocess import os import argparse +import platform import shutil import shlex @@ -11,7 +12,7 @@ args = parser.parse_args() os.environ["FC"] = args.compiler -builddir = f"builddir_{args.compiler}_{args.buildtype}" +builddir = f"builddir_{platform.system()}_{args.compiler}_{args.buildtype}" arg_parallel = "-Dparallel=false" if os.getenv("BUILD_PARALLEL_MF6") is not None: diff --git a/autotest/test_gwe_sfe_strmbedcond.py b/autotest/test_gwe_sfe_strmbedcond.py index 6c9deb027e8..ce390dab2e7 100644 --- a/autotest/test_gwe_sfe_strmbedcond.py +++ b/autotest/test_gwe_sfe_strmbedcond.py @@ -1,7 +1,6 @@ # Test conduction between an advanced package feature, in this case stream # reaches with varying channel geometries and the host GWE gw cells. # This test should include: -# - no gw-sw interaction # - with gw-sw interaction # - hot gw cell warming an upstream reach # - thermally hot stream water warming host gw cells @@ -623,15 +622,9 @@ def check_output(idx, test): + str(j) + "in stress period 1 does not match explicitly-calculated answer" ) - assert np.isclose(wa, shared_area[0, j], atol=1e-4), msg - - msg = ( - "Wetted streambed area of all reaches should be zero in stess " - "period 2" - ) - for val in list(sfrstg[1])[1:]: - assert val == 0.0, msg + assert np.isclose(wa, shared_area[0, j], atol=1e-4), msg + # Sub-scenario checks # initialize search term srchStr = "SFE-1 BUDGET FOR ENTIRE MODEL AT END OF TIME STEP 1, STRESS PERIOD 1" diff --git a/autotest/test_gwf_sfr_gwdischarge.py b/autotest/test_gwf_sfr_gwdischarge.py new file mode 100644 index 00000000000..16dcf16785f --- /dev/null +++ b/autotest/test_gwf_sfr_gwdischarge.py @@ -0,0 +1,120 @@ +# Test evap in SFR reaches (no interaction with gwf) + +import math +import pathlib as pl + +import flopy +import numpy as np +import pytest + +from framework import TestFramework + +cases = ["sfr-gwfout"] + + +def build_models(idx, test): + # Base simulation and model name and workspace + ws = test.workspace + name = cases[idx] + + length_units = "m" + time_units = "sec" + + nrow = 1 + ncol = 1 + nlay = 1 + delr = delc = 1.0 + + sim = flopy.mf6.MFSimulation( + sim_name=name, sim_ws=ws, exe_name="mf6", version="mf6" + ) + flopy.mf6.ModflowTdis(sim, time_units=time_units) + flopy.mf6.ModflowIms( + sim, + inner_dvclose=1e-5, + inner_hclose=1e-6, + ) + + gwf = flopy.mf6.ModflowGwf( + sim, + modelname=name, + ) + flopy.mf6.ModflowGwfdis( + gwf, + length_units=length_units, + nlay=nlay, + nrow=nrow, + ncol=ncol, + delr=delr, + delc=delc, + top=0.0, + botm=-100.0, + ) + flopy.mf6.ModflowGwfnpf( + gwf, + icelltype=1, # >0 means saturated thickness varies with computed head + ) + flopy.mf6.ModflowGwfic(gwf, strt=1.0) + flopy.mf6.ModflowGwfghb(gwf, stress_period_data=[((0, 0, 0), 1.0, 1e6)]) + + # sfr data + # + package_data = [ + (0, (0, 0, 0), delr, 1.0, 1e-3, 0.0, 1.0, 1.0, 0.001, 0, 0.0, 0) + ] + connection_data = [(0)] + + sfr_obs = { + f"{name}.sfr.obs.csv": [ + ("gwf", "sfr", (0,)), + ("outflow", "ext-outflow", (0,)), + ("depth", "depth", (0,)), + ], + "filename": name + ".sfr.obs", + } + + flopy.mf6.ModflowGwfsfr( + gwf, + save_flows=True, + print_stage=True, + print_flows=True, + print_input=True, + length_conversion=1.0, + time_conversion=1.0, + nreaches=1, + packagedata=package_data, + connectiondata=connection_data, + observations=sfr_obs, + ) + + return sim, None + + +def check_output(idx, test): + answer = np.array( + [ + 1.0, + -0.92094535738673577, + -0.92094535738673577, + 0.79053721667952215e-1, + ] + ) + obs_pth = pl.Path(f"{test.workspace}/{cases[idx]}.sfr.obs.csv") + sim_data = flopy.utils.Mf6Obs(obs_pth).get_data() + data_names = sim_data.dtype.names + for idx, name in enumerate(data_names): + assert np.allclose( + sim_data[name][0], answer[idx] + ), f"simulated sfr {name} results do not match answer" + + +@pytest.mark.parametrize("idx, name", enumerate(cases)) +def test_mf6model(idx, name, function_tmpdir, targets): + test = TestFramework( + name=name, + workspace=function_tmpdir, + targets=targets, + build=lambda t: build_models(idx, t), + check=lambda t: check_output(idx, t), + ) + test.run() diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index b35e255a65b..7af8377f5cf 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -43,6 +43,7 @@ \begin{itemize} \item A divide by zero error would occur in the Streamflow Routing package when reaches were deactivated during a simulation. This bug was fixed by checking if the downstream reach is inactive before calculating the flow to the downstream reach. \item When using the mover transport (MVT) package with UZF and UZT, rejected infiltration was paired with the calculated concentration of the UZF object rather than the user-specified concentration assigned to the infiltration. This bug was fixed by instead pairing the rejected infiltration transferred by the MVR and MVT packages with the user-specified concentration assigned in UZTSETTING with the keyword ``INFILTRATION.'' With this change, MODFLOW 6 simulations that use UZF with the ``SIMULATE\_GWSEEP'' option will not transfer this particular source of water with the correct concentration. Instead, the DRN package should be used to simulate groundwater discharge to land surface. By simulating groundwater discharge to land surface with the DRN package and rejected infiltration with the UZF package, the correct concentrations will be assigned to the water transferred by the MVR package. + \item The Streamflow Routing package would not calculate groundwater discharge to a reach in cases where the groundwater head is above the top of the reach and the inflow to the reach from upstream reaches, specified inflows, rainfall, and runoff is zero. This bug has been fixed by eliminating the requirement that the conductance calculated based on the initial calculated stage is greater than zero in order to solve for groundwater . Instead, .As a result, differences in groundwater and surface-water exchange and groundwater heads in existing models may occur. % \item xxx \end{itemize} diff --git a/src/Model/GroundWaterFlow/gwf-sfr.f90 b/src/Model/GroundWaterFlow/gwf-sfr.f90 index 7c42b81742e..70f4021d499 100644 --- a/src/Model/GroundWaterFlow/gwf-sfr.f90 +++ b/src/Model/GroundWaterFlow/gwf-sfr.f90 @@ -184,6 +184,7 @@ module SfrModule procedure, private :: sfr_solve procedure, private :: sfr_update_flows procedure, private :: sfr_calc_qgwf + procedure, private :: sfr_gwf_conn procedure, private :: sfr_calc_cond procedure, private :: sfr_calc_qman procedure, private :: sfr_calc_qd @@ -3372,7 +3373,6 @@ subroutine sfr_solve(this, n, h, hcof, rhs, update) real(DP) :: tp real(DP) :: bt real(DP) :: hsfr - real(DP) :: cstr real(DP) :: qd real(DP) :: en1 real(DP) :: en2 @@ -3499,15 +3499,11 @@ subroutine sfr_solve(this, n, h, hcof, rhs, update) qd = MAX(qsrc, DZERO) qgwf = DZERO ! - ! -- calculate reach conductance for a unit depth of water - ! if equal to zero will skip iterations - call this%sfr_calc_cond(n, d1, cstr, hsfr, hgwf) - ! ! -- set flag to skip iterations isolve = 1 if (hsfr <= tp .and. hgwf <= tp) isolve = 0 if (hgwf <= tp .and. qc < DEM30) isolve = 0 - if (cstr < DEM30) isolve = 0 + if (this%sfr_gwf_conn(n) == 0) isolve = 0 if (this%iboundpak(n) < 0) isolve = 0 ! ! -- iterate to achieve solution @@ -4100,6 +4096,28 @@ subroutine sfr_calc_qgwf(this, n, depth, hgwf, qgwf, gwfhcof, gwfrhs) return end subroutine sfr_calc_qgwf + !> @brief Determine if a reach is connected to a gwf cell + !! + !! Function to determine if a reach is connected to a gwf cell. If connected, + !! the return value is 1. Otherwise, the return value is 0. + !! + !< + function sfr_gwf_conn(this, n) + ! -- return variable + integer(I4B) :: sfr_gwf_conn !< flag indicating if reach is connected to a gwf cell + ! -- dummy variables + class(SfrType) :: this !< SfrType object + integer(I4B), intent(in) :: n !< reach number + ! -- local variables + integer(I4B) :: node + + sfr_gwf_conn = 0 + node = this%igwfnode(n) + if (node > 0 .and. this%hk(n) > DZERO) then + sfr_gwf_conn = 1 + end if + end function sfr_gwf_conn + !> @brief Calculate reach-aquifer conductance !! !! Method to calculate the reach-aquifer conductance for a SFR package reach. @@ -4125,6 +4143,7 @@ subroutine sfr_calc_cond(this, n, depth, cond, hsfr, htmp) vscratio = DONE ! ! -- calculate conductance if GWF cell is active + ! rch-gwf flow will not occur if reach connected to an constant head cell node = this%igwfnode(n) if (node > 0) then if (this%ibound(node) > 0) then