fix(exg): trap multiple connections per cell pair

autotest/test_gwf_split_faces.py

@@ -0,0 +1,247 @@
+"""
+Test DISV grids where adjacent cells share a face with an extra vertex
+splitting their shared boundary. This situation is arranged and tested
+within a single grid, and in separate grids coupled by an exchange. In
+the exchange scenario, 2 connections are configured on the split face.
+
+Three test cases are included:
+1. GWF-GWF exchange without interface model
+2. GWF-GWF exchange with interface model enabled
+3. Single DISV grid with 2 adjacent cells (no exchange)
+
+All three cases should run successfully. The exchange cases should issue
+a warning about multiple connections between the same cell pair. The
+single-grid case should run silently.
+"""
+
+import flopy
+import pytest
+from framework import TestFramework
+
+cases = [
+    "gwfgwf_sf",
+    "gwfgwf_sf_ifmod",
+    "gwf_sf",
+]
+ifmod = [False, True, None]
+
+
+def build_models(idx, test):
+    nlay = 1
+    delr = 1.0
+    delc = 1.0
+    top = 1.0
+    botm = 0.0
+    hk = 1.0
+    name = cases[idx]
+    sim = flopy.mf6.MFSimulation(
+        sim_name=name,
+        version="mf6",
+        exe_name="mf6",
+        sim_ws=test.workspace,
+    )
+    tdis = flopy.mf6.ModflowTdis(
+        sim, time_units="DAYS", nper=1, perioddata=[(1.0, 1, 1.0)]
+    )
+    ims = flopy.mf6.ModflowIms(
+        sim,
+        print_option="SUMMARY",
+        complexity="SIMPLE",
+        outer_dvclose=1.0e-5,
+        outer_maximum=100,
+        under_relaxation="NONE",
+        inner_maximum=100,
+        inner_dvclose=1.0e-6,
+        rcloserecord=0.1,
+        linear_acceleration="BICGSTAB",
+        scaling_method="NONE",
+        reordering_method="NONE",
+        relaxation_factor=0.99,
+    )
+
+    # single grid case with 2 adjacent cells sharing a face with an extra vertex
+    if idx == 2:
+        gwf = flopy.mf6.ModflowGwf(sim, modelname="gwf", save_flows=True)
+        disv = flopy.mf6.ModflowGwfdisv(
+            gwf,
+            nlay=nlay,
+            ncpl=2,
+            nvert=7,
+            top=top,
+            botm=botm,
+            vertices=[
+                [0, 0.0, 0.0],  # v0: left cell bottom-left
+                [1, 1.0, 0.0],  # v1: shared bottom
+                [2, 1.0, 0.5],  # v2: shared middle (extra vertex)
+                [3, 1.0, 1.0],  # v3: shared top
+                [4, 0.0, 1.0],  # v4: left cell top-left
+                [5, 2.0, 0.0],  # v5: right cell bottom-right
+                [6, 2.0, 1.0],  # v6: right cell top-right
+            ],
+            cell2d=[
+                [0, 0.5, 0.5, 5, 0, 4, 3, 2, 1],  # left
+                [1, 1.5, 0.5, 5, 1, 2, 3, 6, 5],  # right
+            ],
+        )
+        ic = flopy.mf6.ModflowGwfic(gwf, strt=1.0)
+        npf = flopy.mf6.ModflowGwfnpf(
+            gwf,
+            save_flows=True,
+            icelltype=0,
+            k=hk,
+        )
+        chd = flopy.mf6.ModflowGwfchd(
+            gwf, stress_period_data=[[(0, 0), 1.0], [(0, 1), 0.0]]
+        )
+        oc = flopy.mf6.ModflowGwfoc(
+            gwf,
+            head_filerecord="gwf.hds",
+            budget_filerecord="gwf.cbc",
+            saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
+        )
+
+        return sim, None
+
+    # left model
+    gwf1 = flopy.mf6.ModflowGwf(sim, modelname="gwf1", save_flows=True)
+    vertices1 = [
+        [0, 0.0, 0.0],  # v0: bottom-left
+        [1, 1.0, 0.0],  # v1: bottom-right
+        [2, 1.0, 0.5],  # v2: mid-right (extra vertex)
+        [3, 1.0, 1.0],  # v3: top-right
+        [4, 0.0, 1.0],  # v4: top-left
+    ]
+    disv1 = flopy.mf6.ModflowGwfdisv(
+        gwf1,
+        nlay=nlay,
+        ncpl=1,
+        nvert=5,
+        top=top,
+        botm=botm,
+        vertices=vertices1,
+        # [icell2d, xc, yc, nvert, v0, v1, v2, ...]
+        cell2d=[[0, 0.5, 0.5, 5, 0, 4, 3, 2, 1]],
+    )
+    ic1 = flopy.mf6.ModflowGwfic(gwf1, strt=1.0)
+    npf1 = flopy.mf6.ModflowGwfnpf(
+        gwf1,
+        save_flows=True,
+        icelltype=0,
+        k=hk,
+    )
+    chd1 = flopy.mf6.ModflowGwfchd(gwf1, stress_period_data=[[(0, 0), 1.0]])
+    oc1 = flopy.mf6.ModflowGwfoc(
+        gwf1,
+        head_filerecord="gwf1.hds",
+        budget_filerecord="gwf1.cbc",
+        saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
+    )
+
+    # right model
+    gwf2 = flopy.mf6.ModflowGwf(sim, modelname="gwf2", save_flows=True)
+    disv2 = flopy.mf6.ModflowGwfdisv(
+        gwf2,
+        nlay=nlay,
+        ncpl=1,
+        nvert=5,
+        top=top,
+        botm=botm,
+        vertices=[
+            [0, 1.0, 0.0],  # v0: bottom-left
+            [1, 2.0, 0.0],  # v1: bottom-right
+            [2, 2.0, 1.0],  # v2: top-right
+            [3, 1.0, 1.0],  # v3: top-left
+            [4, 1.0, 0.5],  # v4: mid-left (extra vertex)
+        ],
+        # [icell2d, xc, yc, nvert, v0, v1, v2, ...]
+        cell2d=[[0, 1.5, 0.5, 5, 0, 4, 3, 2, 1]],
+    )
+    ic2 = flopy.mf6.ModflowGwfic(gwf2, strt=0.0)
+    npf2 = flopy.mf6.ModflowGwfnpf(
+        gwf2,
+        save_flows=True,
+        icelltype=0,
+        k=hk,
+    )
+    chd2 = flopy.mf6.ModflowGwfchd(gwf2, stress_period_data=[[(0, 0), 0.0]])
+    oc2 = flopy.mf6.ModflowGwfoc(
+        gwf2,
+        head_filerecord="gwf2.hds",
+        budget_filerecord="gwf2.cbc",
+        saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")],
+    )
+
+    # create exchange with 2 connections, 1 for each segment of the split face
+    angldegx = 0.0
+    cdist = delr
+    gwfgwf_data = [
+        [
+            (0, 0),  # cell in model 1 (layer, cellid)
+            (0, 0),  # cell in model 2 (layer, cellid)
+            1,  # ihc (horizontal connection)
+            0.5,  # cl1 (distance from gwf1 cell center to face)
+            0.5,  # cl2 (distance from gwf2 cell center to face)
+            delc / 2.0,  # hwva (flow width = half cell height)
+            angldegx,  # auxiliary variable
+            cdist,  # auxiliary variable
+        ],
+        [
+            (0, 0),  # cell in model 1 (layer, cellid)
+            (0, 0),  # cell in model 2 (layer, cellid)
+            1,  # ihc (horizontal connection)
+            0.5,  # cl1 (distance from gwf1 cell center to face)
+            0.5,  # cl2 (distance from gwf2 cell center to face)
+            delc / 2.0,  # hwva (flow width = half cell height)
+            angldegx,  # auxiliary variable
+            cdist,  # auxiliary variable
+        ],
+    ]
+    gwfgwf = flopy.mf6.ModflowGwfgwf(
+        sim,
+        exgtype="GWF6-GWF6",
+        save_flows=True,
+        print_flows=True,
+        nexg=len(gwfgwf_data),
+        exgmnamea="gwf1",
+        exgmnameb="gwf2",
+        exchangedata=gwfgwf_data,
+        auxiliary=["ANGLDEGX", "CDIST"],
+        dev_interfacemodel_on=ifmod[idx],
+    )
+
+    return sim, None
+
+
+def check_output(idx, test):
+    if idx == 2:
+        # single grid case
+        fpth = test.workspace / "gwf.hds"
+        hds = flopy.utils.HeadFile(fpth)
+        heads = hds.get_data()
+
+        head1 = heads[0, 0, 0]
+        head2 = heads[0, 0, 1]
+
+        assert 0.9 < head1 < 1.1
+        assert -0.1 < head2 < 0.1
+        assert head1 > head2
+    else:
+        # exchange cases: expect a warning about multiple connections
+        with open(test.workspace / "mfsim.lst") as f:
+            lines = f.readlines()
+        assert any(
+            "has multiple connections between cells" in line for line in lines
+        ), "expected warning about multiple connections not found in mfsim.lst"
+
+
+@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),
+        xfail=False,
+    )
+    test.run()

src/Exchange/exg-gwegwe.f90

@@ -10,9 +10,9 @@
 module GweGweExchangeModule
 
   use KindModule, only: DP, I4B, LGP
-  use SimVariablesModule, only: errmsg, model_loc_idx
+  use SimVariablesModule, only: errmsg, warnmsg, model_loc_idx
   use SimModule, only: store_error, store_error_filename, &
-                       count_errors, ustop
+                       count_errors, ustop, store_warning
   use BaseModelModule, only: BaseModelType, GetBaseModelFromList
   use BaseExchangeModule, only: BaseExchangeType, AddBaseExchangeToList
   use ConstantsModule, only: LENBOUNDNAME, NAMEDBOUNDFLAG, LINELENGTH, &
@@ -249,8 +249,13 @@ end subroutine gwe_gwe_df
   subroutine validate_exchange(this)
     ! -- dummy
     class(GweExchangeType) :: this !<  GweExchangeType
+    ! -- local
+    integer(I4B) :: n, m, ncount
+    logical(LGP), dimension(:), allocatable :: checked_conns
+    !
+    ! In parallel, only validate on the owning process
+    if (this%is_datacopy) return
     !
-
     ! Ensure gwfmodel names were entered
     if (this%gwfmodelname1 == '') then
       write (errmsg, '(3a)') 'GWE-GWE exchange ', trim(this%name), &
@@ -297,6 +302,31 @@ subroutine validate_exchange(this)
       call store_error(errmsg)
     end if
     !
+    ! At most one connection is allowed between any two cells.
+    allocate (checked_conns(this%nexg))
+    checked_conns = .false.
+    do n = 1, this%nexg
+      if (checked_conns(n)) cycle
+      ! Skip connections with unresolved nodes (model on remote process)
+      if (this%nodem1(n) == -1 .or. this%nodem2(n) == -1) cycle
+      ncount = 0
+      do m = 1, this%nexg
+        if (this%nodem1(m) == this%nodem1(n) .and. &
+            this%nodem2(m) == this%nodem2(n)) then
+          ncount = ncount + 1
+          checked_conns(m) = .true.
+        end if
+      end do
+      if (ncount > 1) then
+        write (warnmsg, '(a,a,a,i0,a,i0,a,i0,a)') &
+          'GWE-GWE exchange ', trim(this%name), &
+          ' has multiple connections between cells ', this%nodem1(n), &
+          ' (model 1) and ', this%nodem2(n), ' (model 2): ', ncount, &
+          ' connections defined for this cell pair.'
+        call store_warning(warnmsg)
+      end if
+    end do
+    !
     if (count_errors() > 0) then
       call ustop()
     end if

src/Exchange/exg-gwfgwf.f90

@@ -10,9 +10,9 @@
 module GwfGwfExchangeModule
 
   use KindModule, only: DP, I4B, LGP
-  use SimVariablesModule, only: errmsg
+  use SimVariablesModule, only: errmsg, warnmsg
   use SimModule, only: count_errors, store_error, store_error_filename, &
-                       store_error_unit
+                       store_error_unit, store_warning
   use BaseModelModule, only: BaseModelType, GetBaseModelFromList
   use BaseExchangeModule, only: BaseExchangeType, AddBaseExchangeToList
   use BaseDisModule, only: DisBaseType
@@ -30,6 +30,7 @@ module GwfGwfExchangeModule
   use SimVariablesModule, only: errmsg, model_loc_idx
   use TableModule, only: TableType, table_cr
   use MatrixBaseModule
+  use DisvGeom, only: DisvGeomType
 
   implicit none
 
@@ -275,6 +276,11 @@ subroutine validate_exchange(this)
     class(GwfExchangeType) :: this !<  GwfExchangeType
     ! -- local
     logical(LGP) :: has_k22, has_spdis, has_vsc
+    integer(I4B) :: n, m, ncount
+    logical(LGP), dimension(:), allocatable :: checked_conns
+    !
+    ! In parallel, only validate on the owning process
+    if (this%is_datacopy) return
     !
     ! Periodic boundary condition in exchange don't allow XT3D (=interface model)
     if (this%v_model1 == this%v_model2) then
@@ -354,6 +360,31 @@ subroutine validate_exchange(this)
         ' in both of the connected models.'
       call store_error(errmsg, terminate=.TRUE.)
     end if
+    !
+    ! At most one connection is allowed between any two cells.
+    allocate (checked_conns(this%nexg))
+    checked_conns = .false.
+    do n = 1, this%nexg
+      if (checked_conns(n)) cycle
+      ! Skip connections with unresolved nodes (model on remote process)
+      if (this%nodem1(n) == -1 .or. this%nodem2(n) == -1) cycle
+      ncount = 0
+      do m = 1, this%nexg
+        if (this%nodem1(m) == this%nodem1(n) .and. &
+            this%nodem2(m) == this%nodem2(n)) then
+          ncount = ncount + 1
+          checked_conns(m) = .true.
+        end if
+      end do
+      if (ncount > 1) then
+        write (warnmsg, '(a,a,a,i0,a,i0,a,i0,a)') &
+          'GWF-GWF exchange ', trim(this%name), &
+          ' has multiple connections between cells ', this%nodem1(n), &
+          ' (model 1) and ', this%nodem2(n), ' (model 2): ', ncount, &
+          ' connections defined for this cell pair.'
+        call store_warning(warnmsg)
+      end if
+    end do
   end subroutine validate_exchange
 
   !> @ brief Add connections