From 1d8b8458c584e9ef60d1389b249f13d581838515 Mon Sep 17 00:00:00 2001 From: Joseph D Hughes Date: Thu, 30 Sep 2021 18:07:01 -0500 Subject: [PATCH 1/5] feat(CellBudget): add support for full3D keyword --- autotest/t004_test_utilarray.py | 18 +- autotest/t050_test.py | 35 +-- ..._ugridtests.py => t075_test_ugridtests.py} | 0 ..._thick.py => t076_test_modelgrid_thick.py} | 0 ...tions.py => t078_test_lake_connections.py} | 0 autotest/t079_test_cbc_full3D.py | 204 ++++++++++++++++++ .../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, 328 insertions(+), 88 deletions(-) rename autotest/{t075_ugridtests.py => t075_test_ugridtests.py} (100%) 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/t075_ugridtests.py b/autotest/t075_test_ugridtests.py similarity index 100% rename from autotest/t075_ugridtests.py rename to autotest/t075_test_ugridtests.py 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..4fb9352ca0 --- /dev/null +++ b/autotest/t079_test_cbc_full3D.py @@ -0,0 +1,204 @@ +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) + + cobj_mg = flopy.utils.CellBudgetFile( + cbc_pth, + modelgrid=modelgrid, + verbose=True, + ) + cbc_eval_size(cobj_mg, nnodes, shape3d) + cbc_eval_data(cobj_mg, shape3d) + + 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: + continue + mf6_eval(name, ws_in) + else: + mf2005_eval(name, ws_in) + + +if __name__ == "__main__": + main() 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): """ From e21e5e153ad500826afea23b19627c70b5fdd1b9 Mon Sep 17 00:00:00 2001 From: Joseph D Hughes Date: Fri, 1 Oct 2021 11:57:59 -0500 Subject: [PATCH 2/5] feat(CellBudget): add support for full3D keyword --- autotest/t079_test_cbc_full3D.py | 3 +- autotest/t505_test.py | 50 +++++++++++++++++--------------- 2 files changed, 28 insertions(+), 25 deletions(-) diff --git a/autotest/t079_test_cbc_full3D.py b/autotest/t079_test_cbc_full3D.py index 4fb9352ca0..e0025a483d 100644 --- a/autotest/t079_test_cbc_full3D.py +++ b/autotest/t079_test_cbc_full3D.py @@ -115,6 +115,7 @@ 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, @@ -123,6 +124,7 @@ def cbc_eval(cbcobj, nnodes, shape3d, modelgrid): ) cbc_eval_size(cobj_mg, nnodes, shape3d) cbc_eval_data(cobj_mg, shape3d) + cobj_mg.close() return @@ -194,7 +196,6 @@ def test_cbc_full3D(): def main(): for (name, ismf6, ws_in) in zip(names, ismf6_lst, ex_pths): if ismf6: - continue mf6_eval(name, ws_in) else: mf2005_eval(name, ws_in) 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", From 40814ca33521b65bf087b24912332cf7cfcc3b41 Mon Sep 17 00:00:00 2001 From: Joseph D Hughes Date: Fri, 1 Oct 2021 14:07:40 -0500 Subject: [PATCH 3/5] feat(CellBudget): add support for full3D keyword --- autotest/t075_test_ugridtests.py | 6 +++++- flopy/utils/voronoi.py | 1 + 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/autotest/t075_test_ugridtests.py b/autotest/t075_test_ugridtests.py index 14fecb01eb..3836768ebf 100644 --- a/autotest/t075_test_ugridtests.py +++ b/autotest/t075_test_ugridtests.py @@ -225,7 +225,11 @@ def test_triangle_unstructured_grid(): xc, yc = tri.get_xcyc().T ncpl = np.array([len(iverts)]) g = UnstructuredGrid( - vertices=verts, iverts=iverts, ncpl=ncpl, xcenters=xc, ycenters=yc, + vertices=verts, + iverts=iverts, + ncpl=ncpl, + xcenters=xc, + ycenters=yc, ) assert len(g.grid_lines) == 8190 assert g.nnodes == g.ncpl == 2730 diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index 4740f43c8e..14c6e445a3 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -349,3 +349,4 @@ def plot(self, ax=None, plot_title=True, **kwargs): if plot_title: ax.set_title(f"ncells: {self.ncpl}; nverts: {self.nverts}") return ax + \ No newline at end of file From 0f4b3f10f3c7384151efb7c8d4e4a740dfb6eac1 Mon Sep 17 00:00:00 2001 From: Joseph D Hughes Date: Fri, 1 Oct 2021 14:11:56 -0500 Subject: [PATCH 4/5] feat(CellBudget): add support for full3D keyword --- flopy/utils/voronoi.py | 1 - 1 file changed, 1 deletion(-) diff --git a/flopy/utils/voronoi.py b/flopy/utils/voronoi.py index 14c6e445a3..4740f43c8e 100644 --- a/flopy/utils/voronoi.py +++ b/flopy/utils/voronoi.py @@ -349,4 +349,3 @@ def plot(self, ax=None, plot_title=True, **kwargs): if plot_title: ax.set_title(f"ncells: {self.ncpl}; nverts: {self.nverts}") return ax - \ No newline at end of file From 493cdab0de0e691de8360e420eee8d711b38519e Mon Sep 17 00:00:00 2001 From: Joseph D Hughes Date: Mon, 4 Oct 2021 21:29:08 -0500 Subject: [PATCH 5/5] feat(CellBudget): add support for full3D keyword --- autotest/t075_test_ugridtests.py | 541 ------------------------------- 1 file changed, 541 deletions(-) delete mode 100644 autotest/t075_test_ugridtests.py diff --git a/autotest/t075_test_ugridtests.py b/autotest/t075_test_ugridtests.py deleted file mode 100644 index 8cf1ba3847..0000000000 --- a/autotest/t075_test_ugridtests.py +++ /dev/null @@ -1,541 +0,0 @@ -""" -Unstructured grid tests - -""" - -import os -import numpy as np -from flopy.discretization import UnstructuredGrid, VertexGrid -from flopy.utils.triangle import Triangle -from flopy.utils.voronoi import VoronoiGrid - -tpth = os.path.join("temp", "t075") -if not os.path.isdir(tpth): - os.mkdir(tpth) - - -def test_unstructured_grid_shell(): - # constructor with no arguments. incomplete shell should exist - g = UnstructuredGrid() - assert g.nlay is None - assert g.nnodes is None - assert g.ncpl is None - assert not g.grid_varies_by_layer - assert not g.is_valid - assert not g.is_complete - return - - -def test_unstructured_grid_dimensions(): - # constructor with just dimensions - ncpl = [1, 10, 1] - g = UnstructuredGrid(ncpl=ncpl) - assert np.allclose(g.ncpl, ncpl) - assert g.nlay == 3 - assert g.nnodes == 12 - assert not g.is_valid - assert not g.is_complete - assert not g.grid_varies_by_layer - return - - -def test_unstructured_minimal_grid(): - - # pass in simple 2 cell minimal grid to make grid valid - vertices = [ - [0, 0.0, 1.0], - [1, 1.0, 1.0], - [2, 2.0, 1.0], - [3, 0.0, 0.0], - [4, 1.0, 0.0], - [5, 2.0, 0.0], - ] - iverts = [[0, 1, 4, 3], [1, 2, 5, 4]] - xcenters = [0.5, 1.5] - ycenters = [0.5, 0.5] - g = UnstructuredGrid( - vertices=vertices, iverts=iverts, xcenters=xcenters, ycenters=ycenters - ) - assert np.allclose(g.ncpl, np.array([2], dtype=int)) - assert g.nlay == 1 - assert g.nnodes == 2 - assert g.is_valid - assert not g.is_complete - assert not g.grid_varies_by_layer - assert g._vertices == vertices - assert g._iverts == iverts - assert g._xc == xcenters - assert g._yc == ycenters - grid_lines = [ - [(0.0, 0), (0.0, 1.0)], - [(0.0, 1), (1.0, 1.0)], - [(1.0, 1), (1.0, 0.0)], - [(1.0, 0), (0.0, 0.0)], - [(1.0, 0), (1.0, 1.0)], - [(1.0, 1), (2.0, 1.0)], - [(2.0, 1), (2.0, 0.0)], - [(2.0, 0), (1.0, 0.0)], - ] - assert ( - g.grid_lines == grid_lines - ), f"\n{g.grid_lines} \n /= \n{grid_lines}" - assert g.extent == (0, 2, 0, 1) - xv, yv, zv = g.xyzvertices - assert xv == [[0, 1, 1, 0], [1, 2, 2, 1]] - assert yv == [[1, 1, 0, 0], [1, 1, 0, 0]] - assert zv is None - - return - - -def test_unstructured_complete_grid(): - - # pass in simple 2 cell complete grid to make grid valid, and put each - # cell in a different layer - vertices = [ - [0, 0.0, 1.0], - [1, 1.0, 1.0], - [2, 2.0, 1.0], - [3, 0.0, 0.0], - [4, 1.0, 0.0], - [5, 2.0, 0.0], - ] - iverts = [[0, 1, 4, 3], [1, 2, 5, 4]] - xcenters = [0.5, 1.5] - ycenters = [0.5, 0.5] - ncpl = [1, 1] - top = [1, 0] - top = np.array(top) - botm = [0, -1] - botm = np.array(botm) - g = UnstructuredGrid( - vertices=vertices, - iverts=iverts, - xcenters=xcenters, - ycenters=ycenters, - ncpl=ncpl, - top=top, - botm=botm, - ) - assert np.allclose(g.ncpl, np.array([1, 1], dtype=int)) - assert g.nlay == 2 - assert g.nnodes == 2 - assert g.is_valid - assert not g.is_complete - assert g.grid_varies_by_layer - assert g._vertices == vertices - assert g._iverts == iverts - assert g._xc == xcenters - assert g._yc == ycenters - grid_lines = { - 0: [ - [(0.0, 0.0), (0.0, 1.0)], - [(0.0, 1.0), (1.0, 1.0)], - [(1.0, 1.0), (1.0, 0.0)], - [(1.0, 0.0), (0.0, 0.0)], - ], - 1: [ - [(1.0, 0.0), (1.0, 1.0)], - [(1.0, 1.0), (2.0, 1.0)], - [(2.0, 1.0), (2.0, 0.0)], - [(2.0, 0.0), (1.0, 0.0)], - ], - } - assert isinstance(g.grid_lines, dict) - assert ( - g.grid_lines == grid_lines - ), f"\n{g.grid_lines} \n /= \n{grid_lines}" - assert g.extent == (0, 2, 0, 1) - xv, yv, zv = g.xyzvertices - assert xv == [[0, 1, 1, 0], [1, 2, 2, 1]] - assert yv == [[1, 1, 0, 0], [1, 1, 0, 0]] - assert np.allclose(zv, np.array([[1, 0], [0, -1]])) - - return - - -def test_loading_argus_meshes(): - datapth = os.path.join("..", "examples", "data", "unstructured") - fnames = [fname for fname in os.listdir(datapth) if fname.endswith(".exp")] - for fname in fnames: - fname = os.path.join(datapth, fname) - print(f"Loading Argus mesh ({fname}) into UnstructuredGrid") - g = UnstructuredGrid.from_argus_export(fname) - print(f" Number of nodes: {g.nnodes}") - - -def test_create_unstructured_grid_from_verts(): - - datapth = os.path.join("..", "examples", "data", "unstructured") - - # simple functions to load vertices and incidence lists - def load_verts(fname): - print(f"Loading vertices from: {fname}") - verts = np.genfromtxt( - fname, dtype=[int, float, float], names=["iv", "x", "y"] - ) - verts["iv"] -= 1 # zero based - return verts - - def load_iverts(fname): - print(f"Loading iverts from: {fname}") - f = open(fname, "r") - iverts = [] - xc = [] - yc = [] - for line in f: - ll = line.strip().split() - iverts.append([int(i) - 1 for i in ll[4:]]) - xc.append(float(ll[1])) - yc.append(float(ll[2])) - return iverts, np.array(xc), np.array(yc) - - # load vertices - fname = os.path.join(datapth, "ugrid_verts.dat") - verts = load_verts(fname) - - # load the incidence list into iverts - fname = os.path.join(datapth, "ugrid_iverts.dat") - iverts, xc, yc = load_iverts(fname) - - ncpl = np.array(5 * [len(iverts)]) - g = UnstructuredGrid(verts, iverts, xc, yc, ncpl=ncpl) - assert isinstance(g.grid_lines, list) - assert np.allclose(g.ncpl, ncpl) - assert g.extent == (0.0, 700.0, 0.0, 700.0) - assert g._vertices.shape == (156,) - assert g.nnodes == g.ncpl.sum() == 1090 - return - - -def test_triangle_unstructured_grid(): - maximum_area = 30000.0 - extent = (214270.0, 221720.0, 4366610.0, 4373510.0) - domainpoly = [ - (extent[0], extent[2]), - (extent[1], extent[2]), - (extent[1], extent[3]), - (extent[0], extent[3]), - ] - tri = Triangle(maximum_area=maximum_area, angle=30, model_ws=tpth) - tri.add_polygon(domainpoly) - tri.build(verbose=False) - verts = [[iv, x, y] for iv, (x, y) in enumerate(tri.verts)] - iverts = tri.iverts - xc, yc = tri.get_xcyc().T - ncpl = np.array([len(iverts)]) - g = UnstructuredGrid( - vertices=verts, - iverts=iverts, - ncpl=ncpl, - xcenters=xc, - ycenters=yc, - ) - assert len(g.grid_lines) == 8190 - assert g.nnodes == g.ncpl == 2730 - return - - -def test_voronoi_vertex_grid(): - xmin = 0.0 - xmax = 2.0 - ymin = 0.0 - ymax = 1.0 - area_max = 0.05 - tri = Triangle(maximum_area=area_max, angle=30, model_ws=tpth) - poly = np.array(((xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax))) - tri.add_polygon(poly) - tri.build(verbose=False) - - # create vor object and VertexGrid - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - vgrid = VertexGrid(**gridprops, nlay=1) - assert vgrid.is_valid - - # arguments for creating a mf6 disv package - gridprops = vor.get_disv_gridprops() - print(gridprops) - assert gridprops["ncpl"] == 43 - assert gridprops["nvert"] == 127 - assert len(gridprops["vertices"]) == 127 - assert len(gridprops["cell2d"]) == 43 - return - - -def test_voronoi_grid0(plot=False): - - name = "vor0" - answer_ncpl = 3804 - domain = [ - [1831.381546, 6335.543757], - [4337.733475, 6851.136153], - [6428.747084, 6707.916043], - [8662.980804, 6493.085878], - [9350.437333, 5891.561415], - [9235.861245, 4717.156511], - [8963.743036, 3685.971717], - [8691.624826, 2783.685023], - [8047.13433, 2038.94045], - [7416.965845, 578.0953252], - [6414.425073, 105.4689614], - [5354.596258, 205.7230386], - [4624.173696, 363.2651598], - [3363.836725, 563.7733141], - [1330.11116, 1809.788273], - [399.1804436, 2998.515188], - [914.7728404, 5132.494831], - [1831.381546, 6335.543757], - ] - area_max = 100.0 ** 2 - tri = Triangle(maximum_area=area_max, angle=30, model_ws=tpth) - poly = np.array(domain) - tri.add_polygon(poly) - tri.build(verbose=False) - - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - ncpl = gridprops["ncpl"] - assert ( - ncpl == answer_ncpl - ), f"Number of cells should be {answer_ncpl}. Found {ncpl}" - - voronoi_grid = VertexGrid(**gridprops, nlay=1) - - if plot: - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot() - ax.set_aspect("equal") - voronoi_grid.plot(ax=ax) - plt.savefig(os.path.join(tpth, f"{name}.png")) - - return - - -def test_voronoi_grid1(plot=False): - - name = "vor1" - answer_ncpl = 1679 - xmin = 0.0 - xmax = 2.0 - ymin = 0.0 - ymax = 1.0 - area_max = 0.001 - tri = Triangle(maximum_area=area_max, angle=30, model_ws=tpth) - poly = np.array(((xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax))) - tri.add_polygon(poly) - tri.build(verbose=False) - - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - voronoi_grid = VertexGrid(**gridprops, nlay=1) - ncpl = gridprops["ncpl"] - assert ( - ncpl == answer_ncpl - ), f"Number of cells should be {answer_ncpl}. Found {ncpl}" - - if plot: - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot() - ax.set_aspect("equal") - voronoi_grid.plot(ax=ax) - plt.savefig(os.path.join(tpth, f"{name}.png")) - - return - - -def test_voronoi_grid2(plot=False): - - name = "vor2" - answer_ncpl = 538 - theta = np.arange(0.0, 2 * np.pi, 0.2) - radius = 100.0 - x = radius * np.cos(theta) - y = radius * np.sin(theta) - circle_poly = [(x, y) for x, y in zip(x, y)] - tri = Triangle(maximum_area=50, angle=30, model_ws=tpth) - tri.add_polygon(circle_poly) - tri.build(verbose=False) - - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - voronoi_grid = VertexGrid(**gridprops, nlay=1) - ncpl = gridprops["ncpl"] - assert ( - ncpl == answer_ncpl - ), f"Number of cells should be {answer_ncpl}. Found {ncpl}" - - if plot: - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot() - ax.set_aspect("equal") - voronoi_grid.plot(ax=ax) - plt.savefig(os.path.join(tpth, f"{name}.png")) - - return - - -def test_voronoi_grid3(plot=False): - - name = "vor3" - answer_ncpl = 300 - - theta = np.arange(0.0, 2 * np.pi, 0.2) - radius = 100.0 - x = radius * np.cos(theta) - y = radius * np.sin(theta) - circle_poly = [(x, y) for x, y in zip(x, y)] - - theta = np.arange(0.0, 2 * np.pi, 0.2) - radius = 30.0 - x = radius * np.cos(theta) + 25.0 - y = radius * np.sin(theta) + 25.0 - inner_circle_poly = [(x, y) for x, y in zip(x, y)] - - tri = Triangle(maximum_area=100, angle=30, model_ws=tpth) - tri.add_polygon(circle_poly) - tri.add_polygon(inner_circle_poly) - tri.add_hole((25, 25)) - tri.build(verbose=False) - - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - voronoi_grid = VertexGrid(**gridprops, nlay=1) - ncpl = gridprops["ncpl"] - assert ( - ncpl == answer_ncpl - ), f"Number of cells should be {answer_ncpl}. Found {ncpl}" - - if plot: - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot() - ax.set_aspect("equal") - voronoi_grid.plot(ax=ax) - plt.savefig(os.path.join(tpth, f"{name}.png")) - - return - - -def test_voronoi_grid4(plot=False): - - name = "vor4" - answer_ncpl = 410 - active_domain = [(0, 0), (100, 0), (100, 100), (0, 100)] - area1 = [(10, 10), (40, 10), (40, 40), (10, 40)] - area2 = [(60, 60), (80, 60), (80, 80), (60, 80)] - tri = Triangle(angle=30, model_ws=tpth) - tri.add_polygon(active_domain) - tri.add_polygon(area1) - tri.add_polygon(area2) - tri.add_region((1, 1), 0, maximum_area=100) # point inside active domain - tri.add_region((11, 11), 1, maximum_area=10) # point inside area1 - tri.add_region((61, 61), 2, maximum_area=3) # point inside area2 - tri.build(verbose=False) - - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - voronoi_grid = VertexGrid(**gridprops, nlay=1) - ncpl = gridprops["ncpl"] - assert ( - ncpl == answer_ncpl - ), f"Number of cells should be {answer_ncpl}. Found {ncpl}" - - if plot: - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot() - ax.set_aspect("equal") - voronoi_grid.plot(ax=ax) - plt.savefig(os.path.join(tpth, f"{name}.png")) - - return - - -def test_voronoi_grid5(plot=False): - - name = "vor5" - answer_ncpl = 1067 - active_domain = [(0, 0), (100, 0), (100, 100), (0, 100)] - area1 = [(10, 10), (40, 10), (40, 40), (10, 40)] - area2 = [(50, 50), (90, 50), (90, 90), (50, 90)] - - tri = Triangle(angle=30, model_ws=tpth) - - # requirement that active_domain is first polygon to be added - tri.add_polygon(active_domain) - - # requirement that any holes be added next - theta = np.arange(0.0, 2 * np.pi, 0.2) - radius = 10.0 - x = radius * np.cos(theta) + 50.0 - y = radius * np.sin(theta) + 70.0 - circle_poly0 = [(x, y) for x, y in zip(x, y)] - tri.add_polygon(circle_poly0) - tri.add_hole((50, 70)) - - # Add a polygon to force cells to conform to it - theta = np.arange(0.0, 2 * np.pi, 0.2) - radius = 10.0 - x = radius * np.cos(theta) + 70.0 - y = radius * np.sin(theta) + 20.0 - circle_poly1 = [(x, y) for x, y in zip(x, y)] - tri.add_polygon(circle_poly1) - # tri.add_hole((70, 20)) - - # add line through domain to force conforming cells - line = [(x, x) for x in np.linspace(10, 90, 100)] - tri.add_polygon(line) - - # then regions and other polygons should follow - tri.add_polygon(area1) - tri.add_polygon(area2) - tri.add_region((1, 1), 0, maximum_area=100) # point inside active domain - tri.add_region((11, 11), 1, maximum_area=10) # point inside area1 - tri.add_region((70, 70), 2, maximum_area=3) # point inside area2 - - tri.build(verbose=False) - - vor = VoronoiGrid(tri) - gridprops = vor.get_gridprops_vertexgrid() - voronoi_grid = VertexGrid(**gridprops, nlay=1) - ncpl = gridprops["ncpl"] - assert ( - ncpl == answer_ncpl - ), f"Number of cells should be {answer_ncpl}. Found {ncpl}" - - if plot: - import matplotlib.pyplot as plt - - fig = plt.figure(figsize=(10, 10)) - ax = fig.add_subplot() - ax.set_aspect("equal") - voronoi_grid.plot(ax=ax) - plt.savefig(os.path.join(tpth, f"{name}.png")) - - return - - -if __name__ == "__main__": - # test_unstructured_grid_shell() - # test_unstructured_grid_dimensions() - # test_unstructured_minimal_grid() - # test_unstructured_complete_grid() - # test_loading_argus_meshes() - # test_create_unstructured_grid_from_verts() - # test_triangle_unstructured_grid() - # test_voronoi_vertex_grid() - test_voronoi_grid0(plot=True) - test_voronoi_grid1(plot=True) - test_voronoi_grid2(plot=True) - test_voronoi_grid3(plot=True) - test_voronoi_grid4(plot=True) - test_voronoi_grid5(plot=True)