From 45df13abfd0611cc423b01bff341b450064e3b01 Mon Sep 17 00:00:00 2001 From: jdhughes-usgs Date: Tue, 5 Oct 2021 11:29:13 -0500 Subject: [PATCH] feat(CellBudget): add support for full3D keyword (#1254) --- autotest/t004_test_utilarray.py | 18 +- autotest/t050_test.py | 35 ++- ..._thick.py => t076_test_modelgrid_thick.py} | 0 ...tions.py => t078_test_lake_connections.py} | 0 autotest/t079_test_cbc_full3D.py | 205 ++++++++++++++++++ autotest/t505_test.py | 50 +++-- .../data/mf6/test001a_Tharmonic/flow15.ims | 6 +- .../test001e_UZF_3lay/test001e_UZF_3lay.ims | 4 +- examples/data/mf6/test003_gwfs_disv/model.ims | 10 +- .../data/mf6/test003_gwftri_disv/model.ims | 10 +- .../data/mf6/test005_advgw_tidal/model.ims | 4 +- .../mf6/test006_2models_mvr/simulation.ims | 6 +- examples/data/mf6/test006_gwf3/flow.ims | 6 +- .../data/mf6/test027_TimeseriesTest/model.ims | 6 +- .../data/mf6/test036_twrihfb/twrihfb2015.ims | 6 +- .../mf6/test045_lake1ss_table/lakeex1b.ims | 4 +- .../data/mf6/test045_lake2tr/lakeex2a.ims | 12 +- flopy/utils/binaryfile.py | 85 +++++--- 18 files changed, 355 insertions(+), 112 deletions(-) rename autotest/{t076_modelgrid_thick.py => t076_test_modelgrid_thick.py} (100%) rename autotest/{t078_lake_connections.py => t078_test_lake_connections.py} (100%) create mode 100644 autotest/t079_test_cbc_full3D.py diff --git a/autotest/t004_test_utilarray.py b/autotest/t004_test_utilarray.py index 282a7bcd48..c776951d4e 100644 --- a/autotest/t004_test_utilarray.py +++ b/autotest/t004_test_utilarray.py @@ -17,7 +17,7 @@ def test_load_txt_free(): a = np.ones((10,), dtype=np.float32) * 250.0 - fp = StringIO(u"10*250.0") + fp = StringIO("10*250.0") fa = Util2d.load_txt(a.shape, fp, a.dtype, "(FREE)") np.testing.assert_equal(fa, a) assert fa.dtype == a.dtype @@ -25,7 +25,7 @@ def test_load_txt_free(): a = np.arange(10, dtype=np.int32).reshape((2, 5)) fp = StringIO( dedent( - u"""\ + """\ 0 1,2,3, 4 5 6, 7, 8 9 """ @@ -39,7 +39,7 @@ def test_load_txt_free(): a[1, 0] = 2.2 fp = StringIO( dedent( - u"""\ + """\ 5*1.0 2.2 2*1.0, +1E-00 1.0 """ @@ -54,7 +54,7 @@ def test_load_txt_fixed(): a = np.arange(10, dtype=np.int32).reshape((2, 5)) fp = StringIO( dedent( - u"""\ + """\ 01234X 56789 """ @@ -66,7 +66,7 @@ def test_load_txt_fixed(): fp = StringIO( dedent( - u"""\ + """\ 0123X 4 5678 @@ -81,7 +81,7 @@ def test_load_txt_fixed(): a = np.array([[-1, 1, -2, 2, -3], [3, -4, 4, -5, 5]], np.int32) fp = StringIO( dedent( - u"""\ + """\ -1 1-2 2-3 3 -44 -55 """ @@ -96,7 +96,7 @@ def test_load_block(): a = np.ones((2, 5), dtype=np.int32) * 4 fp = StringIO( dedent( - u"""\ + """\ 1 1 2 1 5 4 """ @@ -111,7 +111,7 @@ def test_load_block(): a[0, 2:4] = 6.0 fp = StringIO( dedent( - u"""\ + """\ 3 1 2 1 5 4.0 1 2 2 2 9.0 @@ -127,7 +127,7 @@ def test_load_block(): a[0, 2:4] = 8 fp = StringIO( dedent( - u"""\ + """\ 1 1 1 3 4 8 """ diff --git a/autotest/t050_test.py b/autotest/t050_test.py index 5064f3fd1e..8e5588950f 100644 --- a/autotest/t050_test.py +++ b/autotest/t050_test.py @@ -46,8 +46,9 @@ def test_vtk_export_array2d(): assert nlines == 17615 # with smoothing - m.dis.top.export(output_dir, fmt="vtk", name="top_smooth", - binary=False, smooth=True) + m.dis.top.export( + output_dir, fmt="vtk", name="top_smooth", binary=False, smooth=True + ) filetocheck = os.path.join(output_dir, "top_smooth.vtk") nlines1 = count_lines_in_file(filetocheck) assert nlines1 == 17615 @@ -79,7 +80,11 @@ def test_vtk_export_array3d(): # with point scalars m.upw.hk.export( - output_dir, fmt="vtk", name="hk_points", point_scalars=True, binary=False + output_dir, + fmt="vtk", + name="hk_points", + point_scalars=True, + binary=False, ) filetocheck = os.path.join(output_dir, "hk_points.vtk") nlines1 = count_lines_in_file(filetocheck) @@ -117,7 +122,9 @@ def test_vtk_transient_array_2d(): kpers = [0, 1, 1096] # export and check - m.rch.rech.export(output_dir, fmt="vtk", kpers=kpers, binary=False, xml=True) + m.rch.rech.export( + output_dir, fmt="vtk", kpers=kpers, binary=False, xml=True + ) filetocheck = os.path.join(output_dir, "rech_000001.vtk") nlines = count_lines_in_file(filetocheck) assert nlines == 26837 @@ -179,7 +186,9 @@ def test_vtk_export_packages(): # transient package drain kpers = [0, 1, 1096] output_dir = os.path.join(cpth, "DRN") - m.drn.export(output_dir, fmt="vtk", binary=False, xml=True, kpers=kpers, pvd=True) + m.drn.export( + output_dir, fmt="vtk", binary=False, xml=True, kpers=kpers, pvd=True + ) filetocheck = os.path.join(output_dir, "DRN_000001.vtu") nlines3 = count_lines_in_file(filetocheck) assert nlines3 == 27239 @@ -498,8 +507,9 @@ def test_vtk_vertex(): return # disv test - workspace = os.path.join("..", "examples", "data", "mf6", - "test003_gwfs_disv") + workspace = os.path.join( + "..", "examples", "data", "mf6", "test003_gwfs_disv" + ) # outfile = os.path.join("vtk_transient_test", "vtk_pacakages") outfile = os.path.join("temp", "t050", "vtk_disv", "disv.vtk") sim = flopy.mf6.MFSimulation.load(sim_ws=workspace) @@ -538,8 +548,9 @@ def test_vtk_pathline(): ws = os.path.join("..", "examples", "data", "freyberg") modelpth = os.path.join("temp", "t050") outfile = os.path.join(modelpth, "pathline_test", "pathline.vtk") - ml = flopy.modflow.Modflow.load("freyberg.nam", model_ws=ws, - exe_name="mf2005") + ml = flopy.modflow.Modflow.load( + "freyberg.nam", model_ws=ws, exe_name="mf2005" + ) ml.change_model_ws(new_pth=modelpth) ml.write_input() ml.run_model() @@ -587,8 +598,8 @@ def test_vtk_pathline(): maxtime = 0 for p in plines: - if np.max(p['time']) > maxtime: - maxtime = np.max(p['time']) + if np.max(p["time"]) > maxtime: + maxtime = np.max(p["time"]) if not len(totim) == 12054: raise AssertionError("Array size is incorrect for modpath VTK") @@ -613,4 +624,4 @@ def test_vtk_pathline(): test_vtk_vector() test_vtk_unstructured() test_vtk_vertex() - test_vtk_pathline() \ No newline at end of file + test_vtk_pathline() diff --git a/autotest/t076_modelgrid_thick.py b/autotest/t076_test_modelgrid_thick.py similarity index 100% rename from autotest/t076_modelgrid_thick.py rename to autotest/t076_test_modelgrid_thick.py diff --git a/autotest/t078_lake_connections.py b/autotest/t078_test_lake_connections.py similarity index 100% rename from autotest/t078_lake_connections.py rename to autotest/t078_test_lake_connections.py diff --git a/autotest/t079_test_cbc_full3D.py b/autotest/t079_test_cbc_full3D.py new file mode 100644 index 0000000000..e0025a483d --- /dev/null +++ b/autotest/t079_test_cbc_full3D.py @@ -0,0 +1,205 @@ +import os +import sys +import shutil + +import numpy as np +import flopy + +ex_pths = ( + os.path.join("..", "examples", "data", "freyberg"), + os.path.join("..", "examples", "data", "mf6-freyberg"), + os.path.join("..", "examples", "data", "mf6", "test003_gwfs_disv"), + os.path.join("..", "examples", "data", "mf6", "test003_gwftri_disv"), +) +ismf6_lst = ["mf6" in pth for pth in ex_pths] +names = [os.path.basename(pth) for pth in ex_pths] + +tpth = os.path.join("temp", "t079") +# delete the directory if it exists +if os.path.isdir(tpth): + shutil.rmtree(tpth) +# make the directory +os.makedirs(tpth) + +mf6_exe = "mf6" +mf2005_exe = "mf2005" +if sys.platform == "win32": + mf6_exe += ".exe" + mf2005_exe += ".exe" + + +def load_mf2005(name, ws_in): + name_file = f"{name}.nam" + ml = flopy.modflow.Modflow.load( + name_file, + model_ws=ws_in, + exe_name=mf2005_exe, + check=False, + ) + + # change work space + ws_out = os.path.join(tpth, name) + ml.change_model_ws(ws_out) + + # save all budget data to a cell-by cell file + oc = ml.get_package("OC") + oc.reset_budgetunit() + oc.stress_period_data = {(0, 0): ["save budget"]} + + return ml + + +def load_mf6(name, ws_in): + sim = flopy.mf6.MFSimulation.load( + sim_name=name, + exe_name=mf6_exe, + sim_ws=ws_in, + ) + + # change work space + ws_out = os.path.join(tpth, name) + sim.set_sim_path(ws_out) + + # get the groundwater flow model(s) and redefine the output control + # file to save cell by cell output + models = sim.model_names + for model in models: + gwf = sim.get_model(model) + gwf.name_file.save_flows = True + + gwf.remove_package("oc") + budget_filerecord = f"{model}.cbc" + oc = flopy.mf6.ModflowGwfoc( + gwf, + budget_filerecord=budget_filerecord, + saverecord=[("BUDGET", "ALL")], + ) + + return sim + + +def cbc_eval_size(cbcobj, nnodes, shape3d): + cbc_pth = cbcobj.filename + + assert cbcobj.nnodes == nnodes, ( + f"{cbc_pth} nnodes ({cbcobj.nnodes}) " f"does not equal {nnodes}" + ) + a = np.squeeze(np.ones(cbcobj.shape, dtype=float)) + b = np.squeeze(np.ones(shape3d, dtype=float)) + assert a.shape == b.shape, ( + f"{cbc_pth} shape {cbcobj.shape} " f"does not conform to {shape3d}" + ) + + +def cbc_eval_data(cbcobj, shape3d): + cbc_pth = cbcobj.filename + print(f"{cbc_pth}:\n") + cbcobj.list_unique_records() + + names = cbcobj.get_unique_record_names(decode=True) + times = cbcobj.get_times() + for name in names: + text = name.strip() + arr = np.squeeze( + cbcobj.get_data(text=text, totim=times[0], full3D=True)[0] + ) + if text != "FLOW-JA-FACE": + b = np.squeeze(np.ones(shape3d, dtype=float)) + assert arr.shape == b.shape, ( + f"{cbc_pth} shape {arr.shape} for '{text}' budget item " + f"does not conform to {shape3d}" + ) + + +def cbc_eval(cbcobj, nnodes, shape3d, modelgrid): + cbc_pth = cbcobj.filename + cbc_eval_size(cbcobj, nnodes, shape3d) + cbc_eval_data(cbcobj, shape3d) + cbcobj.close() + + cobj_mg = flopy.utils.CellBudgetFile( + cbc_pth, + modelgrid=modelgrid, + verbose=True, + ) + cbc_eval_size(cobj_mg, nnodes, shape3d) + cbc_eval_data(cobj_mg, shape3d) + cobj_mg.close() + + return + + +def clean_run(name): + ws = os.path.join(tpth, name) + if os.path.isdir(ws): + shutil.rmtree(ws) + + +def mf6_eval(name, ws_in): + + sim = load_mf6(name, ws_in) + + # write the simulation + sim.write_simulation() + + # run the simulation + sim.run_simulation() + + # get the groundwater model and determine the size of the model grid + gwf_name = list(sim.model_names)[0] + gwf = sim.get_model(gwf_name) + nnodes, shape3d = gwf.modelgrid.nnodes, gwf.modelgrid.shape + + # get the cell by cell object + cbc = gwf.output.budget() + + # evaluate the full3D option + cbc_eval(cbc, nnodes, shape3d, gwf.modelgrid) + + # clean the run + clean_run(name) + + return + + +def mf2005_eval(name, ws_in): + ml = load_mf2005(name, ws_in) + + # write the model + ml.write_input() + + # run the model + ml.run_model() + + # determine the size of the model grid + nnodes, shape3d = ml.modelgrid.nnodes, ml.modelgrid.shape + + # get the cell by cell object + fpth = os.path.join(tpth, name, f"{name}.cbc") + cbc = flopy.utils.CellBudgetFile(fpth) + + # evaluate the full3D option + cbc_eval(cbc, nnodes, shape3d, ml.modelgrid) + + # clean the run + clean_run(name) + + +def test_cbc_full3D(): + for (name, ismf6, ws_in) in zip(names, ismf6_lst, ex_pths): + if ismf6: + yield mf6_eval, name, ws_in + else: + yield mf2005_eval, name, ws_in + + +def main(): + for (name, ismf6, ws_in) in zip(names, ismf6_lst, ex_pths): + if ismf6: + mf6_eval(name, ws_in) + else: + mf2005_eval(name, ws_in) + + +if __name__ == "__main__": + main() diff --git a/autotest/t505_test.py b/autotest/t505_test.py index 792d206be7..ec914b26a9 100644 --- a/autotest/t505_test.py +++ b/autotest/t505_test.py @@ -139,11 +139,11 @@ def test_np001(): filename="old_name.ims", print_option="ALL", complexity="SIMPLE", - outer_hclose=0.00001, + outer_dvclose=0.00001, outer_maximum=10, under_relaxation="NONE", inner_maximum=10, - inner_hclose=0.001, + inner_dvclose=0.001, linear_acceleration="CG", preconditioner_levels=2, preconditioner_drop_tolerance=0.00001, @@ -156,11 +156,11 @@ def test_np001(): filename=f"{test_ex_name}.ims", print_option="ALL", complexity="SIMPLE", - outer_hclose=0.00001, + outer_dvclose=0.00001, outer_maximum=50, under_relaxation="NONE", inner_maximum=30, - inner_hclose=0.00001, + inner_dvclose=0.00001, linear_acceleration="CG", preconditioner_levels=7, preconditioner_drop_tolerance=0.01, @@ -415,7 +415,8 @@ def test_np001(): assert pymake.compare_heads( None, None, files1=head_file, files2=head_new, outfile=outfile ) - budget_frf = sim.simulation_data.mfdata[(model_name, "CBC", "RIV")] + # budget_frf = sim.simulation_data.mfdata[(model_name, "CBC", "RIV")] + budget_frf = model.output.budget().get_data(text="RIV", full3D=False) assert array_util.riv_array_comp(budget_frf_valid, budget_frf) # clean up @@ -451,7 +452,8 @@ def test_np001(): None, None, files1=head_file, files2=head_new, outfile=outfile ) - budget_frf = sim.simulation_data.mfdata[(model_name, "CBC", "RIV")] + # budget_frf = sim.simulation_data.mfdata[(model_name, "CBC", "RIV")] + budget_frf = model.output.budget().get_data(text="RIV", full3D=False) assert array_util.riv_array_comp(budget_frf_valid, budget_frf) # clean up @@ -660,11 +662,11 @@ def test_np002(): sim, print_option="ALL", complexity="SIMPLE", - outer_hclose=0.00001, + outer_dvclose=0.00001, outer_maximum=50, under_relaxation="NONE", inner_maximum=30, - inner_hclose=0.00001, + inner_dvclose=0.00001, linear_acceleration="CG", preconditioner_levels=7, preconditioner_drop_tolerance=0.01, @@ -909,11 +911,11 @@ def test021_twri(): ims_package = ModflowIms( sim, print_option="SUMMARY", - outer_hclose=0.0001, + outer_dvclose=0.0001, outer_maximum=500, under_relaxation="NONE", inner_maximum=100, - inner_hclose=0.0001, + inner_dvclose=0.0001, rcloserecord=0.001, linear_acceleration="CG", scaling_method="NONE", @@ -1137,11 +1139,11 @@ def test005_advgw_tidal(): sim, print_option="SUMMARY", complexity="SIMPLE", - outer_hclose=0.0001, + outer_dvclose=0.0001, outer_maximum=500, under_relaxation="NONE", inner_maximum=100, - inner_hclose=0.0001, + inner_dvclose=0.0001, rcloserecord=0.001, linear_acceleration="CG", scaling_method="NONE", @@ -1741,11 +1743,11 @@ def test004_bcfss(): print_option="ALL", csv_output_filerecord="bcf2ss.ims.csv", complexity="SIMPLE", - outer_hclose=0.000001, + outer_dvclose=0.000001, outer_maximum=500, under_relaxation="NONE", inner_maximum=100, - inner_hclose=0.000001, + inner_dvclose=0.000001, rcloserecord=0.001, linear_acceleration="CG", scaling_method="NONE", @@ -1895,11 +1897,11 @@ def test035_fhb(): sim, print_option="SUMMARY", complexity="SIMPLE", - outer_hclose=0.001, + outer_dvclose=0.001, outer_maximum=120, under_relaxation="NONE", inner_maximum=100, - inner_hclose=0.0001, + inner_dvclose=0.0001, rcloserecord=0.1, linear_acceleration="CG", preconditioner_levels=7, @@ -2037,11 +2039,11 @@ def test006_gwf3_disv(): ims_package = ModflowIms( sim, print_option="SUMMARY", - outer_hclose=0.00000001, + outer_dvclose=0.00000001, outer_maximum=1000, under_relaxation="NONE", inner_maximum=1000, - inner_hclose=0.00000001, + inner_dvclose=0.00000001, rcloserecord=0.01, linear_acceleration="BICGSTAB", scaling_method="NONE", @@ -2333,11 +2335,11 @@ def test006_2models_gnc(): ims_package = ModflowIms( sim, print_option="SUMMARY", - outer_hclose=0.00000001, + outer_dvclose=0.00000001, outer_maximum=1000, under_relaxation="NONE", inner_maximum=1000, - inner_hclose=0.00000001, + inner_dvclose=0.00000001, rcloserecord=0.01, linear_acceleration="BICGSTAB", scaling_method="NONE", @@ -2669,11 +2671,11 @@ def test050_circle_island(): ims_package = ModflowIms( sim, print_option="SUMMARY", - outer_hclose=0.000001, + outer_dvclose=0.000001, outer_maximum=500, under_relaxation="NONE", inner_maximum=1000, - inner_hclose=0.000001, + inner_dvclose=0.000001, rcloserecord=0.000001, linear_acceleration="BICGSTAB", relaxation_factor=0.0, @@ -2777,7 +2779,7 @@ def test028_sfr(): ims_package = ModflowIms( sim, print_option="SUMMARY", - outer_hclose=0.00001, + outer_dvclose=0.00001, outer_maximum=100, under_relaxation="DBD", under_relaxation_theta=0.85, @@ -2788,7 +2790,7 @@ def test028_sfr(): backtracking_tolerance=1.1, backtracking_reduction_factor=0.7, backtracking_residual_limit=1.0, - inner_hclose=0.00001, + inner_dvclose=0.00001, rcloserecord=0.1, inner_maximum=100, linear_acceleration="CG", diff --git a/examples/data/mf6/test001a_Tharmonic/flow15.ims b/examples/data/mf6/test001a_Tharmonic/flow15.ims index 72a5dc1b9a..7a03edd57b 100644 --- a/examples/data/mf6/test001a_Tharmonic/flow15.ims +++ b/examples/data/mf6/test001a_Tharmonic/flow15.ims @@ -6,21 +6,21 @@ BEGIN OPTIONS END OPTIONS BEGIN NONLINEAR - OUTER_HCLOSE 0.99999997E-05 + OUTER_DVCLOSE 0.99999997E-05 OUTER_MAXIMUM 50 UNDER_RELAXATION NONE END NONLINEAR BEGIN LINEAR INNER_MAXIMUM 30 - INNER_HCLOSE 0.10000000E-05 + INNER_DVCLOSE 0.10000000E-05 INNER_RCLOSE 0.99999997E-05 LINEAR_ACCELERATION CG END LINEAR BEGIN XMD INNER_MAXIMUM 30 - INNER_HCLOSE 0.10000000E-05 + INNER_DVCLOSE 0.10000000E-05 LINEAR_ACCELERATION CG PRECONDITIONER_LEVELS 7 PRECONDITIONER_DROP_TOLERANCE 0.10000000E-02 diff --git a/examples/data/mf6/test001e_UZF_3lay/test001e_UZF_3lay.ims b/examples/data/mf6/test001e_UZF_3lay/test001e_UZF_3lay.ims index 5004d3f9f6..0d0e8b638c 100644 --- a/examples/data/mf6/test001e_UZF_3lay/test001e_UZF_3lay.ims +++ b/examples/data/mf6/test001e_UZF_3lay/test001e_UZF_3lay.ims @@ -6,14 +6,14 @@ BEGIN Options END Options BEGIN Nonlinear - OUTER_HCLOSE 1.E-06 + OUTER_DVCLOSE 1.E-06 OUTER_MAXIMUM 100 UNDER_RELAXATION dbd END Nonlinear BEGIN LINEAR INNER_MAXIMUM 50 - INNER_HCLOSE 1.E-09 + INNER_DVCLOSE 1.E-09 INNER_RCLOSE 1.E-05 LINEAR_ACCELERATION BICGSTAB END LINEAR diff --git a/examples/data/mf6/test003_gwfs_disv/model.ims b/examples/data/mf6/test003_gwfs_disv/model.ims index 19d342427b..23020b0b3f 100644 --- a/examples/data/mf6/test003_gwfs_disv/model.ims +++ b/examples/data/mf6/test003_gwfs_disv/model.ims @@ -3,13 +3,13 @@ begin options end options begin nonlinear - outer_hclose 1.e-4 - outer_maximum 500 + OUTER_DVCLOSE 1.e-4 + outer_maximum 500 under_relaxation none end nonlinear begin linear - inner_hclose 1.0e-4 + INNER_DVCLOSE 1.0e-4 inner_rclose 0.001 #L2NORM_RCLOSE inner_maximum 100 @@ -28,13 +28,13 @@ end linear 1.0E-4 1.0E-4 500 100 1 0 001 #hclose, hiclose,mxiter,iter1,iprsms,nonmeth,linmeth 2 0 0 2 0 0 0 1e-3 -IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN +IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN 0.0001 0.0001 500 100 1 0 1 2 0 0 2 0 0 0 1e-3 -IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN +IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN 0.0001 0.0001 500 100 1 0 5 3 0 0 100000 diff --git a/examples/data/mf6/test003_gwftri_disv/model.ims b/examples/data/mf6/test003_gwftri_disv/model.ims index 19d342427b..23020b0b3f 100644 --- a/examples/data/mf6/test003_gwftri_disv/model.ims +++ b/examples/data/mf6/test003_gwftri_disv/model.ims @@ -3,13 +3,13 @@ begin options end options begin nonlinear - outer_hclose 1.e-4 - outer_maximum 500 + OUTER_DVCLOSE 1.e-4 + outer_maximum 500 under_relaxation none end nonlinear begin linear - inner_hclose 1.0e-4 + INNER_DVCLOSE 1.0e-4 inner_rclose 0.001 #L2NORM_RCLOSE inner_maximum 100 @@ -28,13 +28,13 @@ end linear 1.0E-4 1.0E-4 500 100 1 0 001 #hclose, hiclose,mxiter,iter1,iprsms,nonmeth,linmeth 2 0 0 2 0 0 0 1e-3 -IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN +IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN 0.0001 0.0001 500 100 1 0 1 2 0 0 2 0 0 0 1e-3 -IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN +IACL NORDER LEVEL NORTH IREDSYS RRCTOL IDROPTOL EPSRN 0.0001 0.0001 500 100 1 0 5 3 0 0 100000 diff --git a/examples/data/mf6/test005_advgw_tidal/model.ims b/examples/data/mf6/test005_advgw_tidal/model.ims index a4a285c6ad..36e547e2a8 100644 --- a/examples/data/mf6/test005_advgw_tidal/model.ims +++ b/examples/data/mf6/test005_advgw_tidal/model.ims @@ -4,14 +4,14 @@ BEGIN options END options BEGIN nonlinear - OUTER_HCLOSE 0.0001 + OUTER_DVCLOSE 0.0001 OUTER_MAXIMUM 500 UNDER_RELAXATION none END nonlinear BEGIN linear INNER_MAXIMUM 100 - INNER_HCLOSE 0.0001 + INNER_DVCLOSE 0.0001 inner_rclose 0.001 LINEAR_ACCELERATION cg RELAXATION_FACTOR 0.97 diff --git a/examples/data/mf6/test006_2models_mvr/simulation.ims b/examples/data/mf6/test006_2models_mvr/simulation.ims index ed52ca4f2a..6a15577531 100644 --- a/examples/data/mf6/test006_2models_mvr/simulation.ims +++ b/examples/data/mf6/test006_2models_mvr/simulation.ims @@ -3,13 +3,13 @@ begin options end options begin nonlinear - outer_hclose 1.e-8 - outer_maximum 1000 + OUTER_DVCLOSE 1.e-8 + outer_maximum 1000 under_relaxation none end nonlinear begin linear - inner_hclose 1.0e-8 + INNER_DVCLOSE 1.0e-8 inner_rclose 0.01 inner_maximum 1000 linear_acceleration bicgstab diff --git a/examples/data/mf6/test006_gwf3/flow.ims b/examples/data/mf6/test006_gwf3/flow.ims index d8045187ee..f01fd7a923 100644 --- a/examples/data/mf6/test006_gwf3/flow.ims +++ b/examples/data/mf6/test006_gwf3/flow.ims @@ -3,13 +3,13 @@ begin options end options begin nonlinear - outer_hclose 1.e-8 - outer_maximum 1000 + OUTER_DVCLOSE 1.e-8 + outer_maximum 1000 under_relaxation none end nonlinear begin linear - inner_hclose 1.0e-8 + INNER_DVCLOSE 1.0e-8 inner_rclose 0.01 inner_maximum 1000 linear_acceleration cg diff --git a/examples/data/mf6/test027_TimeseriesTest/model.ims b/examples/data/mf6/test027_TimeseriesTest/model.ims index b1ab693de6..55c527ebbf 100644 --- a/examples/data/mf6/test027_TimeseriesTest/model.ims +++ b/examples/data/mf6/test027_TimeseriesTest/model.ims @@ -3,8 +3,8 @@ begin options end options begin nonlinear - outer_hclose 1.e-9 - outer_maximum 50 + OUTER_DVCLOSE 1.e-9 + outer_maximum 50 under_relaxation none under_relaxation_theta 0.9 under_relaxation_kappa 0.100000E-03 @@ -18,7 +18,7 @@ end nonlinear begin linear - inner_hclose 1.0e-9 + INNER_DVCLOSE 1.0e-9 inner_rclose 0.001 inner_maximum 100 linear_acceleration cg diff --git a/examples/data/mf6/test036_twrihfb/twrihfb2015.ims b/examples/data/mf6/test036_twrihfb/twrihfb2015.ims index eca52e9732..2ecae27927 100644 --- a/examples/data/mf6/test036_twrihfb/twrihfb2015.ims +++ b/examples/data/mf6/test036_twrihfb/twrihfb2015.ims @@ -6,14 +6,14 @@ BEGIN Options END Options BEGIN Nonlinear - OUTER_HCLOSE 0.10000000E-02 + OUTER_DVCLOSE 0.10000000E-02 OUTER_MAXIMUM 50 UNDER_RELAXATION NONE END Nonlinear BEGIN LINEAR INNER_MAXIMUM 100 - INNER_HCLOSE 0.10000000E-03 - INNER_RCLOSE 0.10000000 + INNER_DVCLOSE 0.10000000E-03 + INNER_RCLOSE 0.10000000 LINEAR_ACCELERATION CG END LINEAR diff --git a/examples/data/mf6/test045_lake1ss_table/lakeex1b.ims b/examples/data/mf6/test045_lake1ss_table/lakeex1b.ims index a842a21936..cc8719a394 100644 --- a/examples/data/mf6/test045_lake1ss_table/lakeex1b.ims +++ b/examples/data/mf6/test045_lake1ss_table/lakeex1b.ims @@ -3,7 +3,7 @@ BEGIN Options END Options BEGIN Nonlinear - OUTER_HCLOSE 0.0001 + OUTER_DVCLOSE 0.0001 OUTER_MAXIMUM 500 UNDER_RELAXATION none UNDER_RELAXATION_THETA 0.0 @@ -14,7 +14,7 @@ END Nonlinear BEGIN LINEAR INNER_MAXIMUM 100 - INNER_HCLOSE 0.0001 + INNER_DVCLOSE 0.0001 INNER_RCLOSE 0.1 LINEAR_ACCELERATION cg RELAXATION_FACTOR 1.0 diff --git a/examples/data/mf6/test045_lake2tr/lakeex2a.ims b/examples/data/mf6/test045_lake2tr/lakeex2a.ims index c67b4b46a3..708a04115c 100644 --- a/examples/data/mf6/test045_lake2tr/lakeex2a.ims +++ b/examples/data/mf6/test045_lake2tr/lakeex2a.ims @@ -6,18 +6,18 @@ BEGIN Options END Options BEGIN Nonlinear - OUTER_HCLOSE 1e-4 + OUTER_DVCLOSE 1e-4 OUTER_MAXIMUM 500 UNDER_RELAXATION NONE - UNDER_RELAXATION_THETA 0.000000 - UNDER_RELAXATION_KAPPA 0.000000 - UNDER_RELAXATION_GAMMA 0.000000 - UNDER_RELAXATION_MOMENTUM 0.000000 + UNDER_RELAXATION_THETA 0.000000 + UNDER_RELAXATION_KAPPA 0.000000 + UNDER_RELAXATION_GAMMA 0.000000 + UNDER_RELAXATION_MOMENTUM 0.000000 END Nonlinear BEGIN LINEAR INNER_MAXIMUM 100 - INNER_HCLOSE 1e-4 + INNER_DVCLOSE 1e-4 INNER_RCLOSE 1e-3 LINEAR_ACCELERATION CG #PRECONDITIONER_LEVELS 7 diff --git a/flopy/utils/binaryfile.py b/flopy/utils/binaryfile.py index f04b2ac296..d1503a4122 100755 --- a/flopy/utils/binaryfile.py +++ b/flopy/utils/binaryfile.py @@ -674,6 +674,14 @@ def __init__(self, filename, precision="auto", verbose=False, **kwargs): else: raise Exception(f"Unknown precision specified: {precision}") + # set shape for full3D option + if self.modelgrid is None: + self.shape = (self.nlay, self.nrow, self.ncol) + self.nnodes = self.nlay * self.nrow * self.ncol + else: + self.shape = self.modelgrid.shape + self.nnodes = self.modelgrid.nnodes + if not success: raise Exception( f"Budget file could not be read using {precision} precision" @@ -788,19 +796,24 @@ def _build_index(self): for i in range(33, 127): asciiset += chr(i) + # read first record header = self._get_header() - self.nrow = header["nrow"] - self.ncol = header["ncol"] - self.nlay = np.abs(header["nlay"]) + nrow = header["nrow"] + ncol = header["ncol"] text = header["text"] if isinstance(text, bytes): text = text.decode() - if self.nrow < 0 or self.ncol < 0: + if nrow < 0 or ncol < 0: raise Exception("negative nrow, ncol") + if not text.endswith("FLOW-JA-FACE"): + self.nrow = nrow + self.ncol = ncol + self.nlay = np.abs(header["nlay"]) self.file.seek(0, 2) self.totalbytes = self.file.tell() self.file.seek(0, 0) self.recorddict = {} + # read the remaining records ipos = 0 while ipos < self.totalbytes: self.iposheader.append(ipos) @@ -865,6 +878,16 @@ def _build_index(self): ): print("") + # set the nrow, ncol, and nlay if they have not been set + if self.nrow == 0: + text = header["text"] + if isinstance(text, bytes): + text = text.decode() + if not text.endswith("FLOW-JA-FACE"): + self.nrow = header["nrow"] + self.ncol = header["ncol"] + self.nlay = np.abs(header["nlay"]) + # store record and byte position mapping self.recorddict[ tuple(header) @@ -1548,13 +1571,16 @@ def get_record(self, idx, full3D=False): dtype = np.dtype([("node", np.int32), ("q", self.realtype)]) if self.verbose: if full3D: - s += f"a numpy masked array of size ({nlay}, {nrow}, {ncol})" + s += ( + f"a numpy masked array of " + f"size ({nlay}, {nrow}, {ncol})" + ) else: s += f"a numpy recarray of size ({nlist}, 2)" print(s) data = binaryread(self.file, dtype, shape=(nlist,)) if full3D: - return self.create3D(data, nlay, nrow, ncol) + return self.__create3D(data) else: return data.view(np.recarray) @@ -1564,23 +1590,27 @@ def get_record(self, idx, full3D=False): data = binaryread(self.file, self.realtype(1), shape=(nrow, ncol)) if self.verbose: if full3D: - s += f"a numpy masked array of size ({nlay}, {nrow}, {ncol})" + s += ( + "a numpy masked array of size " + f"({nlay}, {nrow}, {ncol})" + ) else: s += ( - "a list of two 2D numpy arrays. " - "The first is an integer layer array of shape {}. " - "The second is real data array of shape {}".format( - (nrow, ncol), - (nrow, ncol), - ) + "a list of two 2D numpy arrays. The first is an " + f"integer layer array of shape ({nrow}, {ncol}). The " + f"second is real data array of shape ({nrow}, {ncol})" ) print(s) if full3D: - out = np.ma.zeros((nlay, nrow, ncol), dtype=np.float32) + out = np.ma.zeros(self.nnodes, dtype=np.float32) out.mask = True - vertical_layer = ilayer[0] - 1 # This is always the top layer - out[vertical_layer, :, :] = data - return out + vertical_layer = ilayer.flatten() - 1 + # create the 2D cell index and then move it to + # the correct vertical location + idx = np.arange(0, vertical_layer.shape[0]) + idx += vertical_layer * nrow * ncol + out[idx] = data.flatten() + return out.reshape(self.shape) else: return [ilayer, data] @@ -1608,7 +1638,7 @@ def get_record(self, idx, full3D=False): if self.verbose: s += f"a list array of shape ({nlay}, {nrow}, {ncol})" print(s) - return self.create3D(data, nlay, nrow, ncol) + return self.__create3D(data) else: if self.verbose: s += f"a numpy recarray of size ({nlist}, {2 + naux})" @@ -1631,13 +1661,12 @@ def get_record(self, idx, full3D=False): data = binaryread(self.file, dtype, shape=(nlist,)) if self.verbose: if full3D: - s += f"full 3D arrays not supported for imeth = {imeth}" + s += f"a list array of shape ({nlay}, {nrow}, {ncol})" else: s += f"a numpy recarray of size ({nlist}, 2)" print(s) if full3D: - s += f"full 3D arrays not supported for imeth = {imeth}" - raise ValueError(s) + return self.__create3D(data) else: return data.view(np.recarray) else: @@ -1646,34 +1675,30 @@ def get_record(self, idx, full3D=False): # should not reach this point return - def create3D(self, data, nlay, nrow, ncol): + def __create3D(self, data): """ Convert a dictionary of {node: q, ...} into a numpy masked array. - In most cases this should not be called directly by the user unless - you know what you're doing. Instead, it is used as part of the - full3D keyword for get_data. + Used to create full grid arrays when the full3D keyword is set + to True in get_data. Parameters ---------- data : dictionary Dictionary with node keywords and flows (q) items. - nlay, nrow, ncol : int - Number of layers, rows, and columns of the model grid. - Returns ---------- out : numpy masked array List contains unique simulation times (totim) in binary file. """ - out = np.ma.zeros((nlay * nrow * ncol), dtype=np.float32) + out = np.ma.zeros(self.nnodes, dtype=np.float32) out.mask = True for [node, q] in zip(data["node"], data["q"]): idx = node - 1 out.data[idx] += q out.mask[idx] = False - return np.ma.reshape(out, (nlay, nrow, ncol)) + return np.ma.reshape(out, self.shape) def get_times(self): """