Skip to content

Bug in grid point index calculation #67

@EikHen

Description

@EikHen

Hey,

in several places in the code, indices are computed like

short aStartingGridPoint = gridDimensionA_ * aCoord;

with some precautions to ensure 0 <= aCoord < 1 so that aStartingGridPoint < gridDimensionA_ should always evaluate to true.
If aStartingGridPoint >= gridDimensionA_, the code will run into an out-of-bounds access of an array.

Unfortunately, despite the precautions taken, this happens on extremely rare occasions. I do not completely understand why.
My best guess is some kind of numerical error originating in the way short * Real is rounded on some platforms. However, I couldn't reproduce this problem in a minimal example. My code (hundreds of thousands of charges and coordinates) consistently runs into this issue though.
I could confirm this is happening by inserting if (aStartingGridPoint == gridDimensionA_)-checks into the code.
Something like this fixes the problem, though a little dirtily:

inline short computeClampedIndex(unsigned int dim, Real coord) {
    short startingGridPoint = dim * coord;
    if (startingGridPoint == dim) {
        startingGridPoint -= 1;
    }
    return startingGridPoint;
}
...
short aStartingGridPoint = computeClampedIndex(gridDimensionA_, aCoord);

Any ideas why this is happening, despite the precautions which are already taken in the code to prevent this very problem?

Best,
Eik

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions