Fuerzas y campos de velocidad

El integrador es la mitad del sistema; la otra mitad es la función que devuelve la aceleración. Ahí es donde el sistema deja de ser una demostración técnica y empieza a parecer algo: humo, chispas, un enjambre, polvo arrastrado por el viento. Y hay dos herramientas que producen la mayor parte de los resultados convincentes con muy poco código, siempre que se entiendan sus trampas: los atractores puntuales, que se comportan mal cerca del centro si no se suavizan, y el ruido de rotacional, que produce turbulencia que se ve como un fluido porque cumple una propiedad que ningún ruido corriente cumple.

Acumular fuerzas, no aceleraciones

La distinción entre fuerza y aceleración se ignora en cuanto todas las partículas tienen la misma masa, y vuelve a importar en cuanto no la tienen. La gravedad es una aceleración —no depende de la masa— y un empuje o una repulsión son fuerzas. Mezclarlas mal produce el bug clásico de que las partículas pesadas caen más despacio.

fn aceleracionTotal(q: Particula, t: f32) -> vec3f {
  var fuerza = vec3f(0.0);

  fuerza += vec3f(0.0, -9.81, 0.0) * q.masa;      // gravedad: fuerza = m*g
  fuerza += atractor(q.pos, vec3f(0.0, 4.0, 0.0), 40.0);
  fuerza += vec3f(sin(t * 0.7), 0.0, cos(t * 0.4)) * 1.5;   // viento

  var a = fuerza / q.masa;
  a += curl(q.pos * 0.25 + vec3f(0.0, t * 0.2, 0.0)) * 6.0; // campo de velocidad
  a -= q.vel * 0.15;                                        // arrastre lineal
  return a;
}

Las dos últimas líneas no pasan por la masa a propósito: un campo de velocidad y un arrastre proporcional a la velocidad son efectos cinemáticos, no fuerzas, y aplicarlos igual a todas las partículas es lo que produce el aspecto de “todo se mueve con el fluido”.

Atractores y la singularidad

Un atractor newtoniano tiene aceleración proporcional a la inversa del cuadrado de la distancia. Cerca del centro esa expresión diverge: a distancia cero la aceleración es infinita, y a distancia pequeña es enorme, lo cual con cualquier paso de tiempo finito lanza la partícula a la otra punta del universo.

El remedio estándar es el suavizado de Plummer: se añade una constante al cuadrado de la distancia dentro de la raíz. La fuerza deja de diverger, se hace cero en el centro, y el sistema pasa a ser estable con pasos de tiempo razonables.

const EPS2: f32 = 0.25;      // radio de suavizado al cuadrado

fn atractor(pos: vec3f, centro: vec3f, fuerza: f32) -> vec3f {
  let d = centro - pos;
  let r2 = dot(d, d) + EPS2;             // nunca cero
  let inv = inverseSqrt(r2);             // 1/r
  // fuerza * d / r^3, expresado con inv para evitar divisiones
  return d * (fuerza * inv * inv * inv);
}

inverseSqrt es una instrucción de hardware en todas las GPU y sale más barata que dividir por sqrt. Elevarla al cubo evita las dos divisiones que tendría la forma literal.

El valor de EPS2 es una decisión de diseño, no un detalle: fija el radio a partir del cual el atractor deja de tirar más fuerte. Si es demasiado pequeño el sistema es rígido y necesita pasos cortos; si es demasiado grande el atractor se nota blando. Un valor del orden del tamaño visible del atractor funciona bien.

Un repelente es lo mismo con la fuerza negativa. Y varios atractores se suman sin más, con un bucle sobre un array pequeño en un uniform:

struct Atractor { centro: vec3f, fuerza: f32 };
@group(0) @binding(2) var<uniform> attrs: array<Atractor, 8>;

fn atractores(pos: vec3f, n: u32) -> vec3f {
  var f = vec3f(0.0);
  for (var k: u32 = 0u; k < n; k = k + 1u) {
    f += atractor(pos, attrs[k].centro, attrs[k].fuerza);
  }
  return f;
}

Ocho atractores por partícula y un millón de partículas son ocho millones de evaluaciones, que a una GPU no le cuestan nada porque el array de atractores cabe entero en caché y todas las invocaciones leen lo mismo. Es el patrón contrario al de la interacción entre partículas, que sería un millón al cuadrado.

Arrastre lineal y cuadrático

El arrastre lineal, -c·v, corresponde a un fluido viscoso a baja velocidad. Es el que estabiliza el integrador y el que hace que las partículas alcancen una velocidad terminal suave. Su gran virtud es que es incondicionalmente estable con Euler semi-implícito mientras h·c sea menor que 2.

El arrastre cuadrático, -c·|v|·v, corresponde a resistencia aerodinámica real y es el que produce la sensación de peso: una partícula lanzada rápido frena bruscamente al principio y luego se deja llevar.

fn arrastreLineal(v: vec3f, c: f32) -> vec3f { return -v * c; }

fn arrastreCuadratico(v: vec3f, c: f32) -> vec3f {
  let m = length(v);
  if (m < 1e-5) { return vec3f(0.0); }     // evitar dividir por cero
  return -v * (c * m);
}

El cuadrático tiene un límite de estabilidad que depende de la velocidad, así que con velocidades altas y pasos largos puede invertir el signo y hacer que la partícula acelere hacia atrás. La saturación de velocidad del integrador lo contiene, pero conviene ser consciente.

Una tercera forma de arrastre que se usa mucho y no es física es la amortiguación multiplicativa: v = v * pow(amortiguacion, h). Es incondicionalmente estable para cualquier paso, independiente del framerate gracias al exponente, y da un control directo sobre cuánta velocidad se conserva por segundo. Para partículas visuales suele ser lo más cómodo.

Ruido de rotacional

Un campo de velocidad turbulento hecho con tres ruidos independientes para las tres componentes tiene un defecto que se ve enseguida: tiene divergencia. Hay puntos de los que sale más flujo del que entra, así que las partículas se acumulan en unos sitios y se vacían de otros, y el resultado parece polvo estático en vez de fluido.

El rotacional de cualquier campo vectorial tiene divergencia cero, por identidad matemática. Así que si se construye un campo vectorial de ruido —el potencial— y se toma su rotacional, se obtiene un campo de velocidad que no comprime ni expande el volumen. Es la técnica de Bridson, Houriham y Nordenstam, y es la razón de que la turbulencia procedural de las demos se vea como humo.

fn hash31(c: vec3i) -> f32 {
  var h = (bitcast<u32>(c.x) * 73856093u)
        ^ (bitcast<u32>(c.y) * 19349663u)
        ^ (bitcast<u32>(c.z) * 83492791u);
  h = (h ^ (h >> 13u)) * 1274126177u;
  h = h ^ (h >> 16u);
  return f32(h) * (1.0 / 4294967296.0);
}

fn ruido3(p: vec3f) -> f32 {
  let i = vec3i(floor(p));
  let f = fract(p);
  let u = f * f * (3.0 - 2.0 * f);        // suavizado de Hermite
  let n000 = hash31(i + vec3i(0, 0, 0));
  let n100 = hash31(i + vec3i(1, 0, 0));
  let n010 = hash31(i + vec3i(0, 1, 0));
  let n110 = hash31(i + vec3i(1, 1, 0));
  let n001 = hash31(i + vec3i(0, 0, 1));
  let n101 = hash31(i + vec3i(1, 0, 1));
  let n011 = hash31(i + vec3i(0, 1, 1));
  let n111 = hash31(i + vec3i(1, 1, 1));
  let x00 = mix(n000, n100, u.x);
  let x10 = mix(n010, n110, u.x);
  let x01 = mix(n001, n101, u.x);
  let x11 = mix(n011, n111, u.x);
  return mix(mix(x00, x10, u.y), mix(x01, x11, u.y), u.z);
}

// Potencial vectorial: tres ruidos decorrelacionados por desplazamiento.
fn potencial(p: vec3f) -> vec3f {
  return vec3f(ruido3(p),
               ruido3(p + vec3f(31.41, 17.13,  7.77)),
               ruido3(p + vec3f(62.83, 91.19, 45.45)));
}

// Rotacional por diferencias centradas.
fn curl(p: vec3f) -> vec3f {
  let e = 0.15;
  let dx = vec3f(e, 0.0, 0.0);
  let dy = vec3f(0.0, e, 0.0);
  let dz = vec3f(0.0, 0.0, e);

  let px1 = potencial(p + dx); let px0 = potencial(p - dx);
  let py1 = potencial(p + dy); let py0 = potencial(p - dy);
  let pz1 = potencial(p + dz); let pz0 = potencial(p - dz);

  let inv = 1.0 / (2.0 * e);
  return vec3f(
    ((py1.z - py0.z) - (pz1.y - pz0.y)) * inv,
    ((pz1.x - pz0.x) - (px1.z - px0.z)) * inv,
    ((px1.y - px0.y) - (py1.x - py0.x)) * inv,
  );
}

El coste es alto: seis evaluaciones del potencial, o sea dieciocho ruidos, y cada ruido son ocho hashes y siete interpolaciones. Del orden de mil instrucciones por partícula. Con un millón de partículas eso satura la unidad aritmética de cualquier GPU, y es uno de los pocos kernels de partículas que no está limitado por memoria.

Dos formas de abaratarlo. La primera es precalcular el potencial en una textura 3D y muestrearlo con el hardware, que hace la interpolación trilineal gratis; el rotacional sigue necesitando seis muestras pero cada una cuesta una instrucción. La segunda es evaluar el campo con menos frecuencia que el integrador: si el campo cambia lentamente, actualizarlo cada cuatro pasos de simulación y guardar el resultado en la partícula es visualmente indistinguible y cuesta la cuarta parte.

Animar el campo es sumar tiempo a una de las coordenadas del punto de muestreo, que es lo que hace la llamada del primer bloque de código: desplazar en y produce un campo que asciende, como una columna térmica.