diff --git a/autotest/test_particledata.py b/autotest/test_particledata.py index 29928d184..ce47d8cc7 100644 --- a/autotest/test_particledata.py +++ b/autotest/test_particledata.py @@ -2,7 +2,7 @@ import pytest from autotest.test_grid_cases import GridCases -from flopy.modpath import ParticleData +from flopy.modpath import FaceDataType, LRCParticleData, ParticleData structured_cells = [(0, 1, 1), (0, 1, 2)] structured_dtype = np.dtype( @@ -103,9 +103,60 @@ def test_particledata_structured_get_coord_list_no_ids(): ) -@pytest.mark.skip(reason="todo") +def test_lrcparticledata_get_coord_list_top_and_bottom_faces(): + rd = 1 + cd = 1 + sddata = FaceDataType( + rowdivisions5=rd, + columndivisions5=cd, + rowdivisions6=rd, + columndivisions6=cd, + ) + lrcregions = [[0, 1, 1, 0, 1, 1]] + data = LRCParticleData(subdivisiondata=sddata, lrcregions=lrcregions) + grid = GridCases().structured_small() + coords = data.get_coord_list(grid) + + num_cells = len( + [ + (lrc[3] - lrc[0]) * (lrc[4] - lrc[1]) * (lrc[5] - lrc[2]) + for lrc in lrcregions + ] + ) + assert ( + len(coords) == num_cells * rd * cd * 2 + ) # 1 particle each centered on the top and bottom faces + + def test_lrcparticledata_get_coord_list(): - pass + sddata = FaceDataType( + horizontaldivisions1=1, + verticaldivisions1=1, + horizontaldivisions2=1, + verticaldivisions2=1, + horizontaldivisions3=1, + verticaldivisions3=1, + horizontaldivisions4=1, + verticaldivisions4=1, + rowdivisions5=1, + columndivisions5=1, + rowdivisions6=1, + columndivisions6=1, + ) + lrcregions = [[0, 1, 1, 0, 1, 1]] + data = LRCParticleData(subdivisiondata=sddata, lrcregions=lrcregions) + grid = GridCases().structured_small() + coords = data.get_coord_list(grid) + + num_cells = len( + [ + (lrc[3] - lrc[0]) * (lrc[4] - lrc[1]) * (lrc[5] - lrc[2]) + for lrc in lrcregions + ] + ) + assert ( + len(coords) == num_cells * 6 + ) # 1 particle each centered on each face @pytest.mark.skip(reason="todo") diff --git a/flopy/modpath/mp7particledata.py b/flopy/modpath/mp7particledata.py index 95a9b204e..02acb3860 100644 --- a/flopy/modpath/mp7particledata.py +++ b/flopy/modpath/mp7particledata.py @@ -4,6 +4,7 @@ """ +from itertools import product import numpy as np from numpy.lib.recfunctions import unstructured_to_structured @@ -594,7 +595,6 @@ def __init__( self.columndivisions5 = columndivisions5 self.rowdivisions6 = rowdivisions6 self.columndivisions6 = columndivisions6 - return def write(self, f=None): """ @@ -738,7 +738,7 @@ class LRCParticleData: CellDataTypes that are used to create one or more particle templates in a particle group. If subdivisiondata is None, a default CellDataType with 27 particles per cell will be created (default is None). - lrcregions : list of lists tuples or np.ndarrays + lrcregions : list of lists of tuples or np.ndarrays Layer, row, column (zero-based) regions with particles created using the specified template parameters. A region is defined as a list/tuple of minlayer, minrow, mincolumn, maxlayer, maxrow, maxcolumn values. @@ -751,7 +751,7 @@ class LRCParticleData: -------- >>> import flopy - >>> pg = flopy.modpath.LRCParticleData(lrcregions=[0, 0, 0, 3, 10, 10]) + >>> pg = flopy.modpath.LRCParticleData(lrcregions=[[0, 0, 0, 3, 10, 10]]) """ @@ -792,8 +792,8 @@ def __init__(self, subdivisiondata=None, lrcregions=None): "a list with lists, tuples, or arrays".format(self.name) ) t = [] - for lrcregion in lrcregions: - t.append(np.array(lrcregion, dtype=np.int32)) + for region in lrcregions: + t.append(np.array(region, dtype=np.int32)) lrcregions = t else: raise TypeError( @@ -810,11 +810,11 @@ def __init__(self, subdivisiondata=None, lrcregions=None): ) # validate that there are 6 columns in each lrcregions entry - for idx, lrcregion in enumerate(lrcregions): - shapel = lrcregion.shape + for idx, region in enumerate(lrcregions): + shapel = np.array(region).shape if len(shapel) == 1: - lrcregions[idx] = lrcregion.reshape(1, shapel) - shapel = lrcregion[idx].shape + region = region.reshape(1, shapel[0]) + shapel = region.shape if shapel[1] != 6: raise ValueError( "{}: Each lrcregions entry must " @@ -822,17 +822,15 @@ def __init__(self, subdivisiondata=None, lrcregions=None): "{} columns".format(self.name, shapel[1]) ) - # totalcellregioncount = 0 - for lrcregion in lrcregions: - totalcellregioncount += lrcregion.shape[0] + for region in lrcregions: + totalcellregioncount += region.shape[0] # assign attributes self.particletemplatecount = shape self.totalcellregioncount = totalcellregioncount self.subdivisiondata = subdivisiondata self.lrcregions = lrcregions - return def write(self, f=None): """ @@ -876,7 +874,7 @@ def write(self, f=None): return - def to_list(self): + def get_coord_list(self, grid): """ Convert the particle representation to a list of global coordinates. @@ -884,7 +882,165 @@ def to_list(self): ------- A list of particle tuples (particle ID, x coord, y coord, z coord) """ - pass + + if grid.grid_type != "structured": + raise ValueError( + f"Particle representation is structured but grid is not" + ) + + particles = [] + + def add_particles(subdivisiondata, k, i, j): + cell_particles = [] + + if isinstance(subdivisiondata, FaceDataType): + # get cell face coords and span in each dimension + verts = grid.get_cell_vertices(i, j) + xs, ys = list(zip(*verts)) + minx, maxx = min(xs), max(xs) + miny, maxy = min(ys), max(ys) + minz, maxz = grid.botm[k, i, j], grid.top[i, j] + xspan = maxx - minx + yspan = maxy - miny + zspan = maxz - minz + + # face perpendicular to x (outer) + if ( + subdivisiondata.verticaldivisions1 > 0 + and subdivisiondata.horizontaldivisions1 > 0 + ): + yincr = yspan / subdivisiondata.horizontaldivisions1 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.horizontaldivisions1) + ] + zincr = zspan / subdivisiondata.verticaldivisions1 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions1) + ] + prod = list(product(*[ylocs, zlocs])) + cell_particles = cell_particles + [ + (minx, p[0], p[1]) for p in prod + ] + + # face perpendicular to x (outer) + if ( + subdivisiondata.verticaldivisions2 > 0 + and subdivisiondata.horizontaldivisions2 > 0 + ): + yincr = yspan / subdivisiondata.horizontaldivisions2 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * d)) + for d in range(subdivisiondata.horizontaldivisions2) + ] + zincr = zspan / subdivisiondata.verticaldivisions2 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions2) + ] + prod = list(product(*[ylocs, zlocs])) + cell_particles = cell_particles + [ + (maxx, p[0], p[1]) for p in prod + ] + + # face perpendicular to y (inner) + if ( + subdivisiondata.verticaldivisions3 > 0 + and subdivisiondata.horizontaldivisions3 > 0 + ): + xincr = xspan / subdivisiondata.horizontaldivisions3 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.horizontaldivisions3) + ] + zincr = zspan / subdivisiondata.verticaldivisions3 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions3) + ] + prod = list(product(*[xlocs, zlocs])) + cell_particles = cell_particles + [ + (p[0], miny, p[1]) for p in prod + ] + + # face perpendicular to y (outer) + if ( + subdivisiondata.verticaldivisions4 > 0 + and subdivisiondata.horizontaldivisions4 > 0 + ): + xincr = xspan / subdivisiondata.horizontaldivisions4 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.horizontaldivisions4) + ] + zincr = zspan / subdivisiondata.verticaldivisions4 + zlocs = [ + (minz + (zincr * 0.5) + (zincr * d)) + for d in range(subdivisiondata.verticaldivisions4) + ] + prod = list(product(*[xlocs, zlocs])) + cell_particles = cell_particles + [ + (p[0], maxy, p[1]) for p in prod + ] + + # bottom face + if ( + subdivisiondata.rowdivisions5 > 0 + and subdivisiondata.columndivisions5 > 0 + ): + xincr = xspan / subdivisiondata.columndivisions5 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columndivisions5) + ] + yincr = yspan / subdivisiondata.rowdivisions5 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * rd)) + for rd in range(subdivisiondata.rowdivisions5) + ] + prod = list(product(*[xlocs, ylocs])) + cell_particles = cell_particles + [ + (p[0], p[1], minz) for p in prod + ] + + # top face + if ( + subdivisiondata.rowdivisions6 > 0 + and subdivisiondata.columndivisions6 > 0 + ): + xincr = xspan / subdivisiondata.columndivisions6 + xlocs = [ + (minx + (xincr * 0.5) + (xincr * rd)) + for rd in range(subdivisiondata.columndivisions6) + ] + yincr = yspan / subdivisiondata.rowdivisions5 + ylocs = [ + (miny + (yincr * 0.5) + (yincr * rd)) + for rd in range(subdivisiondata.rowdivisions6) + ] + prod = list(product(*[xlocs, ylocs])) + cell_particles = cell_particles + [ + (p[0], p[1], maxz) for p in prod + ] + elif isinstance(subdivisiondata, CellDataType): + pass + else: + raise ValueError( + f"Unsupported subdivision data type: {type(self.subdivisiondata)}" + ) + + return cell_particles + + for region in self.lrcregions: + mink, mini, minj, maxk, maxi, maxj = region + for k in range(mink, maxk + 1): + for i in range(mini, maxi + 1): + for j in range(minj, maxj + 1): + for sd in self.subdivisiondata: + particles = particles + add_particles(sd, k, i, j) + + return particles class NodeParticleData: