Add CFL limit to timestep constraints

docs/common_plane_contact_mass_timestep_plan.md

@@ -0,0 +1,75 @@
+# Plan: Common Plane Contact Mass-Based Timestep Limit
+
+## Summary
+
+Refactor the `COMMON_PLANE + PENALTY` timestep stability limit from a host-`dt` scaling into an absolute contact timestep estimate based on host-provided **lumped nodal masses**. Implement the core registration path in `tribol.hpp` and add an MFEM convenience wrapper in `mfem_tribol.hpp`.
+
+This work should be implemented **after** higher-order Common Plane quadrature is added, since both changes touch overlap integration, contact weighting, and how overlap contributions are accumulated.
+
+## Key Changes
+
+- Replace the current CFL-like stability vote with an absolute contact timestep for active Common Plane penalty contact:
+  \[
+  \Delta t_{\text{contact}} = \alpha \, 2 \sqrt{\frac{m_{\text{pair}}}{k_{\text{pair}}}}
+  \]
+  where `alpha = timestep_scale`, `k_pair = k_c * overlap_area`, and `k_c` is the existing interface stiffness per area.
+- Compute side effective mass from **lumped nodal masses** using the same contact weights already used to distribute force:
+  \[
+  m_{i,\text{eff}} = \left(\sum_a \frac{\phi_{i,a}^2}{m_{i,a}}\right)^{-1}
+  \]
+  using the overlap-based weak-form weights for each side.
+- Combine the two sides into a pair effective mass with harmonic reduction:
+  \[
+  m_{\text{pair}} = \left(\frac{1}{m_{1,\text{eff}}} + \frac{1}{m_{2,\text{eff}}}\right)^{-1}
+  \]
+- Do not add extra overlap-area fraction scaling for mass in v1. The overlap localization is handled through the contact weights.
+- Keep the existing gap-based and velocity-projection timestep checks unchanged. Only the stability-limit bucket changes from a `dt_in` scaler to an absolute contact vote.
+
+## Public API / Interface Changes
+
+- Add a core registration function in `tribol.hpp` / `tribol.cpp`:
+  - `registerLumpedNodalMass( IndexT mesh_id, const RealT* lumped_nodal_mass )`
+- Add storage and validation for lumped nodal mass in mesh nodal data.
+  - Scalar field, length = number of contact surface nodes on the registered mesh.
+  - All values must be strictly positive.
+- Add an MFEM convenience wrapper in `mfem_tribol.hpp` / `mfem_tribol.cpp`:
+  - `registerMfemLumpedNodalMass( IndexT cs_id, const mfem::ParGridFunction& lumped_mass )`
+- MFEM helper behavior:
+  - Accept a scalar nodal `ParGridFunction` on the parent mesh.
+  - Transfer/project it onto the Tribol contact surface nodes using the existing parent-to-submesh machinery.
+- Validation rule:
+  - If `enableTimestepStabilityLimits()` is on for a Common Plane penalty scheme, lumped nodal mass must be registered and valid on both meshes.
+  - Missing or nonpositive nodal mass makes the coupling scheme invalid during `init()`, with a warning, matching existing invalid-enforcement-data behavior.
+
+## Implementation Notes
+
+- Reuse the existing overlap weights for the effective-mass calculation so the inertial reduction follows the same contact discretization as the force calculation.
+- Use the existing stiffness path:
+  - `KINEMATIC_CONSTANT`: `k_i = pen_scale * const_penalty`
+  - `KINEMATIC_ELEMENT`: `k_i = pen_scale * mat_mod / thickness`
+  - `k_c = (k1 * k2) / (k1 + k2)`
+  - `k_pair = k_c * plane.m_area`
+- Guard all effective-mass calculations against zero or invalid mass before division.
+- Keep the existing invalid-thickness early return in Common Plane penalty.
+- Do not support consistent mass in v1 and do not implement internal mass lumping. The host is responsible for providing already-lumped nodal masses.
+- When higher-order quadrature lands, generalize the single-point weight usage to a quadrature-point accumulation rather than reintroducing a single-point assumption.
+
+## Test Plan
+
+- Add/update core tests for Common Plane penalty:
+  - Constant-penalty contact with registered nodal masses: verify absolute contact timestep matches `2*sqrt(m_pair/k_pair)` times `timestep_scale`.
+  - Element-penalty contact with registered nodal masses: same verification using element-based stiffness.
+  - Stationary contact with stability limits enabled: verify the absolute contact limit is computed from mass/stiffness and is not a function of incoming `dt`.
+  - Missing nodal mass with stability limits enabled: verify the coupling scheme is skipped as invalid.
+  - Nonpositive nodal mass with stability limits enabled: verify invalid scheme and no contact response.
+- Add MFEM-path coverage:
+  - Register MFEM lumped nodal mass through the new helper and verify the resulting timestep vote matches the core formula on a simple contact case.
+- Re-run existing `tribol_timestep_vote` and `tribol_common_plane_penalty` suites to ensure the penetration-based votes still behave the same.
+
+## Assumptions
+
+- V1 supports only **lumped nodal masses** as inertial input.
+- V1 computes contact mass from contact weights and does not apply extra overlap-fraction scaling.
+- V1 adds both the core Tribol API support and an MFEM convenience registration helper.
+- `timestep_scale` remains the user-facing safety factor and multiplies the absolute contact critical timestep.
+- The stability constant is `2`, corresponding to the lumped-mass single-DOF spring stability bound.

src/tests/tribol_common_plane_penalty.cpp

@@ -173,6 +173,26 @@ void checkForceSense( tribol::CouplingScheme const* cs, bool isTied = false )
   }
 }  // end checkForceSense()
 
+void checkZeroResponses( tribol::CouplingScheme const* cs, const RealT tol = 1.e-14 )
+{
+  const auto mesh1 = const_cast<tribol::CouplingScheme*>( cs )->getMesh1().getView();
+  const auto mesh2 = const_cast<tribol::CouplingScheme*>( cs )->getMesh2().getView();
+
+  for ( int i = 0; i < 2; ++i ) {
+    auto& mesh = ( i == 0 ) ? mesh1 : mesh2;
+    for ( tribol::IndexT kf = 0; kf < mesh.numberOfElements(); ++kf ) {
+      for ( tribol::IndexT a = 0; a < mesh.numberOfNodesPerElement(); ++a ) {
+        int node_id = mesh.getGlobalNodeId( kf, a );
+        EXPECT_NEAR( mesh.getResponse()[0][node_id], 0., tol );
+        EXPECT_NEAR( mesh.getResponse()[1][node_id], 0., tol );
+        if ( mesh.spatialDimension() == 3 ) {
+          EXPECT_NEAR( mesh.getResponse()[2][node_id], 0., tol );
+        }
+      }
+    }
+  }
+}  // end checkZeroResponses()
+
 /*!
  * Test fixture class with some setup necessary to test
  * the COMMON_PLANE + PENALTY implementation
@@ -449,6 +469,76 @@ TEST_F( CommonPlaneTest, element_penalty_check )
   tribol::finalize();
 }
 
+TEST_F( CommonPlaneTest, element_penalty_nonpositive_thickness_returns_error_without_forces )
+{
+  this->m_mesh.mortarMeshId = 0;
+  this->m_mesh.nonmortarMeshId = 1;
+
+  int nMortarElems = 4;
+  int nElemsXM = nMortarElems;
+  int nElemsYM = nMortarElems;
+  int nElemsZM = nMortarElems;
+
+  int nNonmortarElems = 5;
+  int nElemsXS = nNonmortarElems;
+  int nElemsYS = nNonmortarElems;
+  int nElemsZS = nNonmortarElems;
+
+  RealT x_min1 = 0.;
+  RealT y_min1 = 0.;
+  RealT z_min1 = 0.;
+  RealT x_max1 = 1.;
+  RealT y_max1 = 1.;
+  RealT z_max1 = 1.05;
+
+  RealT x_min2 = 0.;
+  RealT y_min2 = 0.;
+  RealT z_min2 = 0.95;
+  RealT x_max2 = 1.;
+  RealT y_max2 = 1.;
+  RealT z_max2 = 2.;
+
+  RealT element_thickness2 = ( z_max2 - z_min2 ) / nElemsZS;
+
+  this->m_mesh.setupContactMeshHex( nElemsXM, nElemsYM, nElemsZM, x_min1, y_min1, z_min1, x_max1, y_max1, z_max1,
+                                    nElemsXS, nElemsYS, nElemsZS, x_min2, y_min2, z_min2, x_max2, y_max2, z_max2, 0.,
+                                    0. );
+
+  RealT bulk_mod1 = 1.0;
+  RealT bulk_mod2 = 1.0;
+  RealT velX1 = 0.;
+  RealT velY1 = 0.;
+  RealT velZ1 = 0.;
+  RealT velX2 = 0.;
+  RealT velY2 = 0.;
+  RealT velZ2 = 0.;
+
+  this->m_mesh.allocateAndSetVelocities( m_mesh.mortarMeshId, velX1, velY1, velZ1 );
+  this->m_mesh.allocateAndSetVelocities( m_mesh.nonmortarMeshId, velX2, velY2, -velZ2 );
+
+  this->m_mesh.allocateAndSetElementThickness( m_mesh.mortarMeshId, 0. );
+  this->m_mesh.allocateAndSetBulkModulus( m_mesh.mortarMeshId, bulk_mod1 );
+  this->m_mesh.allocateAndSetElementThickness( m_mesh.nonmortarMeshId, element_thickness2 );
+  this->m_mesh.allocateAndSetBulkModulus( m_mesh.nonmortarMeshId, bulk_mod2 );
+
+  tribol::TestControlParameters parameters;
+  parameters.penalty_ratio = true;
+  parameters.const_penalty = 0.75;
+  parameters.dt = 1.e-3;
+
+  int test_mesh_update_err = this->m_mesh.tribolSetupAndUpdate(
+      tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters );
+
+  EXPECT_EQ( test_mesh_update_err, 0 );
+
+  tribol::CouplingSchemeManager& couplingSchemeManager = tribol::CouplingSchemeManager::getInstance();
+  tribol::CouplingScheme* couplingScheme = &couplingSchemeManager.at( 0 );
+  EXPECT_FALSE( couplingScheme->init() );
+  checkZeroResponses( couplingScheme );
+
+  tribol::finalize();
+}
+
 TEST_F( CommonPlaneTest, tied_contact_check )
 {
   this->m_mesh.mortarMeshId = 0;

src/tests/tribol_timestep_vote.cpp

@@ -305,7 +305,7 @@ TEST_F( CommonPlaneTest, numerically_zero_velocity_small_gap )
 
   EXPECT_EQ( test_mesh_update_err, 0 );
 
-  // expect no change in dt
+  // expect no change in dt because velocity is numerically zero and stability limits are disabled by default
   EXPECT_EQ( parameters.dt, dt );
 
   tribol::finalize();
@@ -390,10 +390,7 @@ TEST_F( CommonPlaneTest, zero_velocity_large_gap )
       tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters );
 
   EXPECT_EQ( test_mesh_update_err, 0 );
-  // note that with very small velocity, the dt estimate will be
-  // very large, and won't change the timestep, even if the gap is large. This
-  // allows for a soft contact response in the presence of a small velocity that
-  // won't actually correct too much interpen with a contact dt vote
+  // expect no change in dt because velocity is zero and stability limits are disabled by default
   EXPECT_EQ( parameters.dt, dt );
 
   tribol::finalize();
@@ -652,7 +649,7 @@ TEST_F( CommonPlaneTest, separation_velocity_small_gap )
 
   EXPECT_EQ( test_mesh_update_err, 0 );
 
-  // no change in dt because of separation velocities
+  // expect no change in dt because velocities are separating and stability limits are disabled by default
   EXPECT_EQ( parameters.dt, dt );
 
   tribol::finalize();
@@ -917,6 +914,100 @@ TEST_F( CommonPlaneTest, large_velocity_small_separation_set_alpha )
   tribol::finalize();
 }
 
+TEST_F( CommonPlaneTest, critical_timestep_cfl_limit )
+{
+  // This test isolates and validates the Courant (CFL) stability limit.
+  // By setting a small timestep_scale (alpha = 0.1), we ensure that the CFL limit
+  // becomes the governing minimum timestep limit.
+  this->m_mesh.mortarMeshId = 0;
+  this->m_mesh.nonmortarMeshId = 1;
+
+  constexpr int nMortarElems = 4;
+  constexpr int nElemsXM = nMortarElems;
+  constexpr int nElemsYM = nMortarElems;
+  constexpr int nElemsZM = nMortarElems;
+
+  constexpr int nNonmortarElems = 5;
+  constexpr int nElemsXS = nNonmortarElems;
+  constexpr int nElemsYS = nNonmortarElems;
+  constexpr int nElemsZS = nNonmortarElems;
+
+  constexpr RealT x_min1 = 0.;
+  constexpr RealT y_min1 = 0.;
+  constexpr RealT z_min1 = 0.;
+  constexpr RealT x_max1 = 1.;
+  constexpr RealT y_max1 = 1.;
+  constexpr RealT z_max1 = 1.005;
+
+  constexpr RealT x_min2 = 0.;
+  constexpr RealT y_min2 = 0.;
+  constexpr RealT z_min2 = 0.95;
+  constexpr RealT x_max2 = 1.;
+  constexpr RealT y_max2 = 1.;
+  constexpr RealT z_max2 = 2.;
+
+  constexpr RealT element_thickness1 = ( z_max1 - z_min1 ) / nElemsZM;
+  constexpr RealT element_thickness2 = ( z_max2 - z_min2 ) / nElemsZS;
+
+  this->m_mesh.setupContactMeshHex( nElemsXM, nElemsYM, nElemsZM, x_min1, y_min1, z_min1, x_max1, y_max1, z_max1,
+                                    nElemsXS, nElemsYS, nElemsZS, x_min2, y_min2, z_min2, x_max2, y_max2, z_max2, 0.,
+                                    0. );
+
+  // Stationary contact to bypass impact limit
+  constexpr RealT dt = 1.0;
+  RealT bulk_mod1 = 1.0;
+  RealT bulk_mod2 = 1.0;
+  RealT velX1 = 0.;
+  RealT velY1 = 0.;
+  RealT velZ1 = 0.;
+  RealT velX2 = 0.;
+  RealT velY2 = 0.;
+  RealT velZ2 = 0.;
+
+  this->m_mesh.allocateAndSetVelocities( m_mesh.mortarMeshId, velX1, velY1, velZ1 );
+  this->m_mesh.allocateAndSetVelocities( m_mesh.nonmortarMeshId, velX2, velY2, -velZ2 );
+
+  this->m_mesh.allocateAndSetElementThickness( m_mesh.mortarMeshId, element_thickness1 );
+  this->m_mesh.allocateAndSetBulkModulus( m_mesh.mortarMeshId, bulk_mod1 );
+  this->m_mesh.allocateAndSetElementThickness( m_mesh.nonmortarMeshId, element_thickness2 );
+  this->m_mesh.allocateAndSetBulkModulus( m_mesh.nonmortarMeshId, bulk_mod2 );
+
+  tribol::TestControlParameters parameters;
+  parameters.penalty_ratio = true;
+  parameters.const_penalty = 0.75;
+  parameters.dt = dt;
+  parameters.enable_timestep_vote = true;
+  parameters.enable_timestep_stability_limits = true;
+  parameters.timestep_pen_frac = 0.3;
+  parameters.timestep_scale = 0.1;  // alpha = 0.1
+
+  int test_mesh_update_err = this->m_mesh.tribolSetupAndUpdate(
+      tribol::COMMON_PLANE, tribol::PENALTY, tribol::FRICTIONLESS, tribol::NO_CASE, false, parameters );
+
+  EXPECT_EQ( test_mesh_update_err, 0 );
+
+  // Expected CFL value:
+  // Stiffness under KINEMATIC_ELEMENT is computed as: bulk_modulus / thickness
+  // for each contacting face, where bulk_modulus is 1.0.
+  // The contact interface spring is modeled as mesh 1 and mesh 2's springs in series:
+  // penalty_stiff_per_area = (stiffness1 * stiffness2) / (stiffness1 + stiffness2)
+  constexpr RealT stiffness1 = 1.0 / element_thickness1;
+  constexpr RealT stiffness2 = 1.0 / element_thickness2;
+  constexpr RealT penalty_stiff_per_area = ( stiffness1 * stiffness2 ) / ( stiffness1 + stiffness2 );
+
+  // CFL scaling factors based on added interface stiffness:
+  // f_scale = 1 / sqrt( 1 + penalty_stiff_per_area / stiffness )
+  const RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) );
+  const RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) );
+
+  // dt_CFL = alpha * dt * min( f_scale1, f_scale2 )
+  const RealT expected_dt = parameters.timestep_scale * dt * std::min( f_scale1, f_scale2 );
+  RealT dt_tol = 1.e-6;
+  EXPECT_NEAR( parameters.dt, expected_dt, dt_tol );
+
+  tribol::finalize();
+}
+
 int main( int argc, char* argv[] )
 {
   int result = 0;

src/tribol/common/Parameters.hpp

@@ -483,6 +483,8 @@ struct Parameters {
   int vis_cycle_incr = 100;           ///! Frequency for visualizations dumps
   VisType vis_type = VIS_NONE;        ///! Type of interface physics visualization output
   bool enable_timestep_vote = false;  ///! True if host-code desires the timestep vote to be calculated and returned
+  bool enable_timestep_stability_limits =
+      false;  ///! True if host-code desires dynamic contact stability limits to be included in timestep vote
 
   bool auto_interpen_check = false;  ///! True if the auto-contact interpenetration check is used for full-overlap pairs
 

src/tribol/interface/tribol.cpp

@@ -350,6 +350,19 @@ void enableTimestepVote( IndexT cs_id, const bool enable )
 
 }  // end enableTimestepVote()
 
+//------------------------------------------------------------------------------
+void enableTimestepStabilityLimits( IndexT cs_id, const bool enable )
+{
+  auto cs = CouplingSchemeManager::getInstance().findData( cs_id );
+
+  // check to see if coupling scheme exists
+  SLIC_ERROR_ROOT_IF( !cs, "tribol::enableTimestepStabilityLimits(): call tribol::registerCouplingScheme() "
+                               << "prior to calling this routine." );
+
+  cs->getParameters().enable_timestep_stability_limits = enable;
+
+}  // end enableTimestepStabilityLimits()
+
 //------------------------------------------------------------------------------
 void enableEnzyme( IndexT cs_id, [[maybe_unused]] bool use_enzyme )
 {

src/tribol/interface/tribol.hpp

@@ -226,6 +226,17 @@ void setBinningProximityScale( IndexT cs_id, RealT binning_proximity_scale );
  */
 void enableTimestepVote( IndexT cs_id, const bool enable );
 
+/*!
+ * \brief Enable CFL-like contact stiffness limits in the timestep vote
+ *
+ * \param [in] cs_id coupling scheme id
+ * \param [in] enable the contact stiffness limit will be included in the timestep vote if true
+ *
+ * \note default behavior is to not enable the contact stiffness limit
+ *
+ */
+void enableTimestepStabilityLimits( IndexT cs_id, const bool enable );
+
 /**
  * @brief Enable Enzyme AD for Jacobian calculations
  *