Hacktoberfest 2026: los issues que los mantenedores marcaron para octubre, abiertos y aptos para principiantes. Explorar issues de Hacktoberfest

`IBMStepper`: the IBM loop is not a fixed point iteration — `ibm_tolerance` cannot fire, and `ibm_relaxation`/`ibm_max_iterations` act as a forcing gain

Abierto
  • ¿Algo no está bien?
#166 0 comentarios 0 reacciones 0 asignados Ver en GitHub

Nadie ha tomado este issue todavía.

Evaluación

Dificultad
4/5
Tiempo estimado
3-5 días
Aptitud para principiantes
50/100
Tipo de issue
Error
Claridad
Bastante claro
Estado de actividad
Tranquilo
Stack tecnológico
python
Área
backend, hpc

Línea de trabajo

Comienza en xlb/operator/stepper/ibm_stepper.py, en warp_implementation, y sigue compute_velocity_and_correct, interpolate_velocity_and_update_force y correct_population_ibm para confirmar cómo cambian f_1 y lag_forces entre sweeps. Reproduce el residual y el comportamiento de las iteraciones reportados; después determina si el resultado previsto es una corrección iterativa o un forzado explícito documentado; se considera terminado cuando el comportamiento elegido y su semántica de tolerancia, relajación y número máximo de iteraciones estén cubiertos por la validación.

Escrito por el modelo de indexación a partir del texto del issue.

Descripción

Summary

In xlb/operator/stepper/ibm_stepper.py, the loop in warp_implementation never updates f_1. The only kernel that writes f_1 — correct_population_ibm — is launched once after the loop exits. As a result every sweep recomputes the same u from the same populations, so the loop does not converge in the sense the docstring describes.

This makes the documented behaviour

  1. Residual-based stopping instead of a fixed iteration count. The loop monitors the maximum incremental change in Lagrangian forces and stops early when it falls below ibm_tolerance [...]

unreachable in practice.

Mechanism

Inside the loop:

  • compute_velocity_and_correct reads f_1 to get u and writes eul_velocities, but does not modify f_1.
  • interpolate_velocity_and_update_force interpolates eul_velocities back to the markers and does lag_forces += v_solid - u_interp.

Since f_1 is unchanged, u is identical on every sweep, hence delta_F = v_solid - u_interp is identical on every sweep. Therefore:

  • lag_forces after sweep k is exactly k * delta_F (a linear ramp)
  • the residual |lag_forces - lag_forces_prev| equals |delta_F| — a constant

A constant residual can never cross a threshold it did not already satisfy on the first sweep.

Reproduction

A faithful 1-D transcription of the kernel sequence (Peskin weights, Voronoi areas, same spread/interpolate operators), omega=1.0, 6 sweeps:

residual : 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
|lag|max : 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000

The residual is bit-for-bit constant; the accumulated force grows linearly.

We also instrumented the real Warp path on GPU and counted iterations: 100% of steps hit the ibm_max_iterations ceiling with the convergence flag still set, across two different resolutions and two different ibm_relaxation values (60 steps each). That matches the analysis.

Consequences

1. ibm_tolerance is effectively dead. The early exit never triggers; the loop always runs ibm_max_iterations times.

2. ibm_relaxation and ibm_max_iterations are a boundary forcing gain, not convergence controls. Working through the algebra, the correction applied after the loop is

eul_forces = omega * ( (n-1) * spread(delta_F)/weights  -  u )

so the gain on the target term is ibm_relaxation * (ibm_max_iterations - 1). Note the - 1: the final sweep spreads the forces accumulated by the previous sweep.

3. ibm_max_iterations = 1 silently drops the wall velocity entirely. With n = 1 the target gain is zero and only -omega * u is applied. For a static body that still drives u -> 0, which is the correct no-slip target, so it looks fine. For a moving body the prescribed surface velocity v_solid never enters the equations at all. This one is easy to miss because the static case masks it.

4. Stability. Because the scheme is an explicit forcing with gain rather than a contraction, the near-surface amplification factor flips sign and exceeds 1 as the gain rises. In our reference implementation, for a static body: omega=1.0, n=4 gives -2.61, omega=0.5 gives -0.805, omega=0.25 gives +0.098. We see the corresponding ordering in practice — at fine grid resolutions the run diverges at the higher gains and survives longer as the gain is lowered.

(The absolute numbers overstate the real amplification, since the reference implementation omits the collide/stream between corrections. The robust findings are the ordering and the sign flip.)

Suggested fix

Apply each sweep's correction to f_1 inside the loop, so the next sweep reads an updated u. Two details matter:

  • spread the per-sweep increment, not the accumulated lag_forces — the accumulated total is already baked into f_1 by the previous sweeps' corrections, so spreading it again double-counts;
  • drop the - u term in the normalization for the same reason.

With that change the residual actually decreases, ibm_tolerance and the early exit become meaningful, and the surface correction is monotonically damping for 0 < ibm_relaxation <= 1 (no sign flip at any of the gains above).

We have implemented this in our fork behind an opt-in flag, keeping the current behaviour as the default so existing results stay reproducible. Happy to open a PR against main if that shape is useful to you — either the flag approach or a straight fix.

Environment

  • xlb main (the code path above is unchanged as of today)
  • Warp backend, D3Q27, warp-lang 1.12.1

Alternative reading

If the current behaviour is intentional — i.e. this is meant to be a single explicit forcing step with a tunable gain rather than a multi-direct forcing iteration — then the fix is documentation rather than code: the ibm_tolerance parameter should be removed or marked non-functional, and the docstring's points (2) and (3) reworded, since ibm_relaxation and ibm_max_iterations are gain knobs. Either way the current docs describe behaviour the code does not have.

Lenguaje dominante
Python
Estrellas
508
Forks
85
Merge medio
4 d 18 h
PR fusionados (30 d)
1

Preparar el entorno

Primeros pasos

  1. Lee el issue completo y luego la guía de contribución del proyecto.
  2. Comenta en el issue que vas a ocuparte — evita que dos personas hagan lo mismo.
  3. Haz un fork del repositorio y trabaja en una rama.
  4. Abre un pull request que haga referencia al número del issue.

Más de Autodesk/XLB

Todos los issues de Autodesk/XLB

Issues similares

Más issues de Python

Recibe los nuevos issues en tu correo

Un resumen breve de issues de GitHub para principiantes.