quad_utils_mod search: replace 3D scratch arrays with compressed sparse row format

[hkershaw-brown]

Description:

Quad_utils_mod box search for fully irregular grids:

• Determine the search storage requirement at run time from the grid. Currently there are 3 hardcoded cases, leading to people adding "big" numbers to the module to get their code running vs erroring out. The hardcoded cases I believe are based on the POP workhorse and high res grids ~2015 (looking at POP/model_mod.f90 blame)
• Replaces the quad_utils_mod search storage (3D arrays) for fully irregular grids with a compressed sparse row format.

The runtime calculation of storage requirements is needed in particular for regional CESM where people are creating their own regional grids, but is applicable to any model using quad_utils.
The per-code memory is quite high for the box storage (see #979 & https://github.com/hkershaw-brown/quad_search_perf/blob/main/README.md). This pull request does not address distributing memory for the search storage (may be needed later), but does reduce the per core memory.

Compressed Sparse Row format:
The GRID_QUAD_FULLY_IRREGULAR coarse-index search-info array setup previously allocated reg_list_lon/lat(nrx, nry, max_reg_list_num) as a 3D temporary, reaching ~5.5 GB at nrx=nry=900. It also errored out if any coarse box accumulated more than max_reg_list_num=800 candidate quads. Noted in issue #979

Replace with a two-pass compressed sparse row build:

• Pass 1: count overlaps into grid_num; peak memory ~11 MB
• Prefix sum: derive grid_start and exact total entry count
• Pass 2: fill flat lists using a cursor array

Replace hardcoded num_reg with calculation based on grid
num_reg = max(10, min(nx_or_ny, int(sqrt(nx*ny / TARGET_CANDIDATES))))
where TARGET_CANDIDATES=8, giving ~8 candidate quads per coarse box regardless of grid resolution.

Remove no longer used max_reg_list_num from quad_irreg_grid_coords

To do:

• regional
• wrap long regional
• pole regional
• global - bitwise

Fixes issue

Fixes #979

Types of changes

• Bug fix (non-breaking change which fixes an issue)
• New feature (non-breaking change which adds functionality)
• Breaking change (fix or feature that would cause existing functionality to not work as expected)
• Documentation update

Documentation changes needed?

• My change requires a change to the documentation.
  • I have updated the documentation accordingly.

Tests

Please describe any tests you ran to verify your changes.
Looking at brute force vs. current method vs new method: https://github.com/hkershaw-brown/quad_search_perf
I've only created global test grids. regional grids now,
Quad_util_developer tests bitwise, but this are fairly limited.

Currently running ROMS_rutgers case to compare memory (still in queue)

Checklist for merging

• Updated changelog entry
• Documentation updated
• Update conf.py

Checklist for release

• Merge into main
• Create release from the main branch with appropriate tag
• Delete feature-branch

Testing Datasets

• Dataset needed for testing available upon request
• Dataset download instructions included
• No dataset needed

[hkershaw-brown]

Q. global vs. regional.
if you initialize the interpolation as global, but the grid is regional how does this affect things?

wrap where there is no wrap. (mad wrapping of region around globe)
global, sparse -> slow search
global, not sparse, too many points in a box. -> error out or slow search

How bad is it to have mad-wrap? <-- -->

DART/models/utilities/quad_utils_mod.f90

Lines 764 to 769 in acf0a14

	! for a global grid, the initial values have already been set in
	! the derived type. otherwise, find the min/max of lons and lats.
	if (.not. h%opt%global_grid) then
	h%ii%min_lon = minval(h%ii%lons_2d)
	h%ii%max_lon = maxval(h%ii%lons_2d)
	h%ii%lon_width = h%ii%max_lon - h%ii%min_lon ! FIXME: wrap?

[hkershaw-brown]

or regional refinements in global grids:

[hkershaw-brown]

@mgharamti this is the quad search I am playing with. I think the main concerns are calculating the course grid from the actual grid rather than the hardcoded 3 cases, and the search structure size. But open to any input.

DART/models/utilities/quad_utils_mod.f90

Lines 730 to 751 in acf0a14

	!> Build the data structure for interpolation for an irregular quad grid
	subroutine init_irreg_interp(h)
	type(quad_interp_handle), intent(inout) :: h
	character(len=*), parameter :: routine = 'init_irreg_interp'
	! Need a temporary data structure to build this.
	! These arrays keep a list of the x and y indices of dipole quads
	! that potentially overlap the regular boxes.
	integer, allocatable :: reg_list_lon(:,:,:)
	integer, allocatable :: reg_list_lat(:,:,:)
	real(r8) :: u_c_lons(4), u_c_lats(4), pole_row_lon
	integer :: i, j, k, pindex, nx, ny, nrx, nry, istatus
	integer :: reg_lon_ind(2), reg_lat_ind(2), u_total, u_index
	logical :: cyclic, pole
	integer :: xlim
	allocate(reg_list_lon(h%ii%num_reg_x, h%ii%num_reg_y, h%ii%max_reg_list_num))
	allocate(reg_list_lat(h%ii%num_reg_x, h%ii%num_reg_y, h%ii%max_reg_list_num))

[mgharamti]

@hkershaw-brown, thanks for the work on this. I was always on the edge with the previous style of 3 hard-coded cases but didn't know how to go around that. I think I understand your design and strategy. CSR clearly avoids working with huge 3D arrays.

My first lazy question is why the 8 targets on average? I guess how did you choose it? It seems to me that search cost per obs is ~O(num_grids) and number of bins/boxes is ~O(N/TARGET_CANDIDATES) and so it's a performance tuning parameter. At least we need to explain it in the docs if you decide to go with it.

Second, could we add a post-fill consistency check that verifies fill_pos(i,j) == grid_start(i,j) + grid_num(i,j) for every coarse bin? This would catch any future divergence between the counting pass and filling pass, especially around edge cases such as longitude wrapping, masked/invalid quads, or boundary-touching cells.

[mgharamti]

Do we need to cycle here rather than just print an error? If a failure produces bad corners, both passes may count/fill bad bins or one failure path could behave differently depending on the state (right?)

[hkershaw-brown]

I will take a closer look - I think I am missing the logic in the original code.

The loop limits are the x,y bounds of the lon and lat arrays.
Is the print just checking if the loop limits are wrong?

If you can go off the end of the arrays,
it seems like both should fail, or both should pass.

If one fails and one passes I think this should be a catastrophic error rather than a print, since lon,lat are for the same point.

[hkershaw-brown]

moved to separate issue #1160

[hkershaw-brown]

My first lazy question is why the 8 targets on average? I guess how did you choose it? It seems to me that search cost per obs is ~O(num_grids) and number of bins/boxes is ~O(N/TARGET_CANDIDATES) and so it's a performance tuning parameter. At least we need to explain it in the docs if you decide to go with it.

My lazy choice: Dr. Google for what was a good choice for this kind of algorithm.
I was looking at brute force (best memory, but the worst case timing ~50,000 slower) vs. current hard coded method (max memory, target is to go as fast as that).
#979 (comment).

I think the equivalent average target (nx*ny/num_reg*num_reg)is ~10ish for the hard coded versions for the POP grids:

grid	nx	ny	num_reg	average quads per coarse box
gx1 (3deg?)	320	384	90	15.17037037
tx0.5 0.5deg	720	480	180	10.66666667
tx0.1 0.1 deg	3600	2400	900	10.66666667

The grids are not uniform so 'average' is probably a disgusting assumption.

I think it is well worth doing some runs with the test code and different (and variable) grid resolutions to count the actual number of quads per box with the new code (and counting the per course box with the 3 hardcoded POP cases). Maybe we have 799 quads in one box, 3 in another. Currently the 8 is just an ok choice where I was finding comparable speeds to the current version, but it is definitely a magic number that needs justifying.

[mjs2369]

Shouldn't we keep this check? Even though this is now just the first pass, I think this check would still work at the end of count_reg_overlaps

Assuming I'm understanding the initial motivation for this check correctly, I believe it is still possible for index_x or index_y to return invalid box numbers due to wraparound shenanigans

[mjs2369]

@hkershaw-brown last outstanding comment before approval

[hkershaw-brown]

Put check back in: b0c90c1

[mjs2369]

I think the updates here look really solid. Two-pass CSR looks clean and is an awesome way to reduce the memory overhead of the quad search. The computation time of get_quad_corners is tiny enough that doing it twice doesn't really make a difference in the runtime. Makes you wonder what possibilities could be out there for similar-ish updates elsewhere in the code

Just checking though that this was tested on sufficiently large grids - I see that the timing / memory comparison on quad_search_perf was done with size Nx = 5000, Ny = 4500
Did that ROMS_Rutgers test you mentioned ever finish? How big can we expect grid size to reasonably get, and would having a huge / fine grid change how well the two-pass CSR build performs time wise?

Dynamic calculation of the dimensions for the coarse grid is also clearly a huge win, but I am also curious about the value for TARGET_CANDIDATES being 8.

Earlier you mentioned you were planning to do more tests: "I think it is well worth doing some runs with the test code and different (and variable) grid resolutions to count the actual number of quads per box with the new code (and counting the per course box with the 3 hardcoded POP cases). "
Any findings from these?

[hkershaw-brown]

I think the updates here look really solid. Two-pass CSR looks clean and is an awesome way to reduce the memory overhead of the quad search. The computation time of get_quad_corners is tiny enough that doing it twice doesn't really make a difference in the runtime. Makes you wonder what possibilities could be out there for similar-ish updates elsewhere in the code

get_close uses a two pass, I am not sure why quad utils did not originally use a two pass. But for sure there are all sorts of places in the code that could be improved.
The non-two pass is in POP and CICE.

Just checking though that this was tested on sufficiently large grids - I see that the timing / memory comparison on quad_search_perf was done with size Nx = 5000, Ny = 4500 Did that ROMS_Rutgers test you mentioned ever finish? How big can we expect grid size to reasonably get, and would having a huge / fine grid change how well the two-pass CSR build performs time wise?

Dynamic calculation of the dimensions for the coarse grid is also clearly a huge win, but I am also curious about the value for TARGET_CANDIDATES being 8.

Its just a heuristic, here is the comment to Moha on the same question. It was ~10 on main for the hardcoded choice.

#1102 (comment)

Earlier you mentioned you were planning to do more tests: "I think it is well worth doing some runs with the test code and different (and variable) grid resolutions to count the actual number of quads per box with the new code (and counting the per course box with the 3 hardcoded POP cases). " Any findings from these?

This is one of the reasons for keeping the diagnostics prints in the log file about how many candidates are in each box.
Note there is also a bug(s) #1124 (in main and this pull request) with longitudes close to 0 and #978
This pull request is mostly to allow bigger grids that the hardcoded values in main, and use the same search algorithm as main but use less memory per core. I think the quad_utils_mod could be improved greatly https://github.com/NCAR/DART/issues?q=is%3Aissue%20state%3Aopen%20quad_utils however I have not made the attempt to do this in this pull request, just made issues as a found them.

[mjs2369]

Nice work Helen

[hkershaw-brown]

@mgharamti any remaining issues with this pull request? I believe I addressed the requested changes.

Cheers,
Helen

[mgharamti]

Looks good and ready to ship!