fix(spdis): add the no-flow cell boundary faces to the specific discharge calculation

autotest/test_gwf_spdis_gwfgwf.py

@@ -0,0 +1,229 @@
+"""
+Test the specific discharge calculation
+
+  1. for a DIS grid with a chain of cells extending in the y direction
+  2. that same grid split in two
+
+i.e.
+
+    [1]         [1]
+    ...         ...
+    [5]         [5]
+           vs    +
+    [6]         [1]
+    ...         ...
+    [10         [5]
+
+Case 2. comes with a non-zero angledegx in the exchange file (=270.0)
+and, as a result from rounding errors, a normal vector of the
+form (epsilon,1.0) with epsilon tiny (~1.e-16) but not zero (as
+it would be for the internal connection of the DIS grid).
+The calculation of specific discharge should be able to deal with
+these roundoff errors.
+"""
+
+import flopy
+import numpy as np
+import pytest
+from flopy.mf6.utils import Mf6Splitter
+from framework import TestFramework
+
+single_case_name = "dis_single"
+split_case_name = "dis_split"
+cases = [single_case_name, split_case_name]
+
+spdis_lookup = {}
+
+model_name = "gwf_model"
+
+# solver criterion
+hclose = 1e-9
+max_inner_it = 50
+nper = 1
+
+# model spatial discretization
+nlay = 1
+nrow = 10
+ncol = 1
+
+# idomain
+idomain = np.ones((nlay, nrow, ncol))
+
+# cell size
+delr = 1.0
+delc = 1.0
+area = delr * delc
+
+# top/bot of the aquifer
+tops = [10.0, 0.0]
+
+# hydraulic conductivity
+hk = 10.0
+
+# boundary stress period data
+h_north = 5.0
+h_south = 0.0
+
+# initial head
+h_start = 0.0
+
+
+# head boundaries
+chd_spd = {0: [[(0, 0, 0), h_north], [(0, nrow - 1, 0), h_south]]}
+
+
+def get_model(idx, dir):
+    name = cases[idx]
+
+    # parameters and spd
+    # tdis
+    tdis_rc = []
+    for i in range(nper):
+        tdis_rc.append((1.0, 1, 1))
+
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name, version="mf6", exe_name="mf6", sim_ws=dir
+    )
+
+    tdis = flopy.mf6.ModflowTdis(sim, time_units="DAYS", nper=nper, perioddata=tdis_rc)
+
+    ims = flopy.mf6.ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        outer_dvclose=hclose,
+        outer_maximum=100,
+        under_relaxation="NONE",
+        inner_maximum=50,
+        inner_dvclose=0.1 * hclose,
+        rcloserecord=0.001,
+        linear_acceleration="CG",
+        scaling_method="NONE",
+        reordering_method="NONE",
+        relaxation_factor=0.0,
+        filename="gwf.ims",
+    )
+
+    gwf = flopy.mf6.ModflowGwf(sim, modelname=model_name, save_flows=True)
+
+    dis = flopy.mf6.ModflowGwfdis(
+        gwf,
+        nlay=nlay,
+        nrow=nrow,
+        ncol=ncol,
+        delr=delr,
+        delc=delc,
+        top=tops[0],
+        botm=tops[1:],
+        idomain=idomain,
+    )
+
+    # initial conditions
+    ic = flopy.mf6.ModflowGwfic(gwf, strt=h_start)
+
+    # node property flow
+    npf = flopy.mf6.ModflowGwfnpf(
+        gwf,
+        save_specific_discharge=True,
+        icelltype=0,
+        k=hk,
+    )
+
+    # chd file
+    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_spd)
+
+    # output control
+    oc = flopy.mf6.ModflowGwfoc(
+        gwf,
+        head_filerecord=f"{model_name}.hds",
+        budget_filerecord=f"{model_name}.cbc",
+        headprintrecord=[("COLUMNS", 10, "WIDTH", 15, "DIGITS", 6, "GENERAL")],
+        saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
+    )
+
+    if cases[idx] == split_case_name:
+        # split
+        splitter = Mf6Splitter(sim)
+        mask = np.zeros(shape=(nrow, ncol))
+        mask[nrow // 2 :, :] = 1
+        split_sim = splitter.split_model(mask)
+        split_sim.set_sim_path(dir)
+
+        # gwfgwf = split_sim.get_package("gwfgwf")
+        # gwfgwf.dev_interfacemodel_on = True
+        sim = split_sim
+
+    return sim
+
+
+def build_models(idx, test):
+    sim = get_model(idx, test.workspace)
+    return sim, None
+
+
+def check_output(idx, test):
+    qx = None
+    qy = None
+    qz = None
+    if cases[idx] == single_case_name:
+        sim = test.sims[0]
+        gwf = sim.get_model(model_name)
+        bud = gwf.output.budget()
+        spdis = bud.get_data(text="DATA-SPDIS")[0]
+        qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
+
+        # flatten
+        qx = qx[0, :, 0]
+        qy = qy[0, :, 0]
+        qz = qz[0, :, 0]
+
+    elif cases[idx] == split_case_name:
+        sim = test.sims[0]
+        gwf0 = sim.get_model(sim.model_names[0])
+        gwf1 = sim.get_model(sim.model_names[1])
+
+        bud0 = gwf0.output.budget()
+        spdis0 = bud0.get_data(text="DATA-SPDIS")[0]
+        qx0, qy0, qz0 = flopy.utils.postprocessing.get_specific_discharge(spdis0, gwf0)
+
+        bud1 = gwf1.output.budget()
+        spdis1 = bud1.get_data(text="DATA-SPDIS")[0]
+        qx1, qy1, qz1 = flopy.utils.postprocessing.get_specific_discharge(spdis1, gwf1)
+
+        # flatten and join
+        qx = np.concatenate((qx0[0, :, 0], qx1[0, :, 0]))
+        qy = np.concatenate((qy0[0, :, 0], qy1[0, :, 0]))
+        qz = np.concatenate((qz0[0, :, 0], qz1[0, :, 0]))
+
+    qx = qx[1:-1]
+    qy = qy[1:-1]
+    qz = qz[1:-1]
+    spdis_lookup[cases[idx]] = (qx, qy, qz)
+
+    qy_theory = -hk * (h_north - h_south) / ((nrow - 1) * delc)
+    assert np.allclose(qx, 0.0), "spdis cannot have x component in this problem"
+    assert np.allclose(qy, qy_theory), "spdis y component should equal theory"
+    assert np.allclose(qz, 0.0), "spdis cannot have z component in this problem"
+
+    if cases[idx] == split_case_name:
+        assert np.allclose(qx, spdis_lookup[single_case_name][0]), (
+            "spdis-x should match single model case"
+        )
+        assert np.allclose(qy, spdis_lookup[single_case_name][1]), (
+            "spdis-y should match single model case"
+        )
+        assert np.allclose(qz, spdis_lookup[single_case_name][2]), (
+            "spdis-z should match single model case"
+        )
+
+
+@pytest.mark.parametrize("idx, name", enumerate(cases))
+@pytest.mark.developmode
+def test_mf6model(idx, name, function_tmpdir, targets):
+    test = TestFramework(
+        name=name,
+        workspace=function_tmpdir,
+        build=lambda t: build_models(idx, t),
+        check=lambda t: check_output(idx, t),
+        targets=targets,
+    )
+    test.run()

autotest/test_gwf_uzf_gwet_opt.py

@@ -231,7 +231,7 @@ def set_connectiondata(n, lay, top, back, left, right, front, bottom):
         cl12.append(delc / 2)  # half the cell width along a column
         hwva.append(delr)  # the width perpendicular to the connection
         # for horizontal connection, 90.0 deg points in the positive y direction
-        angldeg.append(2700.0)
+        angldeg.append(270.0)
 
     if bottom:
         jas.append(n + (nrow * ncol))  # below

msvs/mf6core.vfproj

@@ -355,6 +355,7 @@
 		<File RelativePath="..\src\Model\ModelUtilities\PackageMover.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\SfrCrossSectionManager.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\SfrCrossSectionUtils.f90"/>
+		<File RelativePath="..\src\Model\ModelUtilities\SpdisCell.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\SpdisWorkArray.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\SwfCxsUtils.f90"/>
 		<File RelativePath="..\src\Model\ModelUtilities\TimeSelect.f90"/>

src/Exchange/exg-gwfgwf.f90

@@ -17,7 +17,8 @@ module GwfGwfExchangeModule
   use BaseExchangeModule, only: BaseExchangeType, AddBaseExchangeToList
   use BaseDisModule, only: DisBaseType
   use ConstantsModule, only: LENBOUNDNAME, NAMEDBOUNDFLAG, LINELENGTH, &
-                             TABCENTER, TABLEFT, LENAUXNAME, DNODATA
+                             TABCENTER, TABLEFT, LENAUXNAME, DNODATA, &
+                             DONE
   use ListModule, only: ListType
   use ListsModule, only: basemodellist
   use DisConnExchangeModule, only: DisConnExchangeType
@@ -761,7 +762,7 @@ subroutine gwf_gwf_set_flow_to_npf(this)
     real(DP) :: botn1, botn2
     real(DP) :: satn1, satn2
     real(DP) :: hn1, hn2
-    real(DP) :: nx, ny
+    real(DP) :: nx, ny, nz
     real(DP) :: distance
     real(DP) :: dltot
     real(DP) :: hwva
@@ -798,17 +799,21 @@ subroutine gwf_gwf_set_flow_to_npf(this)
       if (ihc == 0) then
         nx = DZERO
         ny = DZERO
+        nz = -DONE
+
         area = hwva
         if (botn1 < botn2) then
           ! -- n1 is beneath n2, so rate is positive downward.  Flip rate
           !    upward so that points in positive z direction
           rrate = -rrate
+          nz = DONE
         end if
       else
         if (this%ianglex > 0) then
           angle = this%auxvar(this%ianglex, i) * DPIO180
           nx = cos(angle)
           ny = sin(angle)
+          nz = DZERO
         else
           ! error?
           call store_error('error in gwf_gwf_cq', terminate=.TRUE.)
@@ -831,7 +836,7 @@ subroutine gwf_gwf_set_flow_to_npf(this)
       distance = dltot * this%cl1(i) / (this%cl1(i) + this%cl2(i))
       if (this%gwfmodel1%npf%icalcspdis == 1) then
         call this%gwfmodel1%npf%set_edge_properties(n1, ihc, rrate, area, &
-                                                    nx, ny, distance)
+                                                    nx, ny, nz, distance)
       end if
       !
       ! -- Submit this connection and flow information to the npf
@@ -845,7 +850,7 @@ subroutine gwf_gwf_set_flow_to_npf(this)
         distance = dltot * this%cl2(i) / (this%cl1(i) + this%cl2(i))
         if (ihc /= 0) rrate = -rrate
         call this%gwfmodel2%npf%set_edge_properties(n2, ihc, rrate, area, &
-                                                    -nx, -ny, distance)
+                                                    -nx, -ny, -nz, distance)
       end if
       !
     end do

src/Model/Connection/GwfGwfConnection.f90

@@ -667,7 +667,7 @@ subroutine setNpfEdgeProps(this)
           end if
           dist = conLen * cl / (imCon%cl1(isym) + imCon%cl2(isym))
           call this%gwfModel%npf%set_edge_properties(nLoc, ihc, rrate, area, &
-                                                     nx, ny, dist)
+                                                     nx, ny, nz, dist)
         end if
       end do
     end do