Improve precision in particle position updates

[hsinhaoHHuang]

Problem Description

For simulations with particles, GAMER currently stores particle positions relative to the edge of the entire simulation box. This method leads to numerical precision problems when:

• Large Simulation Box: When the simulation box is very big, the stored positions are large numbers.
• Small Timesteps: In high refinement levels, the timestep is very small, so the position update (calculated as velocity × timestep) is tiny.

In these scenarios, particle position updates (calculated as velocity × timestep) are significantly affected by floating-point rounding errors. In the worst cases, particles fail to move at all when the calculated displacement falls below the least significant digit of the position value.

For example, consider a particle with

• Position: ParPosX ~ 1.0e2 (code units)
• Velocity: ParVelX ~ 1.0e1 (code units)
• Timestep: dt ~ 1.0e-7 (code units)

The position update magnitude relative to the coordinate value would be:
dx/ParPosX ~ ParVelX*dt/ParPosX ~ (1.0e1 × 1.0e-7) / 1.0e2 ~ 1.0e-8
This is too small to be accurately represented in single-precision floating-point numbers. As a result, the update might be ignored due to rounding errors, and the particle might not move at all.

Using FLOAT8_PAR (double precision) can help by increasing precision. But it may lead to higher memory usage and slower performance.

Possible Solution

We could store particle positions with coordinates relative to the edge of the patch where each particle belongs, rather than the global simulation box. By keeping positions local, small updates will be captured accurately.

When absolute coordinates are needed (e.g., when writing output data), we can reconstruct them by adding the patch coordinates using double-precision arithmetic.

This change may require modifications to routines for:

• Initialization:
  • Change the routine for assigning initial particle data so that positions are stored relative to their patch.
• Particle Updates:
  • Modify the update routine to handle local coordinates and to manage particles moving between patches.
• Output Data:
  • Adjust the output routine to reconstruct the global coordinates from the patch-relative values when necessary.
  • Introduce a runtime parameter to control the floating-point precision of particles stored in snapshots, which may differ from the runtime precision.
• Collectting particle positions:
  • Such as Par_CollectParticle2OneLevel.cpp and Par_LB_CollectParticle2OneLevel.cpp.
• Whenever the global coordinates are needed.

[hyschive]

@hsinhaoHHuang Thanks for the summary.

Note that routines responsible for collecting particle positions, such as Par_CollectParticle2OneLevel.cpp and Par_LB_CollectParticle2OneLevel.cpp, also need to be revised.

We may introduce a runtime parameter to control the floating-point precision of particles stored in snapshots, which may differ from the runtime precision.

Please add these points to the feature request summary. Thanks!

Before implementing this feature, we should consider whether to adopt double precision for particles by default. I don't think it would significantly impact performance as particles are handled by the CPU. However, it would increase memory usage and snapshot sizes. @jzuhone Do you have any preference on this?

[jzuhone]

Most codes adopt double precision for at least particle positions in outputs by default--I'd be in favor of having it by default for all floating-point fields. In most cases the grid fields dominate the disk space anyway.

[hyschive]

Most codes adopt double precision for at least particle positions in outputs by default--I'd be in favor of having it by default for all floating-point fields.

I assume you mean all floating-point particle fields. I prefer to keep the grid fields in single precision by default, as switching them to double precision would significantly impact both performance and memory usage.

@hsinhaoHHuang Could you file a PR to enable --double_par by default?

[hsinhaoHHuang]

@hsinhaoHHuang Could you file a PR to enable --double_par by default?

No problem.

[jzuhone]

I assume you mean all floating-point particle fields. I prefer to keep the grid fields in single precision by default, as switching them to double precision would significantly impact both performance and memory usage.

Yes, that's what I meant.