Wrong forces input for bulk stress tensor calculation in pair granular/hooke/history

LAMMPS LAMMPS Development
, ,
1

Dear community,

I have found a miscalculated stress tensor when running LAMMPS (version: 17 Apr 2024) using KOKKOS. After running an Isotropic compression of an FCC granular sample (see Ref. 1 for more information about the input file and system details), I obtained unphysical differences for the diagonal stress components and high values of the off-diagonal ones. The diagonal component must be almost equal, and the off-diagonal ones must be nearly zero.
By delving into the file pair_gran_hooke_history_kokkos.cpp, I identified that the problem is related to the force input in procedures ev_tally_xyz_atom and ev_tally_xyz in PairGranHookeHistoryKokkos::operator(). Both receive the accumulated forces fx_i, fy_i, and fz_i as input, following the code notation; therefore, there is an overaccumulation of stresses. I have addressed the issue using the pair forces fx, fy, and fz, and the outcome results are the expected ones. This is because the stress for every pair force must be computed before the summation is correctly done. This bug does not affect the force calculations but the resulting stress and barostat simulations. Can anyone look at it to see if what I am commenting is right?

Cheers,
Dariel Hernandez

1 - Salomon, J., O’Sullivan, C., & Patino-Ramirez, F. (2024). On data benchmarking and verification of discrete granular simulations. Data in Brief , 53 , 110252.
Data: Validation and Benchmark Dataset for Discrete Element Method Simulations

2

@ddelfin please provide a simple reproducer input deck. @stamoor will have a closer look at that then. Please also consider submitting a bug-report issue on GitHub instead.

2 Likes
3

I am going to submit a bug report issue on GitHub as suggested. You can find my input attached.

Since I wanted to reproduce the results of the cited paper, I also modified the Hooke model to have the Hertz model instead, creating a new pair: pair granular/hertz/history/kk. That is done by multiplying tangential and normal forces by the polyhertz similar to pair_gran_hertz_history.cpp from the GRANULAR package. The .h and .cpp files of the hertz pair are also shared here. Perhaps this pair interaction or a version of it can be deployed in LAMMPS.
In my current version, I have addressed the previously mentioned issue as follows:

my current version:

 if (EVFLAG == 2)
      ev_tally_xyz_atom<NEIGHFLAG, NEWTON_PAIR>(ev, i, j, fx, fy, fz, delx, dely, delz);
    if (EVFLAG == 1)
      ev_tally_xyz<NEWTON_PAIR>(ev, i, j, fx, fy, fz, delx, dely, delz);

Before: 17 Apr 2024 version

 if (EVFLAG == 2)
      ev_tally_xyz_atom<NEIGHFLAG, NEWTON_PAIR>(ev, i, j, fx_i, fy_i, fz_i, delx, dely, delz);
    if (EVFLAG == 1)
      ev_tally_xyz<NEWTON_PAIR>(ev, i, j, fx_i, fy_i, fz_i, delx, dely, delz);

executing LAMMPS

mpirun -np 1 lmp -k on g 1 -sf kk -pk kokkos newton off neigh half comm device gpu/aware on binsize 0.0105-in Isotropic_mu0_0.0.in