bug: GWF stress package boundnames prevent saving MFSplitter node mapping

[Isaac-AGS]

Describe the bug
I am attempting to use MFSplitter to develop a parallel groundwater flow model simulation in MODFLOW6. However, the node mapping procedure gets interrupted on the backend due to the use of boundnames in the stress packages.

To Reproduce
Below is an example model I created to reproduce the error. I use boundnames in the WEL package.

# Imports
import geopandas
import numpy as np
import flopy
from flopy.utils.gridintersect import GridIntersect
from flopy.mf6.utils import Mf6Splitter
from flopy.utils.geometry import Polygon

# Create base simulation
name, ws = 'model', '.'
sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws)
tdis = flopy.mf6.ModflowTdis(sim, time_units='days')
ims = flopy.mf6.ModflowIms(sim)

# Create base GWF model
gwf = flopy.mf6.ModflowGwf(sim, modelname=name)

# DIS package
delr, delc = 200, 200
nlay, nrow, ncol = 3, 30, 30
shapefile = geopandas.read_file('shapefile.shp')
centroid = shapefile.geometry.centroid
xll = centroid.x.iloc[0]-nrow*delc/2
yll = centroid.y.iloc[0]-nrow*delc/2

dis = flopy.mf6.ModflowGwfdis(
    gwf,
    length_units='feet',
    nlay=nlay, nrow=nrow, ncol=ncol,
    delr=delr, delc=delc,
    top=300, botm=[200, 100, 0],
    xorigin=xll, yorigin=yll, 
)

# Initial conditions
ic = flopy.mf6.ModflowGwfic(gwf, strt=220)

# NPF package
npf = flopy.mf6.ModflowGwfnpf(
    gwf,
    k=[50, 0.01, 200],
    k33=[10, 0.01, 20],
    icelltype=[1, 0, 0],
)

# RCH
rch = flopy.mf6.ModflowGwfrcha(gwf, recharge=1e-3)

# WEL
wel = flopy.mf6.ModflowGwfwel(
    gwf, 
    stress_period_data={0: [(2, 10, 9, -1e5, 'my_well')]},
    boundnames=True
)

# OC
oc = flopy.mf6.ModflowGwfoc(
    gwf,
    head_filerecord=f"{name}.hds",
    budget_filerecord=f"{name}.cbc",
    saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
)

# Write simulation and run
sim.write_simulation()
sim.run_simulation()

# Perform a structured grid intersect to create the split array
xy = shapefile.exterior.get_coordinates().to_numpy()
polygons= Polygon(xy)

modelgrid = gwf.modelgrid
ix = GridIntersect(modelgrid)
ix_poly = ix.intersect(polygons, contains_centroid=True)
split_nodes = set(ix_poly.cellids)

split_nodes = np.array(list(split_nodes)).T
split_array = np.zeros(modelgrid.shape[1:])
split_array[*split_nodes] = 1

# Split the simulation
mfsplit = Mf6Splitter(sim)
split_sim = mfsplit.split_model(split_array)

# Save split sim node mapping
mfsplit.save_node_mapping('mfsplit_node_mapping.hdf5')

Running this produces the following error message:

Traceback (most recent call last):

  File ~\miniconda3\envs\flopy_env\Lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\ags\_scripts\flopy\modflow6\mfsplitter\_bug_report\mfsplitter_bug_report.py:112
    mfsplit.save_node_mapping('mfsplit_node_mapping.hdf5')

  File ~\miniconda3\envs\flopy_env\Lib\site-packages\flopy\mf6\utils\model_splitter.py:400 in save_node_mapping
    node_map = {

  File ~\miniconda3\envs\flopy_env\Lib\site-packages\flopy\mf6\utils\model_splitter.py:401 in <dictcomp>
    int(k): (int(v[0]), int(v[1])) for k, v in self._node_map.items()

ValueError: invalid literal for int() with base 10: 'my_well'

Expected behavior
If used, stress package boundnames get appended to the end of the MFSplitter _node_map attribute. Without using boundnames the code runs without errors and the node map attribute contains of the expected node pairs.

Desktop:

• OS: Windows 11 Home, version 25H2
• IDE: Spyder 6
• FloPy Version: 3.10.0