"Nearest neighbor" is different with different number of processors

[samsrabin]

Requirements

• Reviewed ESMF Reference Manual
• Searched GitHub Discussions

Affiliation(s)

NSF-NCAR

ESMF Version

No response

Issue

In CTSM, we use ESMF to read some input files. One particular pair of input files, specifying crop sowing window start and end dates, is at half-degree resolution. We tell ESMF to do nearest-neighbor¹ spatial interpolation as necessary to match the simulation grid.

When I do a run at 10°x15° resolution, some of the simulation gridcell centers are located exactly at the "corners" of four half-degree input pixels, meaning that those four neighbors are equally near. It doesn't matter to me which of those ESMF chooses as the "nearest neighbor," as long as it's consistent.

Unfortunately, it's not: At least one gridcell has a different "nearest neighbor" chosen depending on how many processors the job is split across.

As an example, I've made a figure based on two cases that are identical in setup except that Case 1 used 128 processors and Case 2 used 64. Due to this issue, a certain crop in the gridcell centered at latitude 0, longitude 30°E² gets sowing window of days 7-82 in Case 1 and 336-46 in Case 2.

The white/gray/black in this figure represents the half-degree sowing window files. Gray pixels match the values in Case 1, black pixels match Case 2, and white pixels match neither. The red lines intersect at the center of the 10x15 CTSM gridcell.

It looks like Case 1 reads from the pixel to the southwest, whereas Case 2 reads from the pixel to the northwest.

Some notes:

• I'm not 100% certain this is an ESMF issue as opposed to something weird that CTSM is doing, but I'm at the point where I've done all the troubleshooting I can within CTSM.
• This reproduces every time, over dozens of tests.

Tagging @ekluzek, @billsacks, and @briandobbins, who have expressed interest in this. By the way, I think I mentioned to y'all that I was having an ERP test pass but the equivalent PEM test fail—this is why! The read of sowing windows only happens at the very beginning of the test, so changing processor count halfway through makes no difference.

Autotag

@oehmke

Footnotes

1. It needs to be nearest-neighbor because dates are modulo—interpolating between Jan. 2 [day 2] and Dec. 31 (day 365) should give Jan. 1 (day 1), not July 3-4 (day [2+365]/2 = 183.5)—and that's not something ESMF can do, to my knowledge. ↩

2. There are other crops in this gridcell that also get different sowing windows. There are no crops in any other gridcell that get different sowing windows, but that doesn't necessarily mean different "nearest" neighbors are getting chosen. That might be happening, just with input pixels that don't differ. ↩

[samsrabin]

It looks like the specific nearest-neighbor method being used is ESMF_REGRIDMETHOD_NEAREST_STOD.

else if (trim(sdat%stream(ns)%mapalgo) == 'nn') then
   call ESMF_FieldReGridStore(sdat%pstrm(ns)%field_stream, lfield_dst, &
        routehandle=sdat%pstrm(ns)%routehandle, &
        regridmethod=ESMF_REGRIDMETHOD_NEAREST_STOD, &
        dstMaskValues=(/sdat%stream(ns)%dst_mask_val/), &
        srcMaskValues=(/sdat%stream(ns)%src_mask_val/), &
        srcTermProcessing=srcTermProcessing_Value, ignoreDegenerate=.true., &
        unmappedaction=ESMF_UNMAPPEDACTION_IGNORE, rc=rc)

I see that there's also a ESMF_REGRIDMETHOD_NEAREST_DTOS. Do we think that might work better?

[samsrabin]

Nope. That results in sowing window dates outside the acceptable range, even though the input files contain no such values.

[samsrabin]

From the ESMF_REGRIDMETHOD_NEAREST_STOD documentation, this might be the problem:

If two points are equally close, then the point with the smallest sequence index is arbitrarily used (i.e. the point which would have the smallest index in the weight matrix).

Could the size or shape of the weight matrix differ between processor counts? I feel like that would mean ESMF is using multiple processors for the regridding—is that true?

[samsrabin]

Maybe?

Regridding can be broken into two stages. The first stage is generation of an interpolation weight matrix that describes how points in the source grid contribute to points in the destination grid. The second stage is the multiplication of values on the source grid by the interpolation weight matrix to produce values on the destination grid. This is implemented as a parallel sparse matrix multiplication.

And earlier:

In the case of the regrid weight generation, user data is used to produce interpolation weights following specific numerical schemes. The bit-for-bit reproducibility of the generated weights depends on the implementation details. Section 24.2 provides more details about the bit-for-bit considerations with respect to the regrid weights generated by ESMF.

In the case of the sparse matrix multiplication, which is the typical method that is used to apply the regrid weights, user data is directly manipulated by ESMF. In order to help users with the implementation of their bit-for-bit requirements, while also considering the associated performance impact, the ESMF sparse matrix implementation provides three levels of bit-for-bit support. The strictest level ensures that the numerical results are bit-for-bit identical, even when executing across different numbers of PETs. In the relaxed level, bit-for-bit reproducibility is guaranteed when running across an unchanged number of PETs. The lowest level makes no guarantees about bit-for-bit reproducibility, however, it provides the greatest performance potential for those cases where numerical round-off differences are acceptable. An in-depth discussion of bit-for-bit reproducibility, and the performance aspects of route-based communication methods, such the sparse matrix multiplication, is given in section 36.2.1.

In Sect. 36.2.1:

In order to help users with the implementation of their bfb requirement, ESMF provides different levels of control over the term order in sparse matrix multiplications, while at the same time offering performance optimization options. In all there are three arguments that will be introduced in the following paragraphs: srcTermProcessing, termorderflag, and pipelineDepth.

Use ESMF_TERMORDER_SRCSEQ together with srcTermProcessing=0 or srcTermProcessing=1 when strict bfb reproducibility is required independent of the source Array distribution, e.g. for different number of PETs.

It looks like CESM does use ESMF_TERMORDER_SRCSEQ with srcTermProcessing=0.

pipelineDepth is not specified in the call of ESMF_FieldReGridStore(), which means it's "auto-tuned."

[samsrabin]

Setting pipelineDepth to 64 doesn't have any effect.

[billsacks]

Thanks for all of these details, @samsrabin . It does seem like the presence of equally-close points in the nearest neighbor, together with the statement that the choice is arbitrary in this case, is likely to be the issue. How time-critical is this? Specifically: can it wait a few weeks until @oehmke is back from vacation?

One thing you could try to verify that the problem is in ESMF is to try to do an identical regridding using the offline ESMF_RegridWeightGen tool. ESMF does its regridding in parallel, and you can try running that tool in a batch job using the same number of processors as are being used in the CTSM runs. Let us know if you'd like to try this, and if so I could try to help you set this up if you're not clear on how.

That said, this could be a rabbit hole: if you are able to reproduce it, that will clearly point to it being an issue within ESMF, but if you're not able to reproduce the problem in ESMF_RegridWeightGen, there could be questions of whether it's just hard to reproduce there for some reason or if the problem is truly elsewhere. So I understand if you don't want to go there right now, or if you'd like to turn it over to the ESMF team at this point.

[samsrabin]

Thanks Bill! Trying it outside of CESM with ESMF_RegridWeightGen sounds like a great idea. Do you know of a destination grid file I could use for this?

Unfortunately this is rather time-critical—if I can't resolve this by the end of the month, I can't bring in the new crop calendars as default.

[billsacks]

I'm not 100% sure, but I think the destination grid / mesh file should be the same as the mesh_lnd file specified in nuopc.runconfig. e.g., for a 10x15 run that I have sitting around, this is share/meshes/10x15_nomask_c110308_ESMFmesh.nc, though you should confirm that.

You could also print the arguments to the call that sets up the stream stuff to verify that this is the mesh that's being used here.

Thanks for the info on the time-critical nature of this. I'm concerned that, if this is an issue in ESMF, it seems unlikely that we'll have a fix in place and a new ESMF version ready in time for this timeline. It makes me wonder if there might be a workaround you could do to get this working well enough while waiting for a robust fix. One thing that comes to mind, for example, is adding some tiny value to all of the coordinates on the input data file to (generally) ensure that there is a unique nearest neighbor. I know that feels really kludgey, but I'm thinking it would resolve this issue for now without changing the science in any significant way. Thoughts?

[samsrabin]

Great thinking—that did it! I made a new version of the sowing window start and end dates with lat/lon values shifted by 1e-4 or 1e-5, and now the sowing windows are consistent across 128 or 64 processors.

I also had a go at ESMF_RegridWeightGen, but was confused about the --weight argument. Since I need to buckle down in order to get the crop calendars in, I think it's best if I leave that testing to the ESMF team to get to whenever.

Thanks again for your help!

[billsacks]

Great, I'm glad that resolved your problem for now. I'll leave this open so we can circle back to it and consider whether a more robust fix is feasible longer-term.

[billsacks]

Closing this discussion; we have opened #276 to fix the underlying issue.