Skip to content

Commit

Permalink
fix(ModflowSfr2.load): don't try to read datasets 6b,c in cases where…
Browse files Browse the repository at this point in the history
… they're missing (modflowpy#671)
  • Loading branch information
aleaf committed Oct 17, 2019
1 parent 9ea47a4 commit 6e59d43
Show file tree
Hide file tree
Showing 2 changed files with 50 additions and 7 deletions.
36 changes: 36 additions & 0 deletions autotest/t009_test.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
import os
import glob
import shutil
import io
import numpy as np
from flopy.utils.recarray_utils import create_empty_recarray

Expand Down Expand Up @@ -49,6 +50,7 @@
'sfrfile': 'TL2009.sfr'}
}


def create_sfr_data():
dtype = np.dtype([('k', int),
('i', int),
Expand Down Expand Up @@ -82,6 +84,7 @@ def create_sfr_data():
d['outseg'] = [4, 0, 6, 8, 3, 8, 1, 2, 8]
return r, d


def sfr_process(mfnam, sfrfile, model_ws, outfolder=outpath):
m = flopy.modflow.Modflow.load(mfnam, model_ws=model_ws, verbose=False)
sfr = m.get_package('SFR')
Expand Down Expand Up @@ -252,6 +255,7 @@ def isequal(v1, v2):
assert isequal(sfr.reach_data.slope[21], 0.2)
assert isequal(sfr.reach_data.slope[-1], default_slope)


def test_const():

fm = flopy.modflow
Expand All @@ -273,6 +277,7 @@ def test_const():
assert sfr.const == 1.486 * 86400.
assert True


def test_export():
fm = flopy.modflow
m = fm.Modflow()
Expand Down Expand Up @@ -310,6 +315,7 @@ def test_export():
assert ra[(ra.iseg0 == 2) & (ra.ireach0 == 1)]['geometry'][0].bounds \
== (6.0, 2.0, 7.0, 3.0)


def test_example():
m = flopy.modflow.Modflow.load('test1ss.nam', version='mf2005',
exe_name='mf2005',
Expand Down Expand Up @@ -383,6 +389,35 @@ def test_example():
for i in range(1, nper):
assert sfr2.dataset_5[i][0] == -1


def test_no_ds_6bc():
"""Test case where datasets 6b and 6c aren't read
(e.g., see table at https://water.usgs.gov/ogw/modflow-nwt/MODFLOW-NWT-Guide/sfr.htm)
"""
sfrfiletxt = (
'REACHINPUT\n'
'2 2 0 0 128390 0.0001 119 0 3 10 1 30 0 4 0.75 91.54\n'
'1 1 1 1 1 1.0 1.0 0.001 1 1 .3 0.02 3.5 0.7\n'
'1 2 2 2 1 1.0 0.5 0.001 1 1 .3 0.02 3.5 0.7\n'
'2 2 0\n'
'1 2 2 0 0 0 0 0 0.041 0.111\n'
'0 3.64 7.28 10.92 14.55 18.19 21.83 25.47\n'
'2.55 1.02 0.76 0 0.25 0.76 1.02 2.55\n'
'2 2 0 0 0 0 0 0 0.041 0.111\n'
'0 3.96 7.92 11.88 15.83 19.79 23.75 27.71\n'
'2.77 1.11 0.83 0 0.28 0.83 1.11 2.77\n'
)
sfrfile = io.StringIO(sfrfiletxt)
m = fm.Modflow('junk', model_ws='temp')
sfr = fm.ModflowSfr2.load(sfrfile, model=m)
assert len(sfr.segment_data[0]) == 2
assert len(sfr.channel_geometry_data[0]) == 2
assert len(sfr.channel_geometry_data[0][1]) == 2
for i in range(2):
assert len(sfr.channel_geometry_data[0][1][i]) == 8
assert sum(sfr.channel_geometry_data[0][1][i]) > 0.


def test_ds_6d_6e_disordered():
path = os.path.join("..", "examples", "data", "hydmod_test")
path2 = os.path.join(".", "temp", "t009")
Expand Down Expand Up @@ -453,6 +488,7 @@ def test_transient_example():
assert m2.sfr.istcb2 == -49
assert m2.get_output_attribute(unit=abs(m2.sfr.istcb2), attr='binflag')


def test_assign_layers():
m = fm.Modflow()
m.dis = fm.ModflowDis(nrow=1, ncol=6, nlay=7,
Expand Down
21 changes: 14 additions & 7 deletions flopy/modflow/mfsfr2.py
Original file line number Diff line number Diff line change
Expand Up @@ -848,13 +848,20 @@ def load(f, model, nper=None, gwt=False, nsol=1, ext_unit_dict=None):
icalc = dataset_6a[1]
# link dataset 6d, 6e by nseg of dataset_6a
temp_nseg = dataset_6a[0]
dataset_6b = _parse_6bc(f.readline(), icalc, nstrm,
isfropt,
reachinput, per=i)
dataset_6c = _parse_6bc(f.readline(), icalc, nstrm,
isfropt,
reachinput, per=i)

# datasets 6b and 6c aren't read under the conditions below
# see table under description of dataset 6c,
# in the MODFLOW Online Guide for a description
# of this logic
# https://water.usgs.gov/ogw/modflow-nwt/MODFLOW-NWT-Guide/sfr.htm
dataset_6b, dataset_6c = (0,) * 9, (0,) * 9
if not (isfropt in [2, 3] and icalc == 1 and i > 1) and \
not (isfropt in [1, 2, 3] and icalc >= 2):
dataset_6b = _parse_6bc(f.readline(), icalc, nstrm,
isfropt,
reachinput, per=i)
dataset_6c = _parse_6bc(f.readline(), icalc, nstrm,
isfropt,
reachinput, per=i)
current[j] = dataset_6a + dataset_6b + dataset_6c

if icalc == 2:
Expand Down

0 comments on commit 6e59d43

Please sign in to comment.