[Feature Request] Calculate and output the energy lost in normal damper in granular simulations

[vikaskok]

Summary

Energy budgeting is crucial in many simulations. Please provide a way to calculate and obtain energy lost in the damper during granular interactions.

Detailed Description
Granular materials are dissipative. The information of the energy dissipated over time during a simulation will help us better analyse the physics of the problem at hand.

I have looked into the source code (esp. GRANULAR/pair_granular.cpp) and I believe that the following can be done:

1. Keep track of overlap between the particle in the previous time step: deltaprev_ij
   This can possible done by storing delta_ij for all possible combinations (natoms x natoms array) to be used in the next time step.
   Or probably the positions of the particles in previous step are stored double **xold (natoms x 3 array) and deltaprev_ij is calculated in next step using xold array.
2. At each time step, for each particle-particle contact calculate the energy lost in the damper "Fdamp*(delta-deltaprev_ij)" and add them up for all contacts.
   Initiate a global scalar "edamplost" when the simulation begins.
   In the function void PairGranular::compute(int eflag, int vflag) of pair_granular.cpp, we may keep summing up the energy lost as follows
   edamplost += Fdamp*(delta-deltaprev_ij)
3. Energy is lost if particles collide with simulation boundary box or with other walls. Hence, similar changes have to be probably made in fix_wall_gran.cpp and fix_wall_gran_region.cpp
4. Output this information via any fix or compute when the user request.
5. Energy lost due to frictional slider can perhaps also be incorporated.
   However, I think if we can get just the energy lost in the damper, then energy lost due to friction can be found by the equation:
   etotal_init - etotal_final = elost = elostdamper + elostfriction
   (Note: I am ignoring tangential damper for this discussion).

Further Information, Files, and Links
Storing the particle positions in the previous step has been discussed here: https://docs.lammps.org/Developer_write_fix.html
So some part of that code may perhaps be used to store particle positions in previous position. Or perhaps functionality like fix ave/time might help.

[akohlmey]

natoms x natoms array

Already arrays of dimension "natoms" are usually a no-go items for adding to LAMMPS, since they prohibit running very large systems and usually require collective operations which interfere with good parallel scaling. Arrays of dimension "natoms x natoms" are completely out of the question.

You would have to come up with a feature proposal that can properly run with distributed data through the domain decomposition of LAMMPS and avoids such global operations or else it will not be implemented.

[vikaskok]

@akohlmey: Thank you very much for your inputs.
version: lammps-stable-03nov2022
The proposal is as follows:

1. We shall need a fix that is capable of storing positions of all particles in the previous timestep (double **xold; natoms x 3 size). LAMMPS code for such a fix (FixSavePos) already exists [here]((https://docs.lammps.org/Developer_write_fix.html). We need to integrate it into lammps.
2. We shall need a global scalar (edamp) which can store the energy expended in the damper. This shall be initiated with a value 0.0
3. In the pair_granular.cpp file, we calculate overlap delta at each time step (line 329). We also calculate the Force in the damper Fdamp (line 378).
   Here, we can call the fix to retrieve xold array and calculate partilce overlaps in the previous timestep, deltaprev.
   We shall then increment edamp by the energy lost in that timestep for that collision. Since energy is lost in the damper irrespective of whether particles are approaching each other or retreating, we shall use absolute value as follows
   edamp += fabs(Fdamp*(delta - deltaprev));
4. We shall keep accumulating edamp and print it out (perhaps via FixSavePos of step1 above?) as and when the user requests the value.
5. Similar changes have to be made in fix_wall_gran.cpp.
   More specifically, particle wall overlap delta (line 1172), and damping force from the wall, Fdamp (line 1217) is readily available. Below this, we can call call FixSavePos to retrieve xold, calculate deltaprev (overlap with walls), and accumulate the energy lost due to collision with the wall in edamp in the same manner we did before:
   edamp += fabs(Fdamp*(delta - deltaprev));
6. I expect the above can properly run with distributed data through the domain decomposition of LAMMPS since we are handling it just like delta, which I believe is already optimised to exploit the parallel computing capabilites of LAMMPS.

Thank you.