exp contact optimizations (#29732)

framework/include/timeintegrators/ExplicitTimeIntegrator.h

@@ -72,7 +72,7 @@ class ExplicitTimeIntegrator : public TimeIntegrator, public MeshChangedInterfac
   NumericVector<Real> & _solution_update;
 
   /// Diagonal of the lumped mass matrix (and its inversion)
-  NumericVector<Real> & _mass_matrix_diag;
+  NumericVector<Real> & _mass_matrix_diag_inverted;
 
   /// Vector of 1's to help with creating the lumped mass matrix
   NumericVector<Real> * _ones;

framework/src/timeintegrators/ExplicitTimeIntegrator.C

@@ -41,7 +41,7 @@ ExplicitTimeIntegrator::ExplicitTimeIntegrator(const InputParameters & parameter
     _solve_type(getParam<MooseEnum>("solve_type")),
     _explicit_residual(_nl.addVector("explicit_residual", false, PARALLEL)),
     _solution_update(_nl.addVector("solution_update", true, PARALLEL)),
-    _mass_matrix_diag(_nl.addVector("mass_matrix_diag", false, PARALLEL))
+    _mass_matrix_diag_inverted(_nl.addVector("mass_matrix_diag_inverted", true, PARALLEL))
 {
   _Ke_time_tag = _fe_problem.getMatrixTagID("TIME");
 
@@ -50,7 +50,7 @@ ExplicitTimeIntegrator::ExplicitTimeIntegrator(const InputParameters & parameter
   _fe_problem.solverParams()._type = Moose::ST_LINEAR;
 
   if (_solve_type == LUMPED || _solve_type == LUMP_PRECONDITIONED)
-    _ones = &_nl.addVector("ones", false, PARALLEL);
+    _ones = &_nl.addVector("ones", true, PARALLEL);
 }
 
 void
@@ -84,7 +84,7 @@ ExplicitTimeIntegrator::meshChanged()
 
   if (_solve_type == LUMP_PRECONDITIONED)
   {
-    _preconditioner = std::make_unique<LumpedPreconditioner>(_mass_matrix_diag);
+    _preconditioner = std::make_unique<LumpedPreconditioner>(_mass_matrix_diag_inverted);
     _linear_solver->attach_preconditioner(_preconditioner.get());
     _linear_solver->init();
   }
@@ -111,13 +111,13 @@ ExplicitTimeIntegrator::performExplicitSolve(SparseMatrix<Number> & mass_matrix)
       // Computes the sum of each row (lumping)
       // Note: This is actually how PETSc does it
       // It's not "perfectly optimal" - but it will be fast (and universal)
-      mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
+      mass_matrix.vector_mult(_mass_matrix_diag_inverted, *_ones);
 
       // "Invert" the diagonal mass matrix
-      _mass_matrix_diag.reciprocal();
+      _mass_matrix_diag_inverted.reciprocal();
 
       // Multiply the inversion by the RHS
-      _solution_update.pointwise_mult(_mass_matrix_diag, _explicit_residual);
+      _solution_update.pointwise_mult(_mass_matrix_diag_inverted, _explicit_residual);
 
       // Check for convergence by seeing if there is a nan or inf
       auto sum = _solution_update.sum();
@@ -130,8 +130,8 @@ ExplicitTimeIntegrator::performExplicitSolve(SparseMatrix<Number> & mass_matrix)
     }
     case LUMP_PRECONDITIONED:
     {
-      mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
-      _mass_matrix_diag.reciprocal();
+      mass_matrix.vector_mult(_mass_matrix_diag_inverted, *_ones);
+      _mass_matrix_diag_inverted.reciprocal();
 
       converged = solveLinearSystem(mass_matrix);
 

modules/contact/src/constraints/ExplicitDynamicsContactConstraint.C

@@ -318,14 +318,14 @@ ExplicitDynamicsContactConstraint::solveImpactEquations(const Node & node,
   auto & u_dot = *_sys.solutionUDot();
 
   // Get lumped mass value
-  auto & diag = _sys.getVector("mass_matrix_diag");
-  Real mass_q = diag(dof_x);
+  auto & diag = _sys.getVector("mass_matrix_diag_inverted");
+  Real mass_inverted = diag(dof_x);
 
   // Include effects of other forces:
   // Initial guess: v_{n-1/2} + dt * M^{-1} * (F^{ext} - F^{int})
-  Real velocity_x = u_dot(dof_x) + _dt / mass_q * -1 * _residual_copy(dof_x);
-  Real velocity_y = u_dot(dof_y) + _dt / mass_q * -1 * _residual_copy(dof_y);
-  Real velocity_z = u_dot(dof_z) + _dt / mass_q * -1 * _residual_copy(dof_z);
+  Real velocity_x = u_dot(dof_x) + _dt * mass_inverted * -1 * _residual_copy(dof_x);
+  Real velocity_y = u_dot(dof_y) + _dt * mass_inverted * -1 * _residual_copy(dof_y);
+  Real velocity_z = u_dot(dof_z) + _dt * mass_inverted * -1 * _residual_copy(dof_z);
 
   Real n_velocity_x = _neighbor_vel_x[0];
   Real n_velocity_y = _neighbor_vel_y[0];
@@ -337,56 +337,19 @@ ExplicitDynamicsContactConstraint::solveImpactEquations(const Node & node,
       n_velocity_x, n_velocity_y, _mesh.dimension() == 3 ? n_velocity_z : 0.0);
   gap_rate = pinfo->_normal * (secondary_velocity - closest_point_velocity);
 
-  // Prepare equilibrium loop
-  bool is_converged(false);
-  unsigned int iteration_no(0);
-  const unsigned int max_no_iterations(20000);
+  // Compute the force increment needed
+  RealVectorValue impulse_force = pinfo->_normal * (gap_rate / mass_inverted) / _dt;
+  pinfo->_contact_force = impulse_force;
 
-  // Initialize augmented iteration variable
-  Real gap_rate_old(0.0);
-  Real force_increment(0.0);
-  Real force_increment_old(0.0);
-  Real lambda_iteration(0);
-
-  while (!is_converged && iteration_no < max_no_iterations)
-  {
-    // Start a loop until we converge on normal contact forces
-    gap_rate_old = gap_rate;
-    gap_rate = pinfo->_normal * (secondary_velocity - closest_point_velocity);
-    force_increment_old = force_increment;
-
-    force_increment = mass_contact_pressure * gap_rate;
-
-    velocity_x -= _dt / mass_q * (pinfo->_normal(0) * (force_increment));
-    velocity_y -= _dt / mass_q * (pinfo->_normal(1) * (force_increment));
-    velocity_z -= _dt / mass_q * (pinfo->_normal(2) * (force_increment));
-
-    // Let's not modify the neighbor velocity, but apply the corresponding force.
-    // TODO: Update for multi-body impacts
-    // n_velocity_x = n_velocity_x;
-    // n_velocity_y = n_velocity_y;
-    // n_velocity_z = n_velocity_z;
-
-    secondary_velocity = {velocity_x, velocity_y, _mesh.dimension() == 3 ? velocity_z : 0.0};
-    closest_point_velocity = {
-        n_velocity_x, n_velocity_y, _mesh.dimension() == 3 ? n_velocity_z : 0.0};
-
-    // Convergence check
-    lambda_iteration += force_increment;
-
-    const Real relative_error = (force_increment - force_increment_old) / force_increment;
-    const Real absolute_error = std::abs(force_increment);
-
-    if (std::abs(relative_error) < TOLERANCE * TOLERANCE || absolute_error < TOLERANCE ||
-        (gap_rate_old) * (gap_rate) < 0.0)
-      is_converged = true;
-    else
-      iteration_no++;
-  }
+  // recalculate velocity to determine gap rate
+  velocity_x -= _dt * mass_inverted * impulse_force(0);
+  velocity_y -= _dt * mass_inverted * impulse_force(1);
+  velocity_z -= _dt * mass_inverted * impulse_force(2);
 
+  secondary_velocity = {velocity_x, velocity_y, _mesh.dimension() == 3 ? velocity_z : 0.0};
+  gap_rate = pinfo->_normal * (secondary_velocity - closest_point_velocity);
+  // gap rate is now always machine precision
   _gap_rate->setNodalValue(gap_rate);
-
-  pinfo->_contact_force = pinfo->_normal * lambda_iteration;
 }
 
 Real

modules/contact/test/tests/explicit_dynamics/tests

@@ -1,5 +1,5 @@
 [Tests]
-  issues = '#25666 #28008'
+  issues = '#25666 #28008 #29732'
   design = 'ExplicitDynamicsContactConstraint.md'
   [block_penalty]
     type = 'CSVDiff'
@@ -17,7 +17,7 @@
     exodiff = 'test_balance_out.e'
     abs_zero = 1.0e-4
     allow_warnings = true
-    method = '!dbg'                  
+    method = '!dbg'
     requirement = 'The system shall be able to solve a simple few-element normal contact problem '
                   'using explicit dynamics solving uncoupled, local equations of momentum balance.'
   []
@@ -27,7 +27,7 @@
     exodiff = 'test_balance_short_out.e'
     abs_zero = 1.0e-4
     allow_warnings = true
-    cli_args = 'Outputs/file_base=test_balance_short_out Executioner/end_time=0.001'            
+    cli_args = 'Outputs/file_base=test_balance_short_out Executioner/end_time=0.001'
     requirement = 'The system shall be able to solve a simple few-element normal contact problem '
                   'using explicit dynamics solving uncoupled, local equations of momentum balance in debug mode.'
   []
@@ -53,4 +53,15 @@
                   'during an impact-settling under increased gravity acceleration.'
     heavy = true
   []
+  [exp_constant_mass]
+    type = 'Exodiff'
+    input = 'test_balance.i'
+    exodiff = 'test_balance_out.e'
+    abs_zero = 1.0e-4
+    allow_warnings = true
+    method = '!dbg'
+    cli_args ='Executioner/TimeIntegrator/use_constant_mass=true Mesh/patch_update_strategy=always'
+    requirement = 'The system shall be able to solve an explicit contact problem'
+    'with the constant mass option and with updating the contact patch often.'
+  []
 []

modules/solid_mechanics/src/bcs/DirectDirichletBCBase.C

@@ -19,7 +19,7 @@ DirectDirichletBCBase::validParams()
 
 DirectDirichletBCBase::DirectDirichletBCBase(const InputParameters & parameters)
   : NodalBC(parameters),
-    _mass_diag(_sys.getVector("mass_matrix_diag")),
+    _mass_diag(_sys.getVector("mass_matrix_diag_inverted")),
     _u_old(_var.nodalValueOld()),
     _u_dot_old(_var.nodalValueDotOld())
 {
@@ -35,7 +35,7 @@ DirectDirichletBCBase::computeQpResidual()
   // This is the force required to enforce the BC in a central difference scheme
   Real avg_dt = (_dt + _dt_old) / 2;
   Real resid = (computeQpValue() - _u_old) / (avg_dt * _dt) - (_u_dot_old) / avg_dt;
-  resid *= -_mass_diag(dofnum);
+  resid /= -_mass_diag(dofnum);
 
   return resid;
 }

modules/solid_mechanics/src/timeintegrators/DirectCentralDifference.C

@@ -81,17 +81,23 @@ DirectCentralDifference::solve()
   auto & mass_matrix = _nonlinear_implicit_system->get_system_matrix();
 
   // Compute the mass matrix
-  if (!_constant_mass || (_constant_mass && _t_step == 1))
+  if (!_constant_mass || _t_step == 1)
+  {
+    // We only want to compute "inverted" lumped mass matrix once.
     _fe_problem.computeJacobianTag(
         *_nonlinear_implicit_system->current_local_solution, mass_matrix, mass_tag);
 
+    // Calculating the lumped mass matrix for use in residual calculation
+    mass_matrix.vector_mult(_mass_matrix_diag_inverted, *_ones);
+
+    // "Invert" the diagonal mass matrix
+    _mass_matrix_diag_inverted.reciprocal();
+  }
+
   // Set time to the time at which to evaluate the residual
   _fe_problem.time() = _fe_problem.timeOld();
   _nonlinear_implicit_system->update();
 
-  // Calculating the lumped mass matrix for use in residual calculation
-  mass_matrix.vector_mult(_mass_matrix_diag, *_ones);
-
   // Compute the residual
   _explicit_residual.zero();
   _fe_problem.computeResidual(
@@ -130,11 +136,9 @@ DirectCentralDifference::performExplicitSolve(SparseMatrix<Number> &)
 {
   bool converged = false;
 
-  // "Invert" the diagonal mass matrix
-  _mass_matrix_diag.reciprocal();
   // Calculate acceleration
   auto & accel = *_sys.solutionUDotDot();
-  accel.pointwise_mult(_mass_matrix_diag, _explicit_residual);
+  accel.pointwise_mult(_mass_matrix_diag_inverted, _explicit_residual);
 
   // Scaling the acceleration
   auto accel_scaled = accel.clone();