Add `ADR` integration by JBludau · Pull Request #284 · ORNL/CabanaPD

JBludau · 2026-01-28T21:30:26Z

This is based upon the formulation in https://www.sciencedirect.com/science/article/pii/S0020768317304833 but without the rotational stiffness.

One important difference is: in (75) we use (L_u(t-1)/Delta_ii-L_u(t)/Delta_ii). This corresponds to a -1*(75). Otherwise we get negative damping coefficients due to the approximated stiffness term being negative (the force difference points opposite to the velocity).

The implementation tries to be independent of CabanaPD's particle interface. This allows for easier testing (testcase is just a single mass) and potential reuse.

Add adapters for particle interface
Add testcase with particles
~~Add interface in Solver class~~ Not sure if that is a good idea ... we need to discuss this probably. In general the interface is still a bit ugly. The points are marked with TODO in comments
Add convenience functions for creation of an ADRIntegrator
Convert Dogbone case to use this method
MPI support/testing (this probably revolves around bounding the iteration to exclude the frozen particles and end at the offset)
Think about how to guess c for non PMB materials

Related to issue #207

JBludau · 2026-01-28T21:39:50Z

+                    stiffness[i] = ( -forces( index, i ) +
+                                     _forces_last_step( index, i ) ) /
+                                   mass[i] / _delta_t /
+                                   _velocities_last_step( index, i );


This is where it differs from the paper

JBludau · 2026-02-01T00:49:45Z

+#include <fstream>
+#include <iostream>
+
+#include "mpi.h"


hmm currently this does not support mpi ... maybe I should remove it

JBludau · 2026-02-01T00:51:20Z

+                         ParticleType const& particles )
+    {
+        auto forces = particles.sliceForce();
+        _integrator.initialSubStep( exec_space, forces );


Hmm...right now the integrator cycles over all particles ... that is probably false and should be corrected somehow ... I guess the best would be to create proxy objects for the slices

JBludau · 2026-02-02T15:32:24Z

+    integrator.middleSubStep( exec_space, particles, args... );
+
+    // Add non-force boundary condition.
+    if ( !boundary_condition.forceUpdate() )
+        boundary_condition.apply( exec_space, particles, time );
+
+    integrator.finalSubStep( exec_space, particles, args... );


This is the reason why we need a middle step. We calculate the velocity in the middle step but need to overwrite it with the bc if the bc sets a velocity. After that is applied we can go and integrate the displacement

yaowwwwwww · 2026-02-02T21:30:12Z

"Think about how to guess c for non PMB materials". This is a big question for the application. People cannot calibrate the c for each material, since the material behavior is different.

I am trying to define a bond-based JC strain rate-dependent model, and i define the equivalent plastic strain as the weighted root-mean-square norm of norm of bond plastic stretches in PMB model. The theoretical consistency needs further study.

Maybe use a reasonable mapping method of bond deformation to the material point strain, it will provide a easy way for defining different classical material in the future.

JBludau · 2026-02-02T22:01:22Z

"Think about how to guess c for non PMB materials". This is a big question for the application. People cannot calibrate the c for each material, since the material behavior is different.

I am trying to define a bond-based JC strain rate-dependent model, and i define the equivalent plastic strain as the weighted root-mean-square norm of norm of bond plastic stretches in PMB model. The theoretical consistency needs further study.

Maybe use a reasonable mapping method of bond deformation to the material point strain, it will provide a easy way for defining different classical material in the future.

I guess we would be fine with guessing a minimal value for c. Due to the ADR using c in the numerator for the fictitious masses, having a higher c than what the force uses increases the fictitious masses and makes the ADR converge slower (The implementation has a safety-factor of 5 anyways). I mainly expect numerical problems if c is far too low. In this case the ADR uses particles that are too light and thus they are underdamped, allowing them to oscillate. If that gets severe, it might lead to large errors in the time integration and thus to wrong solutions which might be hard to detect.

Maybe guessing c for the ADR as we do for PMB via 6E/(pihorizon^4(1-2nu)) is a good approximation for many cases. Basically as long as we guess the maximal stiffness we might be ok (and we still have the safety factor). Nevertheless, if we want to tune for performance that still might be tricky. But other factors might get us there faster (e.g. using something similar to a multi-scale simulation or even just using reduced precision for the ADR)

JBludau

The recent addition uses the specific c for every material instead of a user given c for multimaterial cases.

JBludau · 2026-04-29T18:41:14Z

-#include <CabanaPD_Integrate.hpp>
-#include <CabanaPD_Particles.hpp>
-#include <CabanaPD_config.hpp>
+#include <CabanaPD.hpp>


Not something I like, but the multimaterial and particle together need almost everything

JBludau · 2026-04-29T18:42:34Z

+          std::enable_if_t<current == ParameterPackType::size - 1, int> = 0>
+KOKKOS_INLINE_FUNCTION auto recursion_functor_for_index_in_pack_with_args(
+    FunctorType const& functor, std::size_t index,
+    ParameterPackType& parameterPack, ARGS... args )
 {
    if ( index == current )
    {
        return functor( Cabana::get<current>( parameterPack ), args... );
    }
    else
    {
-        return recursion_functor_for_index_in_pack_with_args<current + 1>(
-            functor, index, parameterPack, args... );
+        Kokkos::abort( "Requested index not contained in ParameterPack" );
+        return functor( Cabana::get<current>( parameterPack ), args... );


This is basically just a reversal of order of the two functions in the file. The old order was spoofing the compiler to not see the second SFINAE version.

JBludau · 2026-04-29T18:45:32Z

+        for ( unsigned i = 0; i < IndexingType::NumBaseModels; ++i )
+        {
+            _c[i] = CabanaPD::Impl::run_functor_for_index_in_pack_with_args(
+                GetCFunctor{}, i, _models.models );
+        }


This is the critical part where the multimaterial differs from single material. We need to get the correct stiffness for a specific material and store it. But we only get the basic set of materials, not the averaged ones. The reason for that is in the comment line 406

JBludau · 2026-04-29T18:49:01Z

+    double KOKKOS_FUNCTION operator()( IndexType index, int ) const
+    {
+        auto materialIndex =
+            _indexing( _particleType( index ), _particleType( index ) );
+        return _safety_factor *
+               ( _delta_t * _delta_t * Kokkos::numbers::pi * _horizon *
+                 _horizon * _horizon * _c[materialIndex] ) /
+               ( 3.0 * _delta_x );
+    }


During the integration, we iterate over the particles and only have one particle index. That means we can only have one c that we use for the fictitious mass even though the particle we are integrating might only use this c in some of its bonds.
But currently I don't see a way for us to get a better estimate for c with reasonable effort ... we could average c over the neighborhood and store it for every particle, but that would be a lot of new algorithmic lift and data. And it would need to be recomputed every step as bonds can break.

test

…n one step

reached

though we are using adr

JBludau self-assigned this Jan 28, 2026

JBludau requested review from pabloseleson and streeve January 28, 2026 21:31

JBludau commented Jan 28, 2026

View reviewed changes

JBludau force-pushed the ADL_integration branch 2 times, most recently from d5e5fef to 48a5aae Compare January 28, 2026 21:54

JBludau commented Feb 1, 2026

View reviewed changes

JBludau commented Feb 2, 2026

View reviewed changes

JBludau marked this pull request as ready for review February 2, 2026 15:33

JBludau force-pushed the ADL_integration branch from bcba56d to a11ecdf Compare April 29, 2026 18:37

JBludau commented Apr 29, 2026

View reviewed changes

JBludau added 16 commits June 15, 2026 16:14

added adl integrator

aa2b74b

add test with a single mass for adr integration

4170586

use exec_space in rangepolicies

1a85a57

remove stray comment

360b932

rework SOLID comment

bc055b5

renaming of first test into SingleMass in preparation for a particle

0dd2593

test

add a particle interface wrapper

a3eafe2

add a particle based ADR test

25c603f

add reset method to wrapper

84bbd90

add convenience function for creation

088ef2b

add nanchecks for the stiffness and damping terms

88462a1

add control functions to integrate particles outside of the solver

072241c

formatting

10797f2

add a dogbone example with ADR

f6c6505

split steps and add a middle step to get correct bc behavior

3bd55ec

adapt test case to use new integration with middle step

2b90589

JBludau added 21 commits June 15, 2026 16:14

use adr steps that are not just the full linear part of the problem i…

9c91ec2

…n one step

extend check to nan and inf

f77ca29

use full stiffness matrix instead of lumped one

1a0892a

tighten test to 1e-16

339961e

add multimaterial adr integration

da2e003

add multimaterial adr integration test

7b7ac8f

have indexing export the number of base models

6d0221b

reverse order of SFINAE set as compiler was struggling

f20b215

indexing typo

55d5a6e

remove commented includes

c681ed0

add a nofail switch class to allow setting all the domain to nofail

2a78f64

add a l2 norm evaluation to the adr integration

7e98106

allow gridSize to be called from a const

2c2b596

add a method that allows to set nofail and iterates until convergence is

1da93e1

reached

change unit tests to allow iterating unitl converged

94737f6

Use adr with single solver in dogbone_ADR

f7b041a

add comment to explain why we need the small timestep in this case even

694e9f2

though we are using adr

move force extracton to earlier point

1355bf6

fix solver and allow force and displacement tolerance

4378bdb

update tests

8c2892b

lower safety factor in dogbone example from 5 to 2

c4d7798

JBludau force-pushed the ADL_integration branch from 535dc64 to c4d7798 Compare June 15, 2026 20:14

JBludau added 8 commits June 15, 2026 16:17

use exec_space for fencing in NoFailSwitch

deb5c8c

add maybe_unused attributes

4b8b2e8

rename adr mass structs to express more specificity

427897d

add exact particle mass by integrating over the neighborhood

b4f1d0c

actually use safety factor

fbc2a10

use new naming of creator functions

bf66e88

use CLASS_LAMBDA for exact mass

6526639

move to extra function to make cuda happy

b53fcfb

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add `ADR` integration#284

Add `ADR` integration#284
JBludau wants to merge 45 commits into
ORNL:mainfrom
JBludau:ADL_integration

JBludau commented Jan 28, 2026 •

edited

Loading

Uh oh!

Uh oh!

JBludau Jan 28, 2026

Uh oh!

JBludau Feb 1, 2026

Uh oh!

JBludau Feb 1, 2026

Uh oh!

JBludau Feb 2, 2026

Uh oh!

yaowwwwwww commented Feb 2, 2026

Uh oh!

JBludau commented Feb 2, 2026 •

edited

Loading

Uh oh!

JBludau left a comment

Uh oh!

JBludau Apr 29, 2026

Uh oh!

JBludau Apr 29, 2026

Uh oh!

JBludau Apr 29, 2026

Uh oh!

JBludau Apr 29, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

Conversation

JBludau commented Jan 28, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

yaowwwwwww commented Feb 2, 2026

Uh oh!

JBludau commented Feb 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

JBludau left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants

JBludau commented Jan 28, 2026 •

edited

Loading

JBludau commented Feb 2, 2026 •

edited

Loading