Extract water budget using binaryfile

[hjia1005]

I want to use binaryfile to extract budget using the following code. I have two GHB boundaries: GHB2 and GHB3, but rec only extracts GHB2. I tried to change text to 'GHB2' or 'GHB3' then it raised the 'Exception: The specified text string is not in the budget file.'

import flopy.utils.binaryfile as bf
cbb = bf.CellBudgetFile('mymodel.cbb')
cbb.list_records()
rec = cbb.get_data(kstpkper=(0,0), text='GHB')

When I list the records, I can see 'GHB2' and 'GHB3'. Also if I use zbud6.exe it does separate these two. How can I extract and separate them using bf?

Thanks.

[langevin-usgs]

It's true that the multi-package support from CellBudgetFile can be tricky. But the get_data() method with text should give you back all the records for a specific package type.

This example is adapted from the front page of the flopy site to use two constant head packages.

import flopy
ws = './temp'
name = 'temp'
sim = flopy.mf6.MFSimulation(sim_name=name, sim_ws=ws, exe_name='mf6')
tdis = flopy.mf6.ModflowTdis(sim)
ims = flopy.mf6.ModflowIms(sim)
gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)
dis = flopy.mf6.ModflowGwfdis(gwf, nrow=10, ncol=10)
ic = flopy.mf6.ModflowGwfic(gwf)
npf = flopy.mf6.ModflowGwfnpf(gwf, save_specific_discharge=True)
chd1 = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 0, 0), 1.]], pname="chd1")
chd2 = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=[[(0, 9, 9), 0.]], pname="chd2")
budget_file = name + '.bud'
head_file = name + '.hds'
oc = flopy.mf6.ModflowGwfoc(gwf,
                            budget_filerecord=budget_file,
                            head_filerecord=head_file,
                            saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')])
sim.write_simulation()
sim.run_simulation(silent=True)

bud = gwf.output.budget()
bud.get_data(text="chd")

The get_data(text="chd") method returns returns the following, which you can see is the constant head flows for package chd1 and package chd2 provided as two separate entries in the list.

[rec.array([(1, 1, 0.33205402)],
           dtype=[('node', '<i4'), ('node2', '<i4'), ('q', '<f8')]),
 rec.array([(100, 1, -0.33204178)],
           dtype=[('node', '<i4'), ('node2', '<i4'), ('q', '<f8')])]

I verified that this also works for multiple ghb packages.

[hjia1005]

Yes it does get two entries in the list, I just used the first one. Thank you very much!