From 45c3a56bbd6bec1b4beb3e7aaa76e3631bd2ecfa Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D" Date: Sat, 9 Sep 2023 16:19:01 -0500 Subject: [PATCH 1/3] refactor(disv): improve cell area calculation * improve cell area calculation * close #1346 --- src/Model/GroundWaterFlow/gwf3disv8.f90 | 10 ++++++++-- src/Model/ModelUtilities/DisvGeom.f90 | 10 ++++++++-- 2 files changed, 16 insertions(+), 4 deletions(-) diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90 index 17895b55f38..25b1d7d1758 100644 --- a/src/Model/GroundWaterFlow/gwf3disv8.f90 +++ b/src/Model/GroundWaterFlow/gwf3disv8.f90 @@ -1325,13 +1325,19 @@ function get_cell2d_area(this, icell2d) result(area) integer(I4B) :: ivert integer(I4B) :: nvert integer(I4B) :: icount + integer(I4B) :: iv1 real(DP) :: x real(DP) :: y + real(DP) :: x1 + real(DP) :: y1 ! ------------------------------------------------------------------------------ ! area = DZERO nvert = this%iavert(icell2d + 1) - this%iavert(icell2d) icount = 1 + iv1 = this%javert(this%iavert(icell2d)) + x1 = this%vertices(1, iv1) + y1 = this%vertices(2, iv1) do ivert = this%iavert(icell2d), this%iavert(icell2d + 1) - 1 x = this%vertices(1, this%javert(ivert)) if (icount < nvert) then @@ -1339,7 +1345,7 @@ function get_cell2d_area(this, icell2d) result(area) else y = this%vertices(2, this%javert(this%iavert(icell2d))) end if - area = area + x * y + area = area + (x - x1) * (y - y1) icount = icount + 1 end do ! @@ -1351,7 +1357,7 @@ function get_cell2d_area(this, icell2d) result(area) else x = this%vertices(1, this%javert(this%iavert(icell2d))) end if - area = area - x * y + area = area - (x - x1) * (y - y1) icount = icount + 1 end do ! diff --git a/src/Model/ModelUtilities/DisvGeom.f90 b/src/Model/ModelUtilities/DisvGeom.f90 index 581ad424b88..d21a10b80ec 100644 --- a/src/Model/ModelUtilities/DisvGeom.f90 +++ b/src/Model/ModelUtilities/DisvGeom.f90 @@ -335,13 +335,19 @@ function get_area(this) result(area) integer(I4B) :: ivert integer(I4B) :: nvert integer(I4B) :: icount + integer(I4B) :: iv1 real(DP) :: x real(DP) :: y + real(DP) :: x1 + real(DP) :: y1 ! ------------------------------------------------------------------------------ ! area = DZERO nvert = this%iavert(this%j + 1) - this%iavert(this%j) icount = 1 + iv1 = this%javert(this%iavert(this%j)) + x1 = this%vertex_grid(1, iv1) + y1 = this%vertex_grid(2, iv1) do ivert = this%iavert(this%j), this%iavert(this%j + 1) - 1 x = this%vertex_grid(1, this%javert(ivert)) if (icount < nvert) then @@ -349,7 +355,7 @@ function get_area(this) result(area) else y = this%vertex_grid(2, this%javert(this%iavert(this%j))) end if - area = area + x * y + area = area + (x - x1) * (y - y1) icount = icount + 1 end do ! @@ -361,7 +367,7 @@ function get_area(this) result(area) else x = this%vertex_grid(1, this%javert(this%iavert(this%j))) end if - area = area - x * y + area = area - (x - x1) * (y - y1) icount = icount + 1 end do ! From 1ffd54846e8d57cef764c0b8381b31dfc73fbe31 Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D" Date: Mon, 18 Sep 2023 11:29:05 -0500 Subject: [PATCH 2/3] add release notes and disv test --- autotest/test_gwf_disv.py | 106 +++++++++++++++++++++++++++++++++++ doc/ReleaseNotes/develop.tex | 1 + 2 files changed, 107 insertions(+) create mode 100644 autotest/test_gwf_disv.py diff --git a/autotest/test_gwf_disv.py b/autotest/test_gwf_disv.py new file mode 100644 index 00000000000..69938066db3 --- /dev/null +++ b/autotest/test_gwf_disv.py @@ -0,0 +1,106 @@ +""" +Test of GWF DISV Package. Use the flopy disv tool to create +a simple regular grid example, but using DISV instead of DISU. +Use a large offset for x and y vertices to ensure that the area +calculation in MODFLOW 6 is correct. For the second case, set +one of the cells inactive and test to make sure connectivity +in binary grid file is correct. + +""" + +import os + +import flopy +import numpy as np +import pytest +from flopy.utils.gridutil import get_disv_kwargs + +from framework import TestFramework +from simulation import TestSimulation + +ex = ["disv01a", "disv01b"] + + +def build_model(idx, dir, mf6): + name = ex[idx] + ws = dir + nlay = 3 + nrow = 3 + ncol = 3 + delr = 10.0 + delc = 10.0 + top = 0 + botm = [-10, -20, -30] + xoff = 100000000.0 + yoff = 100000000.0 + disvkwargs = get_disv_kwargs( + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + xoff, + yoff, + ) + if idx == 1: + # for the second test, set one cell to idomain = 0 + idomain = np.ones((nlay, nrow * ncol), dtype=int) + idomain[0, 1] = 0 + disvkwargs["idomain"] = idomain + + sim = flopy.mf6.MFSimulation( + sim_name=name, + version="mf6", + exe_name=mf6, + sim_ws=ws, + ) + tdis = flopy.mf6.ModflowTdis(sim) + gwf = flopy.mf6.ModflowGwf(sim, modelname=name) + ims = flopy.mf6.ModflowIms(sim, print_option="SUMMARY") + disv = flopy.mf6.ModflowGwfdisv(gwf, **disvkwargs) + ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0) + npf = flopy.mf6.ModflowGwfnpf(gwf) + spd = {0: [[(0, 0), 1.0], [(0, nrow * ncol - 1), 0.0]]} + chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, stress_period_data=spd) + return sim, None + + +def eval_results(sim): + print("evaluating results...") + + name = sim.name + + fname = os.path.join(sim.simpath, name + ".disv.grb") + grbobj = flopy.mf6.utils.MfGrdFile(fname) + ncpl = grbobj._datadict["NCPL"] + ia = grbobj._datadict["IA"] + ja = grbobj._datadict["JA"] + + if sim.idxsim == 1: + # assert ncpl == disvkwargs["ncpl"] + assert np.array_equal(ia[0:4], np.array([1, 4, 4, 7])) + assert np.array_equal(ja[:6], np.array([1, 4, 10, 3, 6, 12])) + assert ia[-1] == 127 + assert ia.shape[0] == 28, "ia should have size of 28" + assert ja.shape[0] == 126, "ja should have size of 126" + + +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/doc/ReleaseNotes/develop.tex b/doc/ReleaseNotes/develop.tex index ce452f1ad43..e09827fb8bf 100644 --- a/doc/ReleaseNotes/develop.tex +++ b/doc/ReleaseNotes/develop.tex @@ -21,6 +21,7 @@ \underline{BASIC FUNCTIONALITY} \begin{itemize} \item Improve error message if the size of data read from a binary array file is inconsistent with READARRAY control line and variable description keywords. + \item The area calculation for cells in the DISV package was inaccurate for some cases with very large cell vertex coordinates. The area calculation was improved by using transformed cell vertex coordinates prior to making the area calculation. % \item xxx % \item xxx \end{itemize} From b6241ea6d3ec20b1592056d48b081ea5f50d880c Mon Sep 17 00:00:00 2001 From: "Langevin, Christian D" Date: Mon, 18 Sep 2023 13:39:32 -0500 Subject: [PATCH 3/3] more test updates --- autotest/test_gwf_disu.py | 127 ++++++++++++++++++---------------- autotest/test_gwf_disv.py | 2 +- autotest/test_gwf_disv_uzf.py | 122 ++++++++------------------------ 3 files changed, 96 insertions(+), 155 deletions(-) diff --git a/autotest/test_gwf_disu.py b/autotest/test_gwf_disu.py index 85880b95d82..ac21a31ca5d 100644 --- a/autotest/test_gwf_disu.py +++ b/autotest/test_gwf_disu.py @@ -1,13 +1,28 @@ +""" +Test of GWF DISU Package. Use the flopy disu tool to create +a simple regular grid example, but using DISU instead of DIS. +The first case is just a simple test. For the second case, set +one of the cells inactive and test to make sure connectivity +in binary grid file is correct. + +""" + import os import flopy import numpy as np +import pytest from flopy.utils.gridutil import get_disu_kwargs +from framework import TestFramework +from simulation import TestSimulation + +ex = ["disu01a", "disu01b"] + -def test_disu_simple(tmpdir, targets): - mf6 = targets["mf6"] - name = "disu01a" +def build_model(idx, dir, mf6): + name = ex[idx] + ws = dir nlay = 3 nrow = 3 ncol = 3 @@ -15,9 +30,26 @@ def test_disu_simple(tmpdir, targets): delc = 10.0 * np.ones(nrow) top = 0 botm = [-10, -20, -30] - disukwargs = get_disu_kwargs(nlay, nrow, ncol, delr, delc, top, botm) + disukwargs = get_disu_kwargs( + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + ) + if idx == 1: + # for the second test, set one cell to idomain = 0 + idomain = np.ones((nlay, nrow * ncol), dtype=int) + idomain[0, 1] = 0 + disukwargs["idomain"] = idomain + sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name=mf6, sim_ws=str(tmpdir) + sim_name=name, + version="mf6", + exe_name=mf6, + sim_ws=ws, ) tdis = flopy.mf6.ModflowTdis(sim) gwf = flopy.mf6.ModflowGwf(sim, modelname=name) @@ -25,67 +57,44 @@ def test_disu_simple(tmpdir, targets): disu = flopy.mf6.ModflowGwfdisu(gwf, **disukwargs) ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0) npf = flopy.mf6.ModflowGwfnpf(gwf) - spd = {0: [[(0,), 1.0], [(nrow * ncol - 1,), 0.0]]} + spd = {0: [[(0,), 1.0], [(nrow * ncol - 1), 0.0]]} chd = flopy.mf6.modflow.mfgwfchd.ModflowGwfchd(gwf, stress_period_data=spd) - sim.write_simulation() - sim.run_simulation() + return sim, None -def test_disu_idomain_simple(tmpdir, targets): - mf6 = targets["mf6"] - name = "disu01b" - nlay = 3 - nrow = 3 - ncol = 3 - delr = 10.0 * np.ones(ncol) - delc = 10.0 * np.ones(nrow) - top = 0 - botm = [-10, -20, -30] - idomain = np.ones(nlay * nrow * ncol, dtype=int) - idomain[1] = 0 - disukwargs = get_disu_kwargs(nlay, nrow, ncol, delr, delc, top, botm) - disukwargs["idomain"] = idomain - sim = flopy.mf6.MFSimulation( - sim_name=name, version="mf6", exe_name=mf6, sim_ws=str(tmpdir) - ) - tdis = flopy.mf6.ModflowTdis(sim) - gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True) - ims = flopy.mf6.ModflowIms(sim, print_option="SUMMARY") - disu = flopy.mf6.ModflowGwfdisu(gwf, **disukwargs) - ic = flopy.mf6.ModflowGwfic(gwf, strt=0.0) - npf = flopy.mf6.ModflowGwfnpf(gwf) - spd = {0: [[(0,), 1.0], [(nrow * ncol - 1,), 0.0]]} - chd = flopy.mf6.modflow.ModflowGwfchd(gwf, stress_period_data=spd) - oc = flopy.mf6.modflow.ModflowGwfoc( - gwf, - budget_filerecord=f"{name}.bud", - head_filerecord=f"{name}.hds", - saverecord=[("HEAD", "LAST"), ("BUDGET", "LAST")], - ) - sim.write_simulation() - sim.run_simulation() +def eval_results(sim): + print("evaluating results...") - # check binary grid file - fname = os.path.join(str(tmpdir), name + ".disu.grb") + name = sim.name + + fname = os.path.join(sim.simpath, name + ".disu.grb") grbobj = flopy.mf6.utils.MfGrdFile(fname) nodes = grbobj._datadict["NODES"] ia = grbobj._datadict["IA"] ja = grbobj._datadict["JA"] - assert nodes == disukwargs["nodes"] - assert np.array_equal(ia[0:4], np.array([1, 4, 4, 7])) - assert np.array_equal(ja[:6], np.array([1, 4, 10, 3, 6, 12])) - assert ia[-1] == 127 - assert ia.shape[0] == 28, "ia should have size of 28" - assert ja.shape[0] == 126, "ja should have size of 126" - # load head array and ensure nodata value in second cell - fname = os.path.join(str(tmpdir), name + ".hds") - hdsobj = flopy.utils.HeadFile(fname) - head = hdsobj.get_alldata().flatten() - assert head[1] == 1.0e30 + if sim.idxsim == 1: + assert np.array_equal(ia[0:4], np.array([1, 4, 4, 7])) + assert np.array_equal(ja[:6], np.array([1, 4, 10, 3, 6, 12])) + assert ia[-1] == 127 + assert ia.shape[0] == 28, "ia should have size of 28" + assert ja.shape[0] == 126, "ja should have size of 126" + - # load flowja to make sure it is the right size - fname = os.path.join(str(tmpdir), name + ".bud") - budobj = flopy.utils.CellBudgetFile(fname, precision="double") - flowja = budobj.get_data(text="FLOW-JA-FACE")[0].flatten() - assert flowja.shape[0] == 126 +@pytest.mark.parametrize("idx, name", list(enumerate(ex))) +def test_mf6model(idx, name, function_tmpdir, targets): + ws = str(function_tmpdir) + sim, _ = build_model(idx, ws, targets.mf6) + sim.write_simulation() + sim.run_simulation() + test = TestFramework() + test.run( + TestSimulation( + name=name, + exe_dict=targets, + exfunc=eval_results, + idxsim=idx, + make_comparison=False, + ), + ws, + ) diff --git a/autotest/test_gwf_disv.py b/autotest/test_gwf_disv.py index 69938066db3..9075ad5e749 100644 --- a/autotest/test_gwf_disv.py +++ b/autotest/test_gwf_disv.py @@ -1,6 +1,6 @@ """ Test of GWF DISV Package. Use the flopy disv tool to create -a simple regular grid example, but using DISV instead of DISU. +a simple regular grid example, but using DISV instead of DIS. Use a large offset for x and y vertices to ensure that the area calculation in MODFLOW 6 is correct. For the second case, set one of the cells inactive and test to make sure connectivity diff --git a/autotest/test_gwf_disv_uzf.py b/autotest/test_gwf_disv_uzf.py index cf89372267d..d331fd0cf0b 100644 --- a/autotest/test_gwf_disv_uzf.py +++ b/autotest/test_gwf_disv_uzf.py @@ -14,97 +14,38 @@ import flopy.utils.cvfdutil import numpy as np import pytest +from flopy.utils.gridutil import get_disv_kwargs + from framework import TestFramework from simulation import TestSimulation ex = ["disv_with_uzf"] nlay = 5 +nrow = 10 +ncol = 10 +ncpl = nrow * ncol +delr = 1. +delc = 1. nper = 5 perlen = [10] * 5 nstp = [5] * 5 tsmult = len(perlen) * [1.0] +top = 25. botm = [20.0, 15.0, 10.0, 5.0, 0.0] strt = 20 nouter, ninner = 100, 300 hclose, rclose, relax = 1e-9, 1e-3, 0.97 -ghb_ids = [] - - -def create_disv_mesh(): - # Create a grid of verts - nx, ny = (11, 11) - x = np.linspace(0, 10, nx) - y = np.linspace(0, 10, ny) - xv, yv = np.meshgrid(x, y) - yv = np.flipud(yv) - - verts = [] - vid = 0 - vert_lkup = {} - for i in yv[:, 0]: - for j in xv[0, :]: - vert_lkup.update({(float(j), float(i)): vid}) - verts.append([int(vid), float(j), float(i)]) - vid += 1 - - ivert = [] - ivid = 0 - xyverts = [] - xc, yc = [], [] # for storing the cell center location - for i in yv[:-1, 0]: - for j in xv[0, :-1]: - xlst, ylst = [], [] - vid_lst = [] - # Start with upper-left corner and go clockwise - for ct in [0, 1, 2, 3]: - if ct == 0: - iadj = 0.0 - jadj = 0.0 - elif ct == 1: - iadj = 0.0 - jadj = 1.0 - elif ct == 2: - iadj = -1.0 - jadj = 1.0 - elif ct == 3: - iadj = -1.0 - jadj = 0.0 - - vid = vert_lkup[(float(j + jadj), float(i + iadj))] - vid_lst.append(vid) - - xlst.append(float(j + jadj)) - ylst.append(float(i + iadj)) - - xc.append(np.mean(xlst)) - yc.append(np.mean(ylst)) - xyverts.append(list(zip(xlst, ylst))) - - rec = [ivid] + vid_lst - ivert.append(rec) - - # if ivert part of right boundary, store id - if j == 9.0: - ghb_ids.append(ivid) - - ivid += 1 - - # finally, create a cell2d record - cell2d = [] - for ix, iv in enumerate(ivert): - xvt, yvt = np.array(xyverts[ix]).T - if flopy.utils.geometry.is_clockwise(xvt, yvt): - rec = [iv[0], xc[ix], yc[ix], len(iv[1:])] + iv[1:] - else: - iiv = iv[1:][::-1] - rec = [iv[0], xc[ix], yc[ix], len(iiv)] + iiv - - cell2d.append(rec) - - return verts, cell2d - -verts, cell2d = create_disv_mesh() +# use flopy util to get disv arguments +disvkwargs = get_disv_kwargs( + nlay, + nrow, + ncol, + delr, + delc, + top, + botm, + ) # Work up UZF data iuzno = 0 @@ -117,7 +58,7 @@ def create_disv_mesh(): eps = 3.5 for k in np.arange(nlay): - for i in np.arange(0, len(cell2d), 1): + for i in np.arange(0, ncpl, 1): if k == 0: landflg = 1 surfdp = 0.25 @@ -128,7 +69,7 @@ def create_disv_mesh(): if k == nlay - 1: ivertcon = -1 else: - ivertcon = iuzno + len(cell2d) + ivertcon = iuzno + ncpl bndnm = "uzf" + "{0:03d}".format(int(i + 1)) uzf_pkdat.append( @@ -160,7 +101,7 @@ def create_disv_mesh(): spd = [] iuzno = 0 for k in np.arange(nlay): - for i in np.arange(0, len(cell2d), 1): + for i in np.arange(0, ncpl, 1): if k == 0: if t == 0: finf = 0.15 @@ -180,6 +121,7 @@ def create_disv_mesh(): # Work up the GHB boundary +ghb_ids = [(ncol - 1) + i * ncol for i in range(nrow)] ghb_spd = [] cond = 1e4 for k in np.arange(3, 5, 1): @@ -230,19 +172,9 @@ def build_model(idx, dir): ) sim.register_ims_package(ims, [gwf.name]) - ncpl = len(cell2d) - nvert = len(verts) - disv = flopy.mf6.ModflowGwfdisv( - gwf, - nlay=nlay, - ncpl=ncpl, - nvert=nvert, - top=25.0, - botm=botm, - vertices=verts, - cell2d=cell2d, - ) - + # disv + disv = flopy.mf6.ModflowGwfdisv(gwf, **disvkwargs) + # initial conditions ic = flopy.mf6.ModflowGwfic(gwf, strt=strt) @@ -389,7 +321,7 @@ def eval_model(sim): for rw in np.arange(arr.shape[0]): fullrw = arr[rw] for cl in np.arange(len(fullrw) - 1): - assert abs(fullrw[cl]) >= abs(fullrw[cl + 1]), ( + assert abs(fullrw[cl]) + 0.01 >= abs(fullrw[cl + 1]), ( "gwet not decreasing to the right as expected. Stress Period: " + str(tm + 1) + "; Row: " @@ -423,7 +355,7 @@ def eval_model(sim): for rw in np.arange(arr.shape[0]): fullrw = arr[rw] for cl in np.arange(len(fullrw) - 1): - assert abs(fullrw[cl]) <= abs(fullrw[cl + 1]), ( + assert abs(fullrw[cl]) <= abs(fullrw[cl + 1]) + 0.01, ( "gwet not decreasing to the right as expected. Stress Period: " + str(tm + 1) + "; Row: "