fix(get_structured_faceflows): cover edge cases, expand tests

[wpbonelli]

• close bug: get_structured_faceflows outputting wrong flows #1911
• fix edge cases for grids extended (>1 cell) in only 1 or 2 dimensions
• update mf6 freyberg test to compare mf6 and mf2005 face flows
• add nlay/nrow/ncol params (needed if grb_file not provided)

[codecov]

Codecov Report

Merging #1968 (65e0c03) into develop (6e23400) will increase coverage by 0.0%.
The diff coverage is 100.0%.

@@           Coverage Diff           @@
##           develop   #1968   +/-   ##
=======================================
  Coverage     72.6%   72.7%           
=======================================
  Files          257     257           
  Lines        57800   57809    +9     
=======================================
+ Hits         42017   42063   +46     
+ Misses       15783   15746   -37     

Files	Coverage Δ
flopy/mf6/utils/postprocessing.py	96.8% <100.0%> (-1.4%)	⬇️

... and 3 files with indirect coverage changes

[langevin-usgs]

Not sure if this would be any clearer, but scrolling through k, i, j seems a bit more intuitive to me. Just a thought. I didn't fully test, so there could be a sign issue with the flows. I think the convention is positive for leaving the k, i, j cell. Feel free to scrap if this is wrong or unclear.

# requires nlay, nrow, ncol, ia, ja, flowja be in scope
frf = np.zeros((nlay, nrow, ncol))
fff = np.zeros((nlay, nrow, ncol))
flf = np.zeros((nlay, nrow, ncol))

def get_node(k, i, j):
    return j + i * ncol + k * nrow * ncol

def get_flow(n, m, ia, ja, flowja):
    for ipos in range(ia[n] + 1, ia[n + 1]):
        if m == ja[ipos]:
            return flowja[ipos]
    return 0.
        
for k in range(nlay):
    for i in range(nrow):
        for j in range(ncol):
        
            # get node number for k, i, j
            n = get_node(k, i, j)
            
            # fill flow to right (positive to east)
            if j < ncol - 1:
                m = get_node(k, i, j + 1)
                frf[k, i, j] = -get_flow(n, m, ia, ja, flowja)
            
            # fill flow to front (positive to south)
            if i < nrow - 1:
                m = get_node(k, i + 1, j)
                fff[k, i, j] = -get_flow(n, m, ia, ja, flowja)

            # fill flow to lower (positive flow upward)
            if k < nlay - 1:
                m = get_node(k + 1, i, j)
                flf[k, i, j] = -get_flow(n, m, ia, ja, flowja)

[wpbonelli]

@langevin-usgs your proposed version passes all tests. I agree iterating cell indices (k, i, j) is easier to understand than nodes / CSR matrix entries.

I did a quick benchmark. The grid shape (nlay, nrow, ncol) is in the name column, all times are ms.

Name	Method	Mean	StdDev	Median	IQR
test_get_structured_faceflows_benchmark[1-1-1000]	indices	2.4508	0.0773	2.4539	0.0748
test_get_structured_faceflows_benchmark[1-1-1000]	nodes	1.9363	0.2961	1.8995	0.0644
test_get_structured_faceflows_benchmark[1-1000-1]	indices	2.5823	0.0696	2.5811	0.0625
test_get_structured_faceflows_benchmark[1-1000-1]	nodes	1.8413	0.2585	1.7936	0.0890
test_get_structured_faceflows_benchmark[1000-1-1]	indices	2.7662	0.1463	2.7433	0.0706
test_get_structured_faceflows_benchmark[1000-1-1]	nodes	1.8334	0.3278	1.7605	0.0767
test_get_structured_faceflows_benchmark[100-1-100]	indices	28.2745	0.4348	28.3100	0.5492
test_get_structured_faceflows_benchmark[100-1-100]	nodes	17.0061	1.7892	16.2677	0.5343
test_get_structured_faceflows_benchmark[1-100-100]	indices	28.2635	0.3137	28.2699	0.4019
test_get_structured_faceflows_benchmark[1-100-100]	nodes	16.7042	1.2910	16.4433	0.6004
test_get_structured_faceflows_benchmark[100-100-1]	indices	29.9212	0.2994	29.9362	0.4460
test_get_structured_faceflows_benchmark[100-100-1]	nodes	16.0866	0.6630	15.9132	0.3964
test_get_structured_faceflows_benchmark[50-50-50]	indices	532.4554	5.1744	529.2169	8.5408
test_get_structured_faceflows_benchmark[50-50-50]	nodes	262.1846	10.6833	259.3503	10.3721

Iterating nodes is generally faster than indices. I suspect the inner loop checking up to 6 neighbors to find the correct one in get_flow() is the culprit. Not sure what the right maintainability vs. speed tradeoff is here. Maybe there is a way to avoid that inner loop while keeping the more intuitive approach.

[langevin-usgs]

Hey @wpbonelli, thanks for running the tests. I'm not sure we can avoid the inner loop with this k,i,j indexing. If any of the adjacent cells for cell n are "idomained out" then they are not included in ja. So I think we need to explicitly check for their presence. There probably would be some way to speed this up, though I'm not sure speed is that critical for this particular routine, as it's use case is fairly limited. I don't have a strong preference for the indices or node approach, and am fine with whatever you decide to include. We have this discussion for the record if we ever need to revisit.