Refactor MPM Solver for stress update schemes

[kks32]

Refactor MPM solver to handle different stress update schemes

Summary

Replacing multiple if-else branching conditions for solver types and interface with a new StressUpdate and Interface class. This design will use the Factory pattern to handle different StressUpdate types and Interface.

Motivation

Adding multiple USF/USL/MUSL schemes pollute the MPMExpliciti class with a lot of branching conditions. This new design changes how MPM solvers are written.

Design Detail

The current design has a lot of if-else conditions for checking stress update schemes. USF needs to compute stress first. So it's structure looks like:

    // Update stress first
    if (this->stress_update_ == mpm::StressUpdate::USF)
      this->compute_stress_strain(phase);

Current solver supports usf and usl. When adding new schemes like USAVG and MUSL. This leads to multiple checks in each iteration of the loop. This RFC proposes deriving a new StressUpdate class and create derived classes for different StressUpdate schemes. StressUpdate will have the following design.

template <unsigned Tdim>
class StressUpdate {
 public:
  //! Default constructor with mesh class
  StressUpdate(const std::shared_ptr<mpm::Mesh<Tdim>>& mesh, double dt);

  //! Intialize
  void initialise();

  //! Map mass and momentum to nodes
  void momentum_nodes(unsigned phase);

  //! Compute stress and strain
  //! \param[in] phase Phase to smooth pressure
  //! \param[in] pressure_smoothing Enable or disable pressure smoothing
  void compute_stress_strain(unsigned phase, bool pressure_smoothing);

  //! Precompute stress and strain (empty call)
  //! \param[in] phase Phase to smooth pressure
  //! \param[in] pressure_smoothing Enable or disable pressure smoothing
  virtual void precompute_stress_strain(unsigned phase,
                                        bool pressure_smoothing) = 0;

  //! Postcompute stress and strain (empty call)
  //! \param[in] phase Phase to smooth postssure
  //! \param[in] postssure_smoothing Enable or disable postssure smoothing
  virtual void postcompute_stress_strain(unsigned phase,
                                         bool pressure_smoothing) = 0;

  //! Pressure smoothing
  //! \param[in] phase Phase to smooth pressure
  //! \param[in] pressure_smoothing Enable or disable pressure smoothing
  void pressure_smoothing(unsigned phase);

  //! Compute forces
  //! \param[in] gravity Acceleration due to gravity
  //! \param[in] step Number of step in solver
  //! \param[in] concentrated_nodal_forces Boolean for if a concentrated force
  //! is applied or not
  void compute_forces(const Eigen::Matrix<double, Tdim, 1>& gravity,
                      unsigned phase, unsigned step,
                      bool concentrated_nodal_forces);

  //! Compute acceleration velocity position
  //! \param[in] velocity_update Velocity or acceleration update flag
  //! \param[in] phase Phase of particle
  //! \param[in] damping_type Type of damping
  //! \param[in] damping_factor Value of critical damping
  void compute_particle_kinematics(bool velocity_update, unsigned phase,
                                   const std::string& damping_type,
                                   double damping_factor);

  //! Compute particle location
  //! \param[in] locate_particles Flag to enable locate particles, if set to
  //! false, unlocated particles will be removed
  void locate_particles(bool locate_particles);

  //! Stress update scheme
  //! \retval scheme Stress update scheme
  virtual std::string scheme() const = 0;

The precompute_stress_strain and post_compute_stress_strain will be modified in USF and USL methods. In USF for example, the stress update and scheme output will look like:

//! Precompute stresses and strains
template <unsigned Tdim>
void mpm::StressUpdateUSF<Tdim>::precompute_stress_strain(
    unsigned phase, bool pressure_smoothing) {
  this->compute_stress_strain(phase, pressure_smoothing);
}

//! Postcompute stresses and strains
template <unsigned Tdim>
void mpm::StressUpdateUSF<Tdim>::postcompute_stress_strain(
    unsigned phase, bool pressure_smoothing) {}

//! Stress update scheme
template <unsigned Tdim>
std::string mpm::StressUpdateUSF<Tdim>::scheme() const {
  return "USF";
}

Similar to the StressUpdate design, we will also alter the Interface and InterfaceContact class. InterfaceContact needs a more appropriate name - TBD. Unlike the StressUpdate, the Interface class will have empty void functions, so if no interface model is defined, nothing will be performed.

  //! Intialize
  virtual void initialise(){};

  //! Compute contact forces
  virtual void compute_contact_forces(){};

A factory register will be used to create both StressUpdate scheme and Interface model

Drawbacks

• The performance of this design needs to be checked against the basic version of if-else

• Some additional refactoring may be necessary in the future to split certain function calls.

Rationale and Alternatives

Why is this design the best in the space of possible designs? This design avoids branching, unlike a standard if-else and offers a cleaner interface to add different StressUpdate schemes.

What other designs have been considered and what is the rationale for not choosing them? If-else conditions - ugly and unwieldy to maintain as more stress update schemes are added.

What is the impact of not doing this? Maintaining the MPM Solver becomes tedious.

Unresolved questions

What parts of the design do you expect to resolve through the RFC process before this gets merged?

Check the performance of the code with this new design.

What related issues do you consider out of scope for this RFC that could be addressed in the future independently of the solution that comes out of this RFC?

A draft PR is at: https://github.com/cb-geo/mpm/pull/685

Changelog