particle to grid transfer in checked integer arithmetic
Particle methods for continuum physics, the material point method and its affine descendants, move mass and momentum from particles onto a grid, solve there, and move the result back. The scatter onto the grid is a sum, and when the sum is taken in floating point its result depends on the order the terms arrive in. On a GPU the terms arrive in whatever order the hardware schedules, so two runs of the same program on the same data give grids that differ in their last bits, and the difference grows from there. The usual remedies are sorting the particles so the order is fixed, or accumulating in fixed point. This page takes the second remedy to its end: the transfer step written in cint, where every quantity is an integer, the scatter is run in two different particle orders, and the two grids are compared node by node.
where the order gets in
Each particle carries a mass and a velocity, so a momentum, and sits somewhere inside a grid cell. The transfer gives each of the cell's nodes a share of the particle's mass and momentum, weighted by how close the particle is to that node, and a node collects shares from every particle near it. In floating point the share is a product that is rounded, and the node's total is a sum of rounded terms that is rounded again after each addition, so the total depends on the order of the additions. In parallel code the additions are atomic operations whose order is not specified, which is why the same run twice does not give the same bits. None of this is a bug in any one program. It is what the arithmetic is.
integers do not care about order
Integer addition is associative and commutative, so a checked sum of integers that completes is the same in every order, and the bounds line in the program, run before any particle is scattered, makes every order complete; if the shares themselves are integers the grid is the same in every order too. The shares are made integers by splitting rather than weighting. For a quantity q and a node fraction f, the far node's share is q times f divided by one, rounded once to the nearest integer, and the near node's share is q minus that. The two shares sum to q by construction, with no rounding to lose. In two dimensions the split is done along x and then each part along y, so a particle's mass lands on its four nodes in four integers that sum to the mass exactly, and the same for each component of momentum. Conservation is then not a property to be checked within a tolerance. It is an identity, and the program prints both sides of it.
the program
Positions are in units of 1/65,536 of a cell on an 8 by 8 grid, masses are 1 to 1,000, velocities are -1,000 to 1,000, and 256 particles are placed by a small generator so that every run sees the same cloud. The first line of output is the bounds analysis as an executable statement: the largest value any node can hold is the particle count times the largest mass times the largest speed, and if that product does not fit I64 the program stops there, before any particle is scattered. The scatter then runs twice, once in particle order and once in a scrambled order, and the second grid is compared with the first. The last lines gather velocity back for one particle from its four nodes, with the division rounded by name.
transfer.ci Run the first Run loads Python, about 12 MB, once
what the output shows
The particles' total mass is 126,536 and their total momentum is (9,985,365, -2,138,879), and the grid's totals are the same three numbers, not close to them. Between the two particle orders, 0 of 81 nodes differ, in mass or in either momentum component. The gather returns the mean velocity of one particle's cell, which is not the particle's own velocity and is not meant to be; it is the quantity the grid holds. Change the 97 and 13 that scramble the order to any odd multiplier and any offset, or reverse the loop, and the count stays at zero, because there is no order in which integers add up differently. The test block at the bottom checks the splitting identity for two thousand values under cint test.
the same scatter in floating point
transfer_float.c is the same scatter written the usual way, in C: each particle's mass and momentum weighted by its distance to each of the four nodes and the products added into the grid, with the same 256 particles from the same generator. It runs the scatter in particle order and in all 128 scrambled orders with an odd multiplier, and compares every grid with the first, bit for bit. Built with gcc 13.3.0 and clang 18.1.3 at -O2 -ffp-contract=off, both give these numbers:
| arithmetic | positions, in cells | nodes differing, order 97 and 13 | orders of 128 matching particle order | grid momentum, x |
|---|---|---|---|---|
| cint, integer (this program) | multiples of 1/65,536 | 0 of 81 | 128 | 9,985,365 |
| C, float | multiples of 1/65,536 | 72 of 81 | 0 | 9,985,364.8214111328 |
| C, double | multiples of 1/65,536 | 0 of 81 | 128 | 9,985,365 |
| C, double | multiples of 1/65,535 | 72 of 81 | 0 | 9,985,364.9999999963 |
In double the page's own particles show nothing, and that is worth saying plainly. Every position is a whole number of 65,536ths of a cell, so every weight is an exact binary fraction, every product fits in 53 bits and every sum the grid takes is exact, and exact addition does not care about order either. Divide the same positions by 65,535 instead and double differs at 72 of 81 nodes, and in float it differs either way. All of it runs on one thread of one processor, so the orders here are chosen rather than scheduled; on a GPU the hardware picks one per run.
(k * a + 13) % 256 with an odd multiplier a, every grid compared bit for bit with the particle-order grid. In float every order changes it: a = 1, which only rotates the particle list, changes 32 nodes, and every other order changes 65 to 76. The integer program changes 0 in every order. Float from transfer_float.c built with -DLIST; integer from the program above with 97 replaced by each a, run in cint_ref.data
| a | float, nodes differing | integer, nodes differing |
|---|---|---|
| 1 | 32 | 0 |
| 3 | 67 | 0 |
| 5 | 68 | 0 |
| 7 | 74 | 0 |
| 9 | 69 | 0 |
| 11 | 74 | 0 |
| 13 | 70 | 0 |
| 15 | 68 | 0 |
| 17 | 74 | 0 |
| 19 | 72 | 0 |
| 21 | 69 | 0 |
| 23 | 73 | 0 |
| 25 | 67 | 0 |
| 27 | 73 | 0 |
| 29 | 70 | 0 |
| 31 | 72 | 0 |
| 33 | 73 | 0 |
| 35 | 73 | 0 |
| 37 | 69 | 0 |
| 39 | 69 | 0 |
| 41 | 66 | 0 |
| 43 | 68 | 0 |
| 45 | 71 | 0 |
| 47 | 70 | 0 |
| 49 | 73 | 0 |
| 51 | 69 | 0 |
| 53 | 72 | 0 |
| 55 | 72 | 0 |
| 57 | 71 | 0 |
| 59 | 72 | 0 |
| 61 | 72 | 0 |
| 63 | 69 | 0 |
| 65 | 67 | 0 |
| 67 | 70 | 0 |
| 69 | 69 | 0 |
| 71 | 71 | 0 |
| 73 | 68 | 0 |
| 75 | 71 | 0 |
| 77 | 73 | 0 |
| 79 | 73 | 0 |
| 81 | 74 | 0 |
| 83 | 72 | 0 |
| 85 | 74 | 0 |
| 87 | 71 | 0 |
| 89 | 69 | 0 |
| 91 | 69 | 0 |
| 93 | 73 | 0 |
| 95 | 68 | 0 |
| 97 | 72 | 0 |
| 99 | 73 | 0 |
| 101 | 68 | 0 |
| 103 | 74 | 0 |
| 105 | 69 | 0 |
| 107 | 69 | 0 |
| 109 | 75 | 0 |
| 111 | 74 | 0 |
| 113 | 71 | 0 |
| 115 | 72 | 0 |
| 117 | 72 | 0 |
| 119 | 69 | 0 |
| 121 | 66 | 0 |
| 123 | 68 | 0 |
| 125 | 69 | 0 |
| 127 | 73 | 0 |
| 129 | 67 | 0 |
| 131 | 71 | 0 |
| 133 | 69 | 0 |
| 135 | 73 | 0 |
| 137 | 71 | 0 |
| 139 | 71 | 0 |
| 141 | 68 | 0 |
| 143 | 74 | 0 |
| 145 | 72 | 0 |
| 147 | 72 | 0 |
| 149 | 71 | 0 |
| 151 | 71 | 0 |
| 153 | 69 | 0 |
| 155 | 71 | 0 |
| 157 | 65 | 0 |
| 159 | 70 | 0 |
| 161 | 71 | 0 |
| 163 | 70 | 0 |
| 165 | 74 | 0 |
| 167 | 72 | 0 |
| 169 | 75 | 0 |
| 171 | 69 | 0 |
| 173 | 71 | 0 |
| 175 | 65 | 0 |
| 177 | 67 | 0 |
| 179 | 72 | 0 |
| 181 | 73 | 0 |
| 183 | 70 | 0 |
| 185 | 70 | 0 |
| 187 | 73 | 0 |
| 189 | 73 | 0 |
| 191 | 73 | 0 |
| 193 | 73 | 0 |
| 195 | 68 | 0 |
| 197 | 74 | 0 |
| 199 | 70 | 0 |
| 201 | 67 | 0 |
| 203 | 73 | 0 |
| 205 | 72 | 0 |
| 207 | 72 | 0 |
| 209 | 68 | 0 |
| 211 | 72 | 0 |
| 213 | 67 | 0 |
| 215 | 70 | 0 |
| 217 | 75 | 0 |
| 219 | 75 | 0 |
| 221 | 73 | 0 |
| 223 | 75 | 0 |
| 225 | 73 | 0 |
| 227 | 70 | 0 |
| 229 | 71 | 0 |
| 231 | 75 | 0 |
| 233 | 70 | 0 |
| 235 | 69 | 0 |
| 237 | 71 | 0 |
| 239 | 76 | 0 |
| 241 | 72 | 0 |
| 243 | 70 | 0 |
| 245 | 71 | 0 |
| 247 | 72 | 0 |
| 249 | 73 | 0 |
| 251 | 69 | 0 |
| 253 | 74 | 0 |
| 255 | 75 | 0 |
same bytes from the reference and the compiler
The compiled program writes the same five lines as the reference, byte for byte, and the test passes under both. The one emitted C file was built with four compilers on three machines, each at -O2 (/O2 for MSVC), and run once on each; every build exited with status 0, as a program that completes does:
| run | machine | stdout SHA-256 |
|---|---|---|
| cint_ref | this page | 1245bd0e97fc6a41a847d1ee90d0df3267f7bfef7bff05b372dff1b3831fbf07 |
| clang 18.1.3 | x86-64 Linux | 1245bd0e97fc6a41a847d1ee90d0df3267f7bfef7bff05b372dff1b3831fbf07 |
| gcc 13.3.0 | x86-64 Linux | 1245bd0e97fc6a41a847d1ee90d0df3267f7bfef7bff05b372dff1b3831fbf07 |
| MSVC 19.44 | x86-64 Windows 11 | 1245bd0e97fc6a41a847d1ee90d0df3267f7bfef7bff05b372dff1b3831fbf07 |
| Apple Clang 21.0.0 | arm64 macOS 27 | 1245bd0e97fc6a41a847d1ee90d0df3267f7bfef7bff05b372dff1b3831fbf07 |
The box below is the C that cintc emitted for the program as shipped, 1,868 lines for 118 lines of cint, most of it the checked arithmetic written out.
cintc (the cint compiler)what is not claimed
This is the transfer step and nothing else. There are no forces, no constitutive model, no time integration and no grid solve, so it is not a fluid or a granular simulation; it is the two places in one where the order of operations enters, isolated. The weights are linear hat functions, not the quadratic B-splines most material point codes use, and the affine velocity terms of the affine particle-in-cell method are not carried; both would be more splits of the same kind. The program runs on one thread and is not written as a kernel. The compiler now compiles kernels and runs them on the CPU and, through CUDA, on an NVIDIA GPU, but that GPU path is not yet a supported backend, so the parallel case that motivates the paper is not executed here. What is executed is the property the parallel case needs: that the grid does not depend on the order the particles arrive in, shown by arriving in two orders.
limits
The bound in the first line is for the demo's ranges; a real simulation declares its own ranges and the same line checks them. Momentum here is mass times an integer velocity with no further scale; a physical unit system would put a scale on velocity and the products would be wider, and I128 accumulators are available in the reference for that but not yet in the compiler. That floating-point accumulation with atomics on a GPU varies from run to run is the published finding of the GPU vendors cited below, not a measurement made here; what is measured here is the order dependence behind it, on one CPU thread with the orders chosen by hand.
references
- D. Sulsky, Z. Chen and H. L. Schreyer. A particle method for history-dependent materials. Computer Methods in Applied Mechanics and Engineering, 1994.
- C. Jiang, C. Schroeder, A. Selle, J. Teran and A. Stomakhin. The affine particle-in-cell method. ACM Transactions on Graphics, 2015.
- Y. Hu, Y. Fang, Z. Ge, Z. Qu, Y. Zhu, A. Pradhana and C. Jiang. A moving least squares material point method with displacement discontinuity and two-way rigid body coupling. ACM Transactions on Graphics, 2018.
- N. Whitehead and A. Fit-Florea. Precision and Performance: Floating Point and IEEE 754 Compliance for NVIDIA GPUs. NVIDIA, 2011.
- NVIDIA. cuBLAS Library documentation, results reproducibility: routines that use atomics are not guaranteed to give bit-wise reproducible results.
cite as
Harriett Little. Particle to grid transfer in checked integer arithmetic. integerc.dev, 2026. https://integerc.dev/papers/mpm/
Harriett Little
October 2026