Skip to content

Commit

Permalink
feat(LRCParticleData): implement get_coord_list() method with minimal…
Browse files Browse the repository at this point in the history 
… tests
  • Loading branch information
@wpbonelli
wpbonelli committed Mar 26, 2023
1 parent 324ac18 commit a78a7f8
Show file tree
Hide file tree
Showing 2 changed files with 225 additions and 18 deletions.
57 changes: 54 additions & 3 deletions autotest/test_particledata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down Expand Up @@ -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")
Expand Down
186 changes: 171 additions & 15 deletions flopy/modpath/mp7particledata.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
"""
from itertools import product

import numpy as np
from numpy.lib.recfunctions import unstructured_to_structured
Expand Down Expand Up @@ -594,7 +595,6 @@ def __init__(
self.columndivisions5 = columndivisions5
self.rowdivisions6 = rowdivisions6
self.columndivisions6 = columndivisions6
return

def write(self, f=None):
"""
Expand Down Expand Up @@ -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.
Expand All @@ -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]])
"""

Expand Down Expand Up @@ -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(
Expand All @@ -810,29 +810,27 @@ 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 "
"have 6 columns passed lrcregions has "
"{} 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):
"""
Expand Down Expand Up @@ -876,15 +874,173 @@ 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.
Returns
-------
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:
Expand Down

0 comments on commit a78a7f8

Please sign in to comment.