From 7dec7c52db7c7bf3f8bca61de4d4a953ac1317d2 Mon Sep 17 00:00:00 2001 From: langevin-usgs Date: Fri, 26 Jul 2024 14:50:26 -0500 Subject: [PATCH] feat(lgr-disv): add to_disv_gridprops() method to lgr object (#2271) * feat(lgr-disv): add new to_disv_gridprops to lgr object * ruff ruff * modflow6 lgr tutorial was mistakenly marked as modflow-lgr tutorial * add test * add layer information to gridprops and update test * ruff * minor tweaks to require idomain the same for all layers * add deprecation warnings for gridlist_to_disv_gridprops * ruff * ruff --- ...gr_tutorial01.py => mf6_lgr_tutorial01.py} | 0 .docs/tutorials.rst | 10 +- autotest/test_export.py | 89 ++-- autotest/test_grid.py | 40 +- autotest/test_lgrutil.py | 74 ++- flopy/utils/cvfdutil.py | 13 + flopy/utils/lgrutil.py | 425 +++++++++++++++++- 7 files changed, 577 insertions(+), 74 deletions(-) rename .docs/Notebooks/{lgr_tutorial01.py => mf6_lgr_tutorial01.py} (100%) diff --git a/.docs/Notebooks/lgr_tutorial01.py b/.docs/Notebooks/mf6_lgr_tutorial01.py similarity index 100% rename from .docs/Notebooks/lgr_tutorial01.py rename to .docs/Notebooks/mf6_lgr_tutorial01.py diff --git a/.docs/tutorials.rst b/.docs/tutorials.rst index b1aaf3c17..7b2525504 100644 --- a/.docs/tutorials.rst +++ b/.docs/tutorials.rst @@ -35,6 +35,7 @@ MODFLOW 6 Notebooks/mf6_output_tutorial01 Notebooks/mf6_sfr_tutorial01 Notebooks/mf6_tutorial01 + Notebooks/mf6_lgr_tutorial01 MODFLOW-2005 @@ -50,15 +51,6 @@ MODFLOW-2005 Notebooks/mf_tutorial02 -MODFLOW-LGR ------------ - -.. toctree:: - :maxdepth: 2 - - Notebooks/lgr_tutorial01 - - MODFLOW-NWT ----------- diff --git a/autotest/test_export.py b/autotest/test_export.py index 4d9b9f379..0028ada34 100644 --- a/autotest/test_export.py +++ b/autotest/test_export.py @@ -1718,50 +1718,51 @@ def test_vtk_export_disv1_model(function_tmpdir): idomain=np.ones((nlay, nrow, ncol)), ) - from flopy.utils.cvfdutil import gridlist_to_disv_gridprops - - gridprops = gridlist_to_disv_gridprops([mg]) - gridprops["top"] = 0 - gridprops["botm"] = np.zeros((nlay, nrow * ncol), dtype=float) - 1 - gridprops["nlay"] = nlay - - disv = ModflowGwfdisv(gwf, **gridprops) - ic = ModflowGwfic(gwf, strt=10) - npf = ModflowGwfnpf(gwf) - - # Export model without specifying packages_names parameter - # create the vtk output - gwf = sim.get_model() - vtkobj = Vtk(gwf, binary=False) - vtkobj.add_model(gwf) - f = function_tmpdir / "gwf.vtk" - vtkobj.write(f) - - # load the output using the vtk standard library - gridreader = vtkUnstructuredGridReader() - gridreader.SetFileName(str(f)) - gridreader.Update() - grid = gridreader.GetOutput() - - # get the points - vtk_points = grid.GetPoints() - vtk_points = vtk_points.GetData() - vtk_points = vtk_to_numpy(vtk_points) - - # get cell locations (ia format of point to cell relationship) - cell_locations = vtk_to_numpy(grid.GetCellLocationsArray()) - cell_locations_answer = np.array([0, 8, 16, 24, 32, 40, 48, 56, 64]) - print(f"Found cell locations {cell_locations} in vtk file.") - print(f"Expecting cell locations {cell_locations_answer}") - errmsg = "vtk cell locations do not match expected result." - assert np.allclose(cell_locations, cell_locations_answer), errmsg - - cell_types = vtk_to_numpy(grid.GetCellTypesArray()) - cell_types_answer = np.array(9 * [42]) - print(f"Found cell types {cell_types} in vtk file.") - print(f"Expecting cell types {cell_types_answer}") - errmsg = "vtk cell types do not match expected result." - assert np.allclose(cell_types, cell_types_answer), errmsg + with pytest.deprecated_call(): + from flopy.utils.cvfdutil import gridlist_to_disv_gridprops + + gridprops = gridlist_to_disv_gridprops([mg]) + gridprops["top"] = 0 + gridprops["botm"] = np.zeros((nlay, nrow * ncol), dtype=float) - 1 + gridprops["nlay"] = nlay + + disv = ModflowGwfdisv(gwf, **gridprops) + ic = ModflowGwfic(gwf, strt=10) + npf = ModflowGwfnpf(gwf) + + # Export model without specifying packages_names parameter + # create the vtk output + gwf = sim.get_model() + vtkobj = Vtk(gwf, binary=False) + vtkobj.add_model(gwf) + f = function_tmpdir / "gwf.vtk" + vtkobj.write(f) + + # load the output using the vtk standard library + gridreader = vtkUnstructuredGridReader() + gridreader.SetFileName(str(f)) + gridreader.Update() + grid = gridreader.GetOutput() + + # get the points + vtk_points = grid.GetPoints() + vtk_points = vtk_points.GetData() + vtk_points = vtk_to_numpy(vtk_points) + + # get cell locations (ia format of point to cell relationship) + cell_locations = vtk_to_numpy(grid.GetCellLocationsArray()) + cell_locations_answer = np.array([0, 8, 16, 24, 32, 40, 48, 56, 64]) + print(f"Found cell locations {cell_locations} in vtk file.") + print(f"Expecting cell locations {cell_locations_answer}") + errmsg = "vtk cell locations do not match expected result." + assert np.allclose(cell_locations, cell_locations_answer), errmsg + + cell_types = vtk_to_numpy(grid.GetCellTypesArray()) + cell_types_answer = np.array(9 * [42]) + print(f"Found cell types {cell_types} in vtk file.") + print(f"Expecting cell types {cell_types_answer}") + errmsg = "vtk cell types do not match expected result." + assert np.allclose(cell_types, cell_types_answer), errmsg @pytest.mark.mf6 diff --git a/autotest/test_grid.py b/autotest/test_grid.py index f716d20bc..142fed140 100644 --- a/autotest/test_grid.py +++ b/autotest/test_grid.py @@ -943,25 +943,27 @@ def test_tocvfd3(): yoff=200, idomain=idomain, ) - gridprops = gridlist_to_disv_gridprops([sg1, sg2]) - assert "ncpl" in gridprops - assert "nvert" in gridprops - assert "vertices" in gridprops - assert "cell2d" in gridprops - - ncpl = gridprops["ncpl"] - nvert = gridprops["nvert"] - vertices = gridprops["vertices"] - cell2d = gridprops["cell2d"] - assert ncpl == 121 - assert nvert == 148 - assert len(vertices) == nvert - assert len(cell2d) == 121 - - # spot check information for cell 28 (zero based) - answer = [28, 250.0, 150.0, 7, 38, 142, 143, 45, 46, 44, 38] - for i, j in zip(cell2d[28], answer): - assert i == j, f"{i} not equal {j}" + + with pytest.deprecated_call(): + gridprops = gridlist_to_disv_gridprops([sg1, sg2]) + assert "ncpl" in gridprops + assert "nvert" in gridprops + assert "vertices" in gridprops + assert "cell2d" in gridprops + + ncpl = gridprops["ncpl"] + nvert = gridprops["nvert"] + vertices = gridprops["vertices"] + cell2d = gridprops["cell2d"] + assert ncpl == 121 + assert nvert == 148 + assert len(vertices) == nvert + assert len(cell2d) == 121 + + # spot check information for cell 28 (zero based) + answer = [28, 250.0, 150.0, 7, 38, 142, 143, 45, 46, 44, 38] + for i, j in zip(cell2d[28], answer): + assert i == j, f"{i} not equal {j}" @requires_pkg("shapely") diff --git a/autotest/test_lgrutil.py b/autotest/test_lgrutil.py index 439e4d5e6..67d7a5325 100644 --- a/autotest/test_lgrutil.py +++ b/autotest/test_lgrutil.py @@ -1,6 +1,6 @@ import numpy as np -from flopy.utils.lgrutil import Lgr +from flopy.utils.lgrutil import Lgr, LgrToDisv def test_lgrutil(): @@ -155,3 +155,75 @@ def test_lgrutil2(): ] assert np.allclose(lgr.delr, answer), f"{lgr.delr} /= {answer}" assert np.allclose(lgr.delc, answer), f"{lgr.delc} /= {answer}" + + +def test_lgrutil3(): + # Define parent grid information + xoffp = 0.0 + yoffp = 0.0 + nlayp = 3 + nrowp = 3 + ncolp = 3 + + dx = 100.0 + dy = 100.0 + dz = 10.0 + delrp = dx * np.ones(ncolp) + delcp = dy * np.ones(nrowp) + topp = dz * np.ones((nrowp, ncolp), dtype=float) + botmp = np.empty((nlayp, nrowp, ncolp), dtype=float) + for k in range(nlayp): + botmp[k] = -(k + 1) * dz + idomainp = np.ones((nlayp, nrowp, ncolp), dtype=int) + idomainp[:, nrowp // 2, ncolp // 2] = 0 + ncpp = 3 + ncppl = nlayp * [1] + lgr = Lgr( + nlayp, + nrowp, + ncolp, + delrp, + delcp, + topp, + botmp, + idomainp, + ncpp=ncpp, + ncppl=ncppl, + xllp=xoffp, + yllp=yoffp, + ) + + # check to make sure gridprops is accessible from lgr + gridprops = lgr.to_disv_gridprops() + assert "ncpl" in gridprops + assert "nvert" in gridprops + assert "vertices" in gridprops + assert "nlay" in gridprops + assert "top" in gridprops + assert "botm" in gridprops + assert gridprops["ncpl"] == 17 + assert gridprops["nvert"] == 32 + assert gridprops["nlay"] == 3 + + # test the lgr to disv class + lgrtodisv = LgrToDisv(lgr) + + # test guts of LgrToDisv to make sure hanging vertices added correctly + assert lgrtodisv.right_face_hanging[(1, 0)] == [0, 4, 8, 12] + assert lgrtodisv.left_face_hanging[(1, 2)] == [3, 7, 11, 15] + assert lgrtodisv.back_face_hanging[(2, 1)] == [12, 13, 14, 15] + assert lgrtodisv.front_face_hanging[(0, 1)] == [0, 1, 2, 3] + + assert lgrtodisv.iverts[1] == [1, 2, 6, 18, 17, 5] + assert lgrtodisv.iverts[3] == [4, 5, 20, 24, 9, 8] + assert lgrtodisv.iverts[4] == [6, 7, 11, 10, 27, 23] + assert lgrtodisv.iverts[6] == [9, 29, 30, 10, 14, 13] + + assert np.allclose(gridprops["top"], dz * np.ones((17,))) + + assert gridprops["botm"].shape == (3, 17) + b = np.empty((3, 17)) + b[0] = -dz + b[1] = -2 * dz + b[2] = -3 * dz + assert np.allclose(gridprops["botm"], b) diff --git a/flopy/utils/cvfdutil.py b/flopy/utils/cvfdutil.py index fea3b5602..3a59031d2 100644 --- a/flopy/utils/cvfdutil.py +++ b/flopy/utils/cvfdutil.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import pandas as pd @@ -390,6 +392,10 @@ def gridlist_to_disv_gridprops(gridlist): be numbered according to consecutive numbering of active cells in the grid list. + This function is deprecated in 3.8 and will be removed in 3.9. Use the + functionality in flopy.utils.cvfdutil.Lgr() to create a DISV mesh for a + nested grid. + Parameters ---------- gridlist : list @@ -403,6 +409,13 @@ def gridlist_to_disv_gridprops(gridlist): modflow6 disv package. """ + warnings.warn( + "the gridlist_to_disv_gridprops function is deprecated and will be " + "removed in version 3.9. Use flopy.utils.cvfdutil.Lgr() instead, which " + "allows a nested grid to be created and exported to a DISV mesh.", + PendingDeprecationWarning, + ) + verts, iverts = gridlist_to_verts(gridlist) gridprops = get_disv_gridprops(verts, iverts) return gridprops diff --git a/flopy/utils/lgrutil.py b/flopy/utils/lgrutil.py index 043b1a999..467f46a70 100644 --- a/flopy/utils/lgrutil.py +++ b/flopy/utils/lgrutil.py @@ -2,6 +2,7 @@ from ..discretization import StructuredGrid from ..modflow import Modflow +from .cvfdutil import get_disv_gridprops from .util_array import Util2d, Util3d @@ -162,7 +163,7 @@ def __init__( assert idomainp.shape == (nlayp, nrowp, ncolp) self.idomain = idomainp idxl, idxr, idxc = np.asarray(idomainp == 0).nonzero() - assert idxl.shape[0] > 1, "no zero values found in idomain" + assert idxl.shape[0] > 0, "no zero values found in idomain" # child cells per parent and child cells per parent layer self.ncpp = ncpp @@ -585,3 +586,425 @@ def child(self): yorigin, ) return simple_regular_grid + + def to_disv_gridprops(self): + """ + Create and return a gridprops dictionary that can be + used to create a disv grid (instead of a separate parent + and child representation). The gridprops dictionary can + be unpacked into the flopy.mf6.Modflowdisv() constructor + and flopy.discretization.VertexGrid() contructor. + + Note that export capability will only work if the parent + and child models have corresponding layers. + + Returns + ------- + gridprops : dict + Dictionary containing ncpl, nvert, vertices, cell2d, + nlay, top, and botm + + """ + return LgrToDisv(self).get_disv_gridprops() + + +class LgrToDisv: + def __init__(self, lgr): + """ + Helper class used to convert and Lgr() object into + the grid properties needed to create a disv vertex + nested grid. After instantiation, self.verts and + self.iverts are available. + + The primary work of this class is identify hanging + vertices along the shared parent-child boundary and + include these hanging vertices in the vertex indicence + list for parent cells. + + Parameters + ---------- + lgr : Lgr instance + Lgr() object describing a parent-child relation + + """ + + # store information + self.lgr = lgr + self.pgrid = lgr.parent.modelgrid + self.cgrid = lgr.child.modelgrid + + # count active parent and child cells + self.ncpl_parent = np.count_nonzero(self.pgrid.idomain[0] > 0) + self.ncpl_child = np.count_nonzero(self.cgrid.idomain[0] > 0) + self.ncpl = self.ncpl_child + self.ncpl_parent + + # find child vertices that act as hanging vertices on parent + # model cells + self.right_face_hanging = None + self.left_face_hanging = None + self.front_face_hanging = None + self.back_face_hanging = None + self.parent_ij_to_global = None + self.child_ij_to_global = None + self.find_hanging_vertices() + + # build global verts and iverts keeping only idomain > 0 + self.verts = None + self.iverts = None + self.build_verts_iverts() + + # todo: remove unused vertices? + + def find_hanging_vertices(self): + """ + Hanging vertices are vertices that must be included + along the edge of parent cells. These hanging vertices + mark the locations of corners of adjacent child cells. + Hanging vertices are not strictly + necessary to define the shape of a parent cell, but they are + required by modflow to describe connections between + parent and child cells. + + This routine finds hanging vertices parent cells along + a parent-child boundary. These hanging vertices are + stored in 4 member dictionaries, called right_face_hanging, + left_face_hanging, front_face_hanging, and back_face_hanging. + These dictionaries are used subsequently to insert + hanging vertices into the iverts array. + + """ + + # create dictionaries for parent left, right, back, and front + # faces that have a key that is parent (row, col) + # and a value that is a list of child vertex numbers + + # this list of child vertex numbers will be ordered from + # left to right (back/front) and from back to front (left/right) + # so when they are used later, two of them will need to be + # reversed so that clockwise ordering is maintained + + nrowc = self.lgr.nrow + ncolc = self.lgr.ncol + iverts = self.cgrid.iverts + cidomain = self.lgr.get_idomain() + + self.right_face_hanging = {} + self.left_face_hanging = {} + self.front_face_hanging = {} + self.back_face_hanging = {} + + # map (i, j) to global cell number + self.parent_ij_to_global = {} + self.child_ij_to_global = {} + + kc = 0 + nodec = 0 + for ic in range(nrowc): + for jc in range(ncolc): + plist = self.lgr.get_parent_connections(kc, ic, jc) + for (kp, ip, jp), idir in plist: + if cidomain[kc, ic, jc] == 0: + continue + + if ( + idir == -1 + ): # left child face connected to right parent face + # child vertices 0 and 3 added as hanging nodes + if (ip, jp) in self.right_face_hanging: + hlist = self.right_face_hanging.pop((ip, jp)) + else: + hlist = [] + ivlist = iverts[nodec] + for iv in (ivlist[0], ivlist[3]): + if iv not in hlist: + hlist.append(iv) + self.right_face_hanging[(ip, jp)] = hlist + + elif idir == 1: + # child vertices 1 and 2 added as hanging nodes + if (ip, jp) in self.left_face_hanging: + hlist = self.left_face_hanging.pop((ip, jp)) + else: + hlist = [] + ivlist = iverts[nodec] + for iv in (ivlist[1], ivlist[2]): + if iv not in hlist: + hlist.append(iv) + self.left_face_hanging[(ip, jp)] = hlist + + elif idir == 2: + # child vertices 0 and 1 added as hanging nodes + if (ip, jp) in self.front_face_hanging: + hlist = self.front_face_hanging.pop((ip, jp)) + else: + hlist = [] + ivlist = iverts[nodec] + for iv in (ivlist[0], ivlist[1]): + if iv not in hlist: + hlist.append(iv) + self.front_face_hanging[(ip, jp)] = hlist + + elif idir == -2: + # child vertices 3 and 2 added as hanging nodes + if (ip, jp) in self.back_face_hanging: + hlist = self.back_face_hanging.pop((ip, jp)) + else: + hlist = [] + ivlist = iverts[nodec] + for iv in (ivlist[3], ivlist[2]): + if iv not in hlist: + hlist.append(iv) + self.back_face_hanging[(ip, jp)] = hlist + + nodec += 1 + + def build_verts_iverts(self): + """ + Build the verts and iverts members. self.verts is a 2d + numpy array of size (nvert, 2). Column 1 is x and column 2 + is y. self.iverts is a list of size ncpl (number of cells + per layer) with each entry being the list of vertex indices + that define the cell. + + """ + + # stack vertex arrays; these will have more points than necessary, + # because parent and child vertices will overlap at corners, but + # duplicate vertices will be filtered later + pverts = self.pgrid.verts + cverts = self.cgrid.verts + nverts_parent = pverts.shape[0] + nverts_child = cverts.shape[0] + verts = np.vstack((pverts, cverts)) + + # build iverts list first with active parent cells + iverts = [] + iglo = 0 + for i in range(self.pgrid.nrow): + for j in range(self.pgrid.ncol): + if self.pgrid.idomain[0, i, j] > 0: + ivlist = self.pgrid._build_structured_iverts(i, j) + + # merge hanging vertices if they exist + ivlist = self.merge_hanging_vertices(i, j, ivlist) + + iverts.append(ivlist) + self.parent_ij_to_global[(i, j)] = iglo + iglo += 1 + + # now add active child cells + for i in range(self.cgrid.nrow): + for j in range(self.cgrid.ncol): + if self.cgrid.idomain[0, i, j] > 0: + ivlist = [ + iv + nverts_parent + for iv in self.cgrid._build_structured_iverts(i, j) + ] + iverts.append(ivlist) + self.child_ij_to_global[(i, j)] = iglo + iglo += 1 + self.verts = verts + self.iverts = iverts + + def merge_hanging_vertices(self, ip, jp, ivlist): + """ + Given a list of vertices (ivlist) for parent row and column + (ip, jp) merge hanging vertices from adjacent child cells + into ivlist. + + Parameters + ---------- + ip : int + parent cell row number + + jp : int + parent cell column number + + ivlist : list of ints + list of vertex indices that define the parent + cell (ip, jp) + + Returns + ------- + ivlist : list of ints + modified list of vertices that now also contains + any hanging vertices needed to properly define + a parent cell adjacent to child cells + + """ + assert len(ivlist) == 4 + child_ivlist_offset = self.pgrid.verts.shape[0] + + # construct back edge + idx = 0 + reverse = False + face_hanging = self.back_face_hanging + back_edge = [ivlist[idx]] + if (ip, jp) in face_hanging: + hlist = face_hanging[(ip, jp)] + if len(hlist) > 2: + hlist = hlist[1:-1] # do not include two ends + hlist = [h + child_ivlist_offset for h in hlist] + if reverse: + hlist = hlist[::-1] + else: + hlist = [] + back_edge = [ivlist[idx]] + hlist + + # construct right edge + idx = 1 + reverse = False + face_hanging = self.right_face_hanging + right_edge = [ivlist[idx]] + if (ip, jp) in face_hanging: + hlist = face_hanging[(ip, jp)] + if len(hlist) > 2: + hlist = hlist[1:-1] # do not include two ends + hlist = [h + child_ivlist_offset for h in hlist] + if reverse: + hlist = hlist[::-1] + else: + hlist = [] + right_edge = [ivlist[idx]] + hlist + + # construct front edge + idx = 2 + reverse = True + face_hanging = self.front_face_hanging + front_edge = [ivlist[idx]] + if (ip, jp) in face_hanging: + hlist = face_hanging[(ip, jp)] + if len(hlist) > 2: + hlist = hlist[1:-1] # do not include two ends + hlist = [h + child_ivlist_offset for h in hlist] + if reverse: + hlist = hlist[::-1] + else: + hlist = [] + front_edge = [ivlist[idx]] + hlist + + # construct left edge + idx = 3 + reverse = True + face_hanging = self.left_face_hanging + left_edge = [ivlist[idx]] + if (ip, jp) in face_hanging: + hlist = face_hanging[(ip, jp)] + if len(hlist) > 2: + hlist = hlist[1:-1] # do not include two ends + hlist = [h + child_ivlist_offset for h in hlist] + if reverse: + hlist = hlist[::-1] + else: + hlist = [] + left_edge = [ivlist[idx]] + hlist + + ivlist = back_edge + right_edge + front_edge + left_edge + + return ivlist + + def get_xcyc(self): + """ + Construct a 2d array of size (nvert, 2) that + contains the cell centers. + + Returns + ------- + xcyc : ndarray + 2d array of x, y positions for cell centers + + """ + xcyc = np.empty((self.ncpl, 2)) + pidx = self.pgrid.idomain[0] > 0 + cidx = self.cgrid.idomain[0] > 0 + px = self.pgrid.xcellcenters[pidx].flatten() + cx = self.cgrid.xcellcenters[cidx].flatten() + xcyc[:, 0] = np.vstack( + (np.atleast_2d(px).T, np.atleast_2d(cx).T) + ).flatten() + py = self.pgrid.ycellcenters[pidx].flatten() + cy = self.cgrid.ycellcenters[cidx].flatten() + xcyc[:, 1] = np.vstack( + (np.atleast_2d(py).T, np.atleast_2d(cy).T) + ).flatten() + return xcyc + + def get_top(self): + """ + Construct a 1d array of size (ncpl) that + contains the cell tops. + + Returns + ------- + top : ndarray + 1d array of top elevations + + """ + top = np.empty((self.ncpl,)) + pidx = self.pgrid.idomain[0] > 0 + cidx = self.cgrid.idomain[0] > 0 + pa = self.pgrid.top[pidx].flatten() + ca = self.cgrid.top[cidx].flatten() + top[:] = np.hstack((pa, ca)) + return top + + def get_botm(self): + """ + Construct a 2d array of size (nlay, ncpl) that + contains the cell bottoms. + + Returns + ------- + botm : ndarray + 2d array of bottom elevations + + """ + botm = np.empty((self.lgr.nlay, self.ncpl)) + pidx = self.pgrid.idomain[0] > 0 + cidx = self.cgrid.idomain[0] > 0 + for k in range(self.lgr.nlay): + pa = self.pgrid.botm[k, pidx].flatten() + ca = self.cgrid.botm[k, cidx].flatten() + botm[k, :] = np.hstack((pa, ca)) + return botm + + def get_disv_gridprops(self): + """ + Create and return a gridprops dictionary that can be + used to create a disv grid (instead of a separate parent + and child representation). The gridprops dictionary can + be unpacked into the flopy.mf6.Modflowdisv() constructor + and flopy.discretization.VertexGrid() contructor. + + Note that export capability will only work if the parent + and child models have corresponding layers. + + Returns + ------- + gridprops : dict + Dictionary containing ncpl, nvert, vertices, cell2d, + nlay, top, and botm + + """ + + # check + assert ( + self.lgr.ncppl.min() == self.lgr.ncppl.max() + ), "Exporting disv grid properties requires ncppl to be 1." + assert ( + self.lgr.nlayp == self.lgr.nlay + ), "Exporting disv grid properties requires parent and child models to have the same number of layers." + for k in range(self.lgr.nlayp - 1): + assert np.allclose( + self.lgr.idomain[k], self.lgr.idomain[k + 1] + ), "Exporting disv grid properties requires parent idomain is same for all layers." + + # get information and build gridprops + xcyc = self.get_xcyc() + top = self.get_top() + botm = self.get_botm() + gridprops = get_disv_gridprops(self.verts, self.iverts, xcyc=xcyc) + gridprops["nlay"] = self.lgr.nlay + gridprops["top"] = top + gridprops["botm"] = botm + return gridprops