Draft: implement super-ellipsoid particle shapes for use in granular simulations

[jibril-b-coulibaly]

Summary

This PR allows super-ellipsoid particle shapes and, contact detection and force calculation for granular simulations.

Related Issue(s)

follow-up #3999

Author(s)

Jibril B. Coulibaly (Sandia National Laboratories)

Licensing

By submitting this pull request, I agree, that my contribution will be included in LAMMPS and redistributed under either the GNU General Public License version 2 (GPL v2) or the GNU Lesser General Public License version 2.1 (LGPL v2.1).

Backward Compatibility

• Super-ellipsoids require 2 additional exponents that will be read from and written to the data file. For read_data, I believe they can be made optional and exponents values default to 2. write_data would always write these 2 additional values, which shouldn't be a problem for LAMMPS. This could change the behavior of some post-processing routine that is not expecting the 2 additional values on the ellipsoid lines
• more TBD

Implementation Notes

• Super-ellipsoid particles are defined by adding blockiness exponents to the ellipsoid atom style.
• The algorithms and methods are mostly based on the paper by Podlozhnyuk et al. 2017 and can/will leverage some of the existing LIGGGHTS implementation.

Post Submission Checklist

• The feature or features in this pull request is complete
• Licensing information is complete
• Corresponding author information is complete
• The source code follows the LAMMPS formatting guidelines
• Suitable new documentation files and/or updates to the existing docs are included
• The added/updated documentation is integrated and tested with the documentation build system
• The feature has been verified to work with the conventional build system
• The feature has been verified to work with the CMake based build system
• Suitable tests have been added to the unittest tree.
• A package specific README file has been included or updated
• One or more example input decks are included

Further Information, Files, and Links

TBD

[jibril-b-coulibaly]

Hello, I'm running into an obstacle I'd like your insight about.

The volume and moments of inertia of the super-ellipsoid are defined in terms of the Beta function.

The STL special functions are implementation-dependent and std::beta() does not compile on my end: I need to include <tr1/cmath> and call std::tr1::beta() as suggested by the doc. TR1 is from C++03 so it seems like it should work on most implementation, but it does not pass CI. For comparison, the LIGGGHTS implementation checks:

#if __STDCPP_MATH_SPEC_FUNCS__ >= 201003L
    #define BETA_NAMESPACE std
#elif defined(HAVE_TR1_CMATH)
    #include <tr1/cmath>
    #define BETA_NAMESPACE std::tr1
#elif defined(BOOST_INCLUDED)
    #include "boost/math/special_functions/beta.hpp"
    #define BETA_NAMESPACE boost::math

Would a similar check for __STDCPP_MATH_SPEC_FUNCS__ >= 201003L be appropriate inatom_vec_ellipsoid.cpp ? Or would special functions in LAMMPS need to be managed at a higher level, e.g., math_special?

Thank you

[akohlmey]

Would a similar check for __STDCPP_MATH_SPEC_FUNCS__ >= 201003L be appropriate inatom_vec_ellipsoid.cpp ? Or would special functions in LAMMPS need to be managed at a higher level, e.g., math_special?

As far as I can tell std::beta() is part of the C++17 standard. But core code has to be compatible with C++11. For that we have to provide our own implementation. We prefer to keep the code itself as clean as possible, so a wrapper function is the way to go. The MathSpecial namespace would be the best match. Please also note, that we (very deliberately) do not use Boost.

[jibril-b-coulibaly]

For that we have to provide our own implementation. We prefer to keep the code itself as clean as possible, so a wrapper function is the way to go.

For the wrapper function, do I understand correctly that I would need to write a copy of an existing implementation, e.g., the way erfcx_y100 is implemented in math_special ?

[stanmoore1]

@jibril-b-coulibaly seeing some compile errors with this PR:

/mnt/lammps/workspace/dev/pull_requests/compilation_tests/lammps/src/atom_vec_ellipsoid.cpp: In member function 'virtual void LAMMPS_NS::AtomVecEllipsoid::pack_data_pre(int)':
/mnt/lammps/workspace/dev/pull_requests/compilation_tests/lammps/src/atom_vec_ellipsoid.cpp:513:5: error: 'block' was not declared in this scope
  513 |     block = bonus[ellipsoid_flag].block;
      |     ^~~~~
/mnt/lammps/workspace/dev/pull_requests/compilation_tests/lammps/src/atom_vec_ellipsoid.cpp:514:5: error: 'flag_super' was not declared in this scope
  514 |     flag_super = bonus[ellipsoid_flag].flag_super;
      |     ^~~~~~~~~~
/mnt/lammps/workspace/dev/pull_requests/compilation_tests/lammps/src/atom_vec_ellipsoid.cpp: In member function 'double LAMMPS_NS::AtomVecEllipsoid::compute_radcirc(double*, double*, bool)':