bug: issues with get_ts support for aux vars

[wpbonelli]

Hey @wpbonelli ,
I've been testing this functionality and seem to run into two issues (I'm using flopy version 3.10.0.dev5 on Windows):

• For the use-case attached below, zeros are returned for a time series of qx whereas these values are non-zero when obtaining them from the get_data() method.
• For DISV grids, the bud.get_ts() method fails on the following line. This is also true when text points to a different record, e.g. 'RCHA').
• Also, could the idx argument in the docstrings of get_ts() also include the required format for DISV grids, i.e. (layer, dummy, cellid2d)

Example code

import numpy as np
import os
import flopy
from flopy.utils.triangle import Triangle

gridtype = 'dis' # disv or dis for triangle or rectilinear grid


# create MODFLOW 6 model
sim = flopy.mf6.MFSimulation(sim_name='evt-test', sim_ws='model')
tdis = flopy.mf6.ModflowTdis(sim, nper=1, perioddata=[(100.0, 10, 1.0)])
ims = flopy.mf6.ModflowIms(sim)

gwf = flopy.mf6.ModflowGwf(sim, modelname='evt-test', save_flows=True)

if gridtype == 'disv':
  # create triangular grid
  triangle_ws = 'triangle'
  if not os.path.exists(triangle_ws):
      os.mkdir(triangle_ws)
      
  active_area = [(0, 0), (0, 1000), (1000, 1000), (1000, 0)]
  tri = Triangle(model_ws=triangle_ws, angle=30)
  tri.add_polygon(active_area)
  tri.add_region((1, 1), maximum_area=50**2)

  tri.build()

  # build vertex grid object
  vgrid = flopy.discretization.VertexGrid(
      vertices=tri.get_vertices(),
      cell2d=tri.get_cell2d(),
  )


  disv = flopy.mf6.ModflowGwfdisv(
      gwf,
      nlay=1,
      top=0.0,
      botm=-10.0,
      ncpl=vgrid.ncpl,
      nvert=vgrid.nvert,
      cell2d=vgrid.cell2d,
      vertices=tri.get_vertices() # this is not stored in the Vertex grid object?
  )


elif gridtype == 'dis':
  dis = flopy.mf6.ModflowGwfdis(gwf, nlay=1, nrow=20, ncol=20, delr=50, delc=50, top=0.0, botm=-10)
else:
  raise ValueError('gridtype should be "dis" or "disv"')

ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0)

grid = gwf.modelgrid
chd_id0 = grid.intersect(0.1, 0.1)
chd_id1 = grid.intersect(999.9, 999.9)
pb = grid.intersect(800, 800)

if gridtype == 'disv':
  pbidx_full = (0, 0, pb) # requires a dummy index in the second position
  pbidx = (0, pb)
  chdspd = [
  [(0, chd_id0), 1.0], 
  [(0, chd_id1), 0.0],
  ]
else:
  pbidx = (0,) + pb
  pbidx_full = pbidx
  chdspd = [
  [(0,) + chd_id0, 1.0], 
  [(0,) + chd_id1, 0.0],
  ]

chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chdspd)
rcha = flopy.mf6.ModflowGwfrcha(gwf, recharge=0.4 / 365)
# evta = flopy.mf6.ModflowGwfevta(gwf, rate=0.2 / 365, surface=0.0, depth=1)

npf = flopy.mf6.ModflowGwfnpf(gwf, k=10.0, save_specific_discharge=True, icelltype=1)
sto = flopy.mf6.ModflowGwfsto(gwf, sy=0.15, ss=1e-5, transient={0: True})
oc = flopy.mf6.ModflowGwfoc(gwf, 
                          head_filerecord=f'{gwf.name}.hds',
                          budget_filerecord=f'{gwf.name}.bud',
                          saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')],
                          printrecord=[('BUDGET', 'ALL')],
                          )

# run
sim.write_simulation(silent=True)
succes, pbuff = sim.run_simulation(silent=True)
assert succes, pbuff

# get head time series
hds = gwf.output.head()
h = hds.get_ts(idx=pbidx_full)
print(h)

# get qx time series
bud = gwf.output.budget()
qx = bud.get_ts(idx=pbidx_full, text='DATA-SPDIS', variable='qx') # fails with DISV
print(qx) # zeros

# check with spdis array
spdis = bud.get_data(text='DATA-SPDIS')[0]
qx, qy, qz = flopy.utils.postprocessing.get_specific_discharge(spdis, gwf)
print(qx[pbidx]) # non-zero, at the first time step

Originally posted by @cneyens in #2270