La rejilla uniforme y el hash espacial

En cuanto las partículas dejan de ser independientes y empiezan a interactuar —chocar, repelerse, comportarse como un fluido— el coste salta de lineal a cuadrático. Un millón de partículas comparándose con todas las demás son un billón de pares por paso de tiempo, y ninguna GPU del planeta hace eso a sesenta cuadros por segundo. La estructura que rompe la cuadraticidad es más simple que cualquier árbol: dividir el espacio en cajas del tamaño del radio de interacción y aceptar que un vecino solo puede estar en la caja propia o en las adyacentes.

La cuenta que obliga a hacer algo

Con n partículas y una interacción entre pares, el número de pares es n(n-1)/2. Para diez mil partículas son cincuenta millones, que una GPU hace sin despeinarse. Para cien mil son cinco mil millones, que ya cuesta. Para un millón son quinientos mil millones, y con una evaluación de unos veinte ciclos por par estamos hablando de horas por paso.

La observación que lo arregla es que la interacción tiene un radio de corte h: más allá de esa distancia, la fuerza es cero. Una repulsión de contacto tiene el radio del diámetro de la partícula; un kernel de SPH tiene un soporte compacto por definición. Así que la inmensa mayoría de esos quinientos mil millones de pares se evalúan para dar cero.

Si la densidad de partículas es aproximadamente uniforme y hay d partículas por unidad de volumen, el número de vecinos dentro de una esfera de radio h es d·(4/3)π h³, una constante que no depende de n. El coste total pasa a ser n por esa constante: lineal. Todo el problema se reduce a encontrar esos vecinos sin recorrer el resto.

El tamaño de celda

La rejilla divide el espacio en cubos de lado L. Un vecino a distancia menor que h está, necesariamente, en la celda propia o en una de las adyacentes, siempre que L sea mayor o igual que h.

Con L = h, hay que mirar 3 por 3 por 3 celdas, o sea 27, y el volumen inspeccionado es 27h³. Con L = 2h, bastan 2 por 2 por 2 celdas —las que tocan el punto—, o sea 8, pero cada una tiene volumen 8h³, así que el volumen inspeccionado es 64h³.

La conclusión sorprende a mucha gente: la celda del tamaño del radio inspecciona menos de la mitad de volumen que la del doble, aunque haya que visitar más celdas. Y como el coste real es proporcional al número de partículas inspeccionadas, que es proporcional al volumen, L = h gana. El precio es más celdas y por tanto más memoria de índice.

De esas 540 inspeccionadas, solo las que caen dentro de la esfera de radio h interactúan de verdad: (4/3)π h³ · d, unas 84. O sea que se inspeccionan seis veces más de las necesarias. Ese factor 27/(4π/3), aproximadamente 6,4, es el precio inevitable de usar cajas para aproximar una esfera, y no hay estructura razonable que lo baje mucho.

Índice lineal para un dominio acotado

Si la simulación ocurre dentro de una caja conocida —un acuario, un recinto, el volumen visible— lo mejor es un índice directo, sin hash. La coordenada de celda sale de una división entera y el índice lineal de una fórmula de tres dimensiones.

struct Rejilla {
  origen:  vec3f,      // esquina minima del dominio
  celda:   f32,        // lado de la celda, igual al radio de interaccion
  dims:    vec3u,      // numero de celdas en cada eje
  total:   u32,        // dims.x * dims.y * dims.z
};
@group(0) @binding(0) var<uniform> g: Rejilla;

fn coordCelda(pos: vec3f) -> vec3i {
  return vec3i(floor((pos - g.origen) / g.celda));
}

fn indiceCelda(c: vec3i) -> u32 {
  // Saturar en vez de descartar: una particula que se escapa acaba
  // en la celda del borde en lugar de corromper el indice.
  let d = vec3i(g.dims);
  let s = clamp(c, vec3i(0), d - vec3i(1));
  return u32(s.x) + u32(s.y) * g.dims.x + u32(s.z) * g.dims.x * g.dims.y;
}

El clamp no es opcional. Sin él, una partícula que sale del dominio produce una coordenada negativa, la conversión a u32 da un número enorme, y el atomicAdd del conteo por celda escribe fuera de rango. WebGPU descarta la escritura por robustez, así que no revienta, pero esa partícula desaparece de la estructura de vecinos y su interacción se pierde de forma silenciosa. Saturar es peor físicamente y muchísimo mejor a la hora de depurar, porque el error se ve.

El coste de memoria es el que decide hasta dónde se puede llegar:

Con 4 bytes por celda, la de 256 al cubo ya son 64 MiB solo para el contador, y hacen falta al menos dos arrays de ese tamaño. La de 512 al cubo se sale del límite por defecto de maxStorageBufferBindingSize. En la práctica, 128 al cubo o 256 al cubo son los tamaños razonables, y eso fija la relación entre el radio de interacción y el tamaño del dominio.

Hash espacial para un dominio no acotado

Cuando el dominio no tiene límites —un sistema que se expande, un mundo abierto, una simulación con escalas muy dispares— no se puede reservar una celda por posición posible. La salida es una tabla hash: se calcula un número a partir de la coordenada de celda y se toma el módulo del tamaño de la tabla.

La función clásica es la de Teschner y colaboradores, que combina tres primos grandes con o exclusivo:

const P1: u32 = 73856093u;
const P2: u32 = 19349663u;
const P3: u32 = 83492791u;

fn hashCelda(c: vec3i, tabla: u32) -> u32 {
  let h = (bitcast<u32>(c.x) * P1)
        ^ (bitcast<u32>(c.y) * P2)
        ^ (bitcast<u32>(c.z) * P3);
  return h % tabla;
}

El bitcast en vez de una conversión de tipo es deliberado: convierte el patrón de bits del entero con signo, incluidas las coordenadas negativas, sin comportamiento indefinido. Y el módulo funciona mejor si el tamaño de tabla es primo; con una potencia de dos y una máscara es más rápido pero se pierden bits altos del hash y las colisiones se estructuran.

El precio del hash son las colisiones: dos celdas alejadas pueden caer en la misma entrada, y entonces la búsqueda de vecinos encuentra partículas que no lo son. Eso obliga a comprobar la distancia real de todas formas, cosa que ya se hacía, así que las colisiones no producen errores, solo trabajo desperdiciado.

La tasa de colisión sigue el problema del cumpleaños. Con m celdas ocupadas y una tabla de tamaño T, la fracción de entradas con más de una celda es aproximadamente 1 - (1 + m/T)·exp(-m/T). Un factor de carga m/T de 0,5 da alrededor de un 9 por ciento de entradas con colisión, y de 2 da alrededor de un 59 por ciento. La regla práctica es dimensionar la tabla al doble del número esperado de celdas ocupadas.

Con el dominio acotado no hay colisiones y el índice es más barato de calcular. Con el hash, la memoria no depende del tamaño del mundo. Elegir es fácil: si conoces la caja, índice directo; si no la conoces, hash.