[REQ] Efficient Persistent Neighbor List for Bonded Particle Simulations (DEM)

[jack123-gif]

Description

Suggestion: Implement native, efficient support within Warp for managing persistent neighbor relationships, specifically tailored for bonded particle contacts in Discrete Element Method (DEM) simulations.

Context

Motivation:

Warp's current hash grid implementation (wp.hash_grid and wp.hash_grid_query) is exceptionally fast and efficient for simulating granular materials (non-bonded particles) where contacts are transient and primarily based on spatial proximity each timestep.

However, many critical engineering and scientific applications involve cohesive or bonded materials. Examples include:

• Geomechanics: Rock masses, cohesive soils, cemented sands.
• Materials Science: Sintered powders, ceramics, composites, particle agglomeration.
• Civil Engineering: Asphalt, concrete degradation.
• Biomechanics: Tissue modeling, cell clusters.

The Problem with Current Methods:

In simulations with bonded materials, particle bonds (representing chemical cementation, sintering necks, physical cohesion, etc.) form and must persist based on physical criteria (e.g., stress/strain limits), even if the bonded particles temporarily move slightly outside the standard spatial query radius (max_dist in wp.hash_grid_query).

The current approach, relying solely on wp.hash_grid_query within the force calculation kernel each timestep, faces significant challenges for bonded materials:

1. Incorrect Bond Dropping: If bonded particles i and j separate slightly beyond max_dist due to vibrations or relative motion, wp.hash_grid_query will no longer report them as neighbors. Consequently, the bond force calculation is skipped, effectively breaking the bond prematurely and leading to physically incorrect simulation results.
2. Performance Penalty of Manual Management: Implementing bond management manually within user kernels (e.g., creating explicit neighbor lists per particle, using atomic counters to build global pair lists (i, j) where i < j, searching these lists to check for existing bonds before adding new ones, managing bond breakage flags, and performing list compaction) introduces severe performance overhead. This involves:
   • High atomic contention when updating global lists.
   • Complex data management logic within kernels.
   • Potentially poor memory access patterns.
   • Significant reduction in the performance advantage that Warp otherwise provides.

Proposed Solution:

We request the introduction of a dedicated, optimized mechanism within Warp for managing persistent neighbor relationships (bonds). This could potentially take the form of:

• A new Warp data structure (e.g., wp.BondList, wp.PersistentContacts).
• Associated functions/methods optimized for GPU execution, potentially including:
  • wp.bond_list_build(...): A function to update the bond list based on particle positions and user-defined criteria kernels (for both bond formation and breakage). This should handle duplicate prevention internally (e.g., only storing (i, j) if i < j).
  • wp.bond_list_query_pairs(...): A function that efficiently returns arrays (bond_indices_i, bond_indices_j) or an iterable object representing the currently active bond pairs, suitable for direct use in force calculation kernels.
  • wp.bond_list_compact(...): (Optional but highly useful) An efficient parallel compaction function to remove broken bonds and keep the list dense for optimal iteration performance.

Why Native Support is Crucial:

• Performance: A native C++/CUDA implementation managed by Warp can leverage highly optimized parallel algorithms (prefix sums, stream compaction, specialized hash tables or graph structures) far beyond what users can efficiently implement purely at the Warp kernel level. This is essential for maintaining high performance.
• Correctness & Robustness: Provides a reliable, tested mechanism for handling bond persistence, eliminating a major source of errors in user implementations.
• Ease of Use: Greatly simplifies the complexity for users wanting to simulate bonded materials, making Warp more accessible and powerful for a broader range of scientific and engineering problems.

Additional Considerations for Implementation:

• Bond Data Association: How can users associate custom data (e.g., initial bond length, stiffness, strength, damage variables, bond type ID) with each bond? Allowing linkage to user-defined wp.struct arrays indexed by the bond pair would be powerful.
• Criteria Flexibility: The bond formation and breakage checks should ideally accept user-defined wp.func or kernel pointers for maximum flexibility.
• Efficient Existence Check: A fast mechanism to check if a bond between i and j already exists would be beneficial, especially during the bond formation step.
• Memory Management: Clear strategies for allocation, potential resizing, or handling overflow if the number of bonds exceeds the initial capacity.
• Workflow: How would this system integrate with wp.hash_grid_query? (e.g., using wp.hash_grid_query for standard repulsive/frictional contacts and potential new bond detection, then using the bond system for persistent force calculation and breakage checks).

Benefits:

• Enables efficient and physically accurate simulation of bonded particulate systems, significantly expanding Warp's applicability.
• Reduces user code complexity and potential for errors.
• Aligns Warp's capabilities more closely with features found in established (but often slower) CPU-based DEM packages.

Thank you for considering this request. Adding robust support for bonded contacts would be a highly valuable addition to the Warp framework.

[mmacklin]

Hi @jack123-gif , many thanks for the message - one question regarding this comment:

Incorrect Bond Dropping: If bonded particles i and j separate slightly beyond max_dist due to vibrations or relative motion, wp.hash_grid_query will no longer report them as neighbors. Consequently, the bond force calculation is skipped, effectively breaking the bond prematurely and leading to physically incorrect simulation results.

Would one solution be to query with a conservative bounds (e.g.: max_dist + epsilon) so you have some margin in which to handle the break event?

I am concerned that things like BondList may be too specific for the core library, but perhaps we could build these on top of the existing data structures as part of a domain specific library or extension.

I will also ask some folks internally who work with these types of representations to comment.

Cheers,
Miles

[jack123-gif]

Hi @jack123-gif , many thanks for the message - one question regarding this comment:

Incorrect Bond Dropping: If bonded particles i and j separate slightly beyond max_dist due to vibrations or relative motion, wp.hash_grid_query will no longer report them as neighbors. Consequently, the bond force calculation is skipped, effectively breaking the bond prematurely and leading to physically incorrect simulation results.

Would one solution be to query with a conservative bounds (e.g.: max_dist + epsilon) so you have some margin in which to handle the break event?

I am concerned that things like BondList may be too specific for the core library, but perhaps we could build these on top of the existing data structures as part of a domain specific library or extension.

I will also ask some folks internally who work with these types of representations to comment.

Cheers, Miles

Thank you for the thoughtful reply and for taking the time to consider this. The solution you proposed (querying with conservative bounds like max_dist + epsilon) is definitely an interesting way to handle break events.

Regarding the concern about data structures like BondList being too specific for the core library, I'd like to offer a bit more context.

From the perspective of the Discrete Element Method (DEM), calculations can be broadly divided into two main categories:

Particle-based (Category 1): This is the method represented by the current Warp examples, which effectively "skips" the explicit concept of a contact and calculates forces directly between neighboring particles.

Contact-based (Category 2): This is a far more common method in the field. In this approach, all calculations (e.g., forces, torques, bond breakage) are first performed at the contact interface that forms between two particles. These results are then applied back to the two respective particles.

From a data structure standpoint, this second method does require building and maintaining a "neighbor contact list" (or "contact list"), which, as you noted, leads to higher memory consumption.

However, the significant advantage of this contact-centric approach is that it guarantees a high degree of computational accuracy. This level of precision is often essential for both scientific research and robust industrial applications, where simulating the exact physics of the contact is critical.

This is why I believe such a feature isn't too domain-specific for the core library; rather, it's a fundamental method for a large and important class of physics-based simulations that Warp is perfectly suited to accelerate.

Of course, this perspective is just based on my experience, having worked in the DEM field for about five or six years now.

Thank you again for the discussion and for escalating this to your colleagues. I appreciate it.

Best regards,