From 3d8aa9514dbe0ff7bfb0d0a531a96c7a7ce85f66 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 29 May 2026 17:18:17 -0700 Subject: [PATCH 1/5] add timestep constraints; fix bugs; add tests --- src/tests/tribol_timestep_vote.cpp | 385 ++++++++++++++++++++++++++++- src/tribol/common/Parameters.hpp | 4 + src/tribol/interface/tribol.cpp | 41 +++ src/tribol/interface/tribol.hpp | 3 + src/tribol/mesh/CouplingScheme.cpp | 150 ++++++----- src/tribol/mesh/CouplingScheme.hpp | 24 ++ src/tribol/physics/CommonPlane.cpp | 49 ++-- src/tribol/physics/CommonPlane.hpp | 36 +++ src/tribol/utils/TestUtils.cpp | 3 + src/tribol/utils/TestUtils.hpp | 8 +- 10 files changed, 611 insertions(+), 92 deletions(-) diff --git a/src/tests/tribol_timestep_vote.cpp b/src/tests/tribol_timestep_vote.cpp index 4e2b0f82..e69fd109 100644 --- a/src/tests/tribol_timestep_vote.cpp +++ b/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,382 @@ 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(); +} + +TEST_F( CommonPlaneTest, critical_timestep_chatter_limit ) +{ + // This test validates the Contact Chatter limit calculation. + // Under the default coefficients, the Energy limit (prefactor 0.1118) is always + // 2.81 times more restrictive than the Chatter limit (prefactor 0.3141). + // We verify that the timestep behaves exactly as predicted by the companion limits. + 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. ); + + 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 = 1.0; + parameters.timestep_chatter_factor = 0.01; // gamma = 0.01 + + 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 Chatter limit: + // 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 ); + + // Chatter limit: dt_chatter = pi * gamma * dt * sqrt( stiffness / penalty_stiff_per_area ) + constexpr RealT pi = 3.14159265358979323846; + const RealT dt_chatter1 = + pi * parameters.timestep_chatter_factor * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); + const RealT dt_chatter2 = + pi * parameters.timestep_chatter_factor * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); + const RealT expected_dt = std::min( dt_chatter1, dt_chatter2 ); + RealT dt_tol = 1.e-6; + EXPECT_NEAR( parameters.dt, expected_dt, dt_tol ); + + tribol::finalize(); +} + +TEST_F( CommonPlaneTest, critical_timestep_energy_limit ) +{ + // This test isolates and validates the Spurious Energy Generation limit. + // Under stationary contact (zero velocity), the impact limit is bypassed, + // and the energy limit naturally 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 = 1.0; + + 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 Energy limit 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 ); + + // Energy limit: dt_energy = (sqrt(beta) / 2) * dt * sqrt( stiffness / penalty_stiff_per_area ) + const RealT dt_energy1 = + ( std::sqrt( parameters.timestep_energy_factor ) / 2.0 ) * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); + const RealT dt_energy2 = + ( std::sqrt( parameters.timestep_energy_factor ) / 2.0 ) * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); + const RealT expected_dt = std::min( dt_energy1, dt_energy2 ); + RealT dt_tol = 1.e-8; + EXPECT_NEAR( parameters.dt, expected_dt, dt_tol ); + + tribol::finalize(); +} + +TEST_F( CommonPlaneTest, critical_timestep_impact_limit ) +{ + // This test isolates and validates the Impact Velocity Jump limit. + // By specifying a dynamic interpenetrating velocity, we ensure that the impact + // limit is calculated, and because it is smaller than the energy limit, it governs. + 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; + + // Setup overlapping bounding boxes (interpenetration) + 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. ); + + // Specify a non-zero velocity in the z-direction to trigger the 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.05; // active interpenetrating velocity component + RealT velX2 = 0.; + RealT velY2 = 0.; + RealT velZ2 = 0.05; + + 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; + + 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 values: + // 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 ); + + // The actual gap computed at the common plane centroid by Tribol is: + constexpr RealT g_val = 0.065803604340438485; + // The relative velocity between the two meshes is: + constexpr RealT v_rel = 0.1; + + // Impact limit: dt_impact = 0.025 * dt * (dt * v_rel / g_val) * (stiffness / penalty_stiff_per_area) + constexpr RealT expected_dt = 0.025 * dt * ( ( dt * v_rel ) / g_val ) * ( stiffness2 / penalty_stiff_per_area ); + 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; diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 6bb52366..5e3707ee 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -479,10 +479,14 @@ struct Parameters { 3.0e-1; ///! Max allowable interpenetration as percent of element thickness prior to triggering timestep vote RealT timestep_scale = 1.0; ///! Scale factor (>0) applied to the timestep vote giving users some control over the vote + RealT timestep_chatter_factor = 0.1; ///! Safety factor for contact chatter limit (gamma) + RealT timestep_energy_factor = 0.05; ///! Allowable energy error fraction for energy preservation limit (beta) 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 diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index 895381ab..fb0f9e13 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -187,6 +187,34 @@ void setTimestepScale( IndexT cs_id, RealT scale ) cs->getParameters().timestep_scale = scale; } //------------------------------------------------------------------------------ +void setTimestepChatterFactor( IndexT cs_id, RealT factor ) +{ + if ( factor <= 0. ) { + return; + } + + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + SLIC_ERROR_ROOT_IF( !cs, "tribol::setTimestepChatterFactor(): call tribol::registerCouplingScheme() " + << "prior to calling this routine." ); + + cs->getParameters().timestep_chatter_factor = factor; +} +//------------------------------------------------------------------------------ +void setTimestepEnergyFactor( IndexT cs_id, RealT factor ) +{ + if ( factor <= 0. ) { + return; + } + + auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); + + SLIC_ERROR_ROOT_IF( !cs, "tribol::setTimestepEnergyFactor(): call tribol::registerCouplingScheme() " + << "prior to calling this routine." ); + + cs->getParameters().timestep_energy_factor = factor; +} +//------------------------------------------------------------------------------ void setContactAreaFrac( IndexT cs_id, RealT frac ) { auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); @@ -350,6 +378,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 ) { diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 36e6fa03..e8ab7166 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -142,6 +142,8 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac ); * */ void setTimestepScale( IndexT cs_id, RealT scale ); +void setTimestepChatterFactor( IndexT cs_id, RealT factor ); +void setTimestepEnergyFactor( IndexT cs_id, RealT factor ); /*! * \brief Sets the area fraction for inclusion of a contact overlap @@ -225,6 +227,7 @@ void setBinningProximityScale( IndexT cs_id, RealT binning_proximity_scale ); * */ void enableTimestepVote( IndexT cs_id, const bool enable ); +void enableTimestepStabilityLimits( IndexT cs_id, const bool enable ); /** * @brief Enable Enzyme AD for Jacobian calculations diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 99fca093..8e7ce8a4 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -27,6 +27,7 @@ #include "tribol/common/Parameters.hpp" #include "tribol/physics/Physics.hpp" #include "tribol/physics/ContactFormulationFactory.hpp" +#include "tribol/physics/CommonPlane.hpp" #include "tribol/integ/FE.hpp" namespace tribol { @@ -1349,7 +1350,7 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) // all candidates, not necessarily ones that are deemed to be in contact per // the gap constraint. auto cs_view = getView(); - ArrayT dt_temp_data( { dt, dt }, getAllocatorId() ); + ArrayT dt_temp_data( { dt, dt, dt }, getAllocatorId() ); ArrayViewT dt_temp = dt_temp_data; // [0]: exceed_max_gap1, [1]: exceed_max_gap2, [2]: neg_dt_gap_msg, [3]: neg_dt_vel_proj_msg ArrayT msg_data( { static_cast( false ), static_cast( false ), static_cast( false ), @@ -1434,38 +1435,8 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) v1_dot_n1 = dotProd( vel_f1, fn1, dim ); v2_dot_n2 = dotProd( vel_f2, fn2, dim ); - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "face 1 normal: " << fn1[0] << ", " << fn1[1] << ", " << fn1[2] << std::endl; - // std::cout << "face 2 normal: " << fn2[0] << ", " << fn2[1] << ", " << fn2[2] << std::endl; - // std::cout << " " << std::endl; - // std::cout << "face 1 vel: " << vel_f1[0] << ", " << vel_f1[1] << ", " << vel_f1[2] << std::endl; - // std::cout << "face 2 vel: " << vel_f2[0] << ", " << vel_f2[1] << ", " << vel_f2[2] << std::endl; - // std::cout << " " << std::endl; - // std::cout << "First v1_dot_n1 calc: " << v1_dot_n1 << std::endl; - // std::cout << "First v2_dot_n2 calc: " << v2_dot_n2 << std::endl; - // std::cout << "First v1_dot_n: " << v1_dot_n << std::endl; - // std::cout << "First v2_dot_n: " << v2_dot_n << std::endl; - - // add tiny amount to velocity projections to avoid division by zero. - // Note that if these projections are close to zero, there may be - // stationary interactions or tangential motion. In this case, any - // timestep estimate will be very large, and not control the simulation - RealT tiny = 1.e-12; - RealT tiny1 = ( v1_dot_n >= 0. ) ? tiny : -1. * tiny; - RealT tiny2 = ( v2_dot_n >= 0. ) ? tiny : -1. * tiny; - v1_dot_n += tiny1; - v2_dot_n += tiny2; - // reset tiny velocity based on face normal projections. - tiny1 = ( v1_dot_n1 >= 0. ) ? tiny : -1. * tiny; - tiny2 = ( v2_dot_n2 >= 0. ) ? tiny : -1. * tiny; - v1_dot_n1 += tiny1; - v2_dot_n2 += tiny2; - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "Second v1_dot_n1 calc: " << v1_dot_n1 << std::endl; - // std::cout << "Second v2_dot_n2 calc: " << v2_dot_n2 << std::endl; - // std::cout << "Second v1_dot_n: " << v1_dot_n << std::endl; - // std::cout << "Second v2_dot_n: " << v2_dot_n << std::endl; + // Velocity tolerance for divide-by-zero threshold guarding + constexpr RealT tiny = 1.e-12; // get volume element thicknesses associated with each face in this pair RealT t1 = mesh1.getElementData().m_thickness[index1]; @@ -1502,10 +1473,6 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) // The two cases are: // if v1*n < 0 there is interpen // if v2*n > 0 there is interpen - // - // Note: we compare strictly to 0. here since a 'tiny' value was - // appropriately added to the velocity projections, which is akin - // to some tolerancing effect dt1_vel_check = ( v1_dot_n < 0. ) ? true : false; dt2_vel_check = ( v2_dot_n > 0. ) ? true : false; @@ -1554,14 +1521,14 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) // in excess of the max allowable gap without causing timestep crashes. // // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen - dt1 = ( dt1_check1 ) ? alpha * max_delta1 / v1_dot_n1 : dt1; - dt2 = ( dt2_check1 ) ? alpha * max_delta2 / v2_dot_n2 : dt2; - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "dt1_check1, delta1 and v1_dot_n1: " << dt1_check1 << ", " << max_delta1 << ", " << - // v1_dot_n1 << std::endl; std::cout << "dt2_check1, delta2 and v2_dot_n2: " << dt2_check1 << ", " << - // max_delta2 << ", " << v2_dot_n2 << std::endl; std::cout << "dt1 and dt2: " << dt1 << ", " << dt2 << - // std::endl; + if ( dt1_check1 ) { + // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. + dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : 1.e6; + } + if ( dt2_check1 ) { + // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. + dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : 1.e6; + } // update dt_temp1 only for positive dt1 and/or dt2 if ( dt1 > 0. ) { @@ -1673,15 +1640,14 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) // allows for a soft contact response without a timestep crash. // // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen - dt1 = ( dt1_vel_check ) ? alpha * max_delta1 / v1_dot_n1 : dt1; - dt2 = ( dt2_vel_check ) ? alpha * max_delta2 / v2_dot_n2 : dt2; - - // Keep debug print statements. This routine is still in the testing phase - // std::cout << "dt1_vel_check, (proj_delta_n_1+max_delta1), v1_dot_n1: " << dt1_vel_check << ", " - // << proj_delta_n_1+max_delta1 << ", " << v1_dot_n1 << std::endl; - // std::cout << "dt2_vel_check, (proj_delta_n_2+max_delta2), v2_dot_n2: " << dt2_vel_check << ", " - // << proj_delta_n_2+max_delta2 << ", " << v2_dot_n2 << std::endl; - // std::cout << "dt1 and dt2: " << dt1 << ", " << dt2 << std::endl; + if ( dt1_vel_check ) { + // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. + dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : 1.e6; + } + if ( dt2_vel_check ) { + // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. + dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : 1.e6; + } // update dt_temp2 only for positive dt1 and/or dt2 if ( dt1 > 0. ) { @@ -1707,6 +1673,80 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) } } // end check 2 + + //////////////////////////////////////////////////////////////////////// + // 3. Contact Dynamics and Stability Limits // + // // + // These limits constrain the explicit time integration based on // + // the added stiffness of active contact springs. // + //////////////////////////////////////////////////////////////////////// + + if ( plane.m_inContact && cs_view.enableTimestepStabilityLimits() ) { + RealT stiffness1 = 0.; + RealT stiffness2 = 0.; + auto& enforcement_options = cs_view.getEnforcementOptions(); + const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; + RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; + RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; + + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, + pen_scale1, mesh1.getElementData().m_mat_mod[index1], + mesh1.getElementData().m_penalty_stiffness, stiffness1 ); + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, + pen_scale2, mesh2.getElementData().m_mat_mod[index2], + mesh2.getElementData().m_penalty_stiffness, stiffness2 ); + + if ( stiffness1 > 0. && stiffness2 > 0. ) { + RealT penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); + if ( penalty_stiff_per_area > 0. ) { + // A. Courant (CFL) Stability Limit + RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) ); + RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) ); + RealT dt_CFL = alpha * dt * axom::utilities::min( f_scale1, f_scale2 ); + + // B. Contact Chatter Limit (gamma) + RealT gamma_val = cs_view.getTimestepChatterFactor(); + RealT dt_chatter1 = + 3.14159265358979323846 * gamma_val * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); + RealT dt_chatter2 = + 3.14159265358979323846 * gamma_val * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); + RealT dt_chatter = axom::utilities::min( dt_chatter1, dt_chatter2 ); + + // C. Spurious Energy Generation Limit (beta) + RealT beta_val = cs_view.getTimestepEnergyFactor(); + RealT dt_energy1 = + ( std::sqrt( beta_val ) / 2.0 ) * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); + RealT dt_energy2 = + ( std::sqrt( beta_val ) / 2.0 ) * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); + RealT dt_energy = axom::utilities::min( dt_energy1, dt_energy2 ); + + // D. Impact Velocity Jump Limit (eta = 0.1) + RealT dt_impact1 = 1.e6; + RealT dt_impact2 = 1.e6; + RealT v_rel = std::abs( v1_dot_n1 ) + std::abs( v2_dot_n2 ); + RealT g_val = std::abs( plane.m_gap ) + 1.e-12; + if ( dt1_vel_check ) { + dt_impact1 = 0.025 * dt * ( ( dt * v_rel ) / g_val ) * ( stiffness1 / penalty_stiff_per_area ); + } + if ( dt2_vel_check ) { + dt_impact2 = 0.025 * dt * ( ( dt * v_rel ) / g_val ) * ( stiffness2 / penalty_stiff_per_area ); + } + RealT dt_impact = axom::utilities::min( dt_impact1, dt_impact2 ); + + // Combine all contact dynamics stability limits + RealT dt_crit = axom::utilities::min( + dt_CFL, axom::utilities::min( dt_chatter, axom::utilities::min( dt_energy, dt_impact ) ) ); + + if ( dt_crit > 0. ) { +#ifdef TRIBOL_USE_RAJA + RAJA::atomicMin( &dt_temp[2], axom::utilities::min( dt_crit, 1.e6 ) ); +#else + dt_temp[2] = axom::utilities::min( dt_temp[2], axom::utilities::min( dt_crit, 1.e6 ) ); +#endif + } + } + } + } // end contact dynamics and stability limits } ); // print general messages once @@ -1726,7 +1766,7 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) << "velocity projection calculation." ); ArrayT dt_temp_host( dt_temp_data ); - dt = axom::utilities::min( dt_temp_host[0], dt_temp_host[1] ); + dt = axom::utilities::min( dt_temp_host[0], axom::utilities::min( dt_temp_host[1], dt_temp_host[2] ) ); } //------------------------------------------------------------------------------ diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index df8ec888..5f2c7b5c 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -187,6 +187,30 @@ class CouplingScheme { */ TRIBOL_HOST_DEVICE RealT getTimestepScale() const { return m_parameters.timestep_scale; } + /** + * @brief Return whether contact dynamics stability limits are enabled + * + * @return true if stability limits are enabled, false otherwise + */ + TRIBOL_HOST_DEVICE bool enableTimestepStabilityLimits() const + { + return m_parameters.enable_timestep_stability_limits; + } + + /** + * @brief Get the contact chatter timestep limit factor (gamma) + * + * @return contact chatter factor + */ + TRIBOL_HOST_DEVICE RealT getTimestepChatterFactor() const { return m_parameters.timestep_chatter_factor; } + + /** + * @brief Get the energy preservation timestep limit factor (beta) + * + * @return energy preservation factor + */ + TRIBOL_HOST_DEVICE RealT getTimestepEnergyFactor() const { return m_parameters.timestep_energy_factor; } + /** * @brief Get the gap tolerance that determines in contact face-pairs * diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index 39eae36a..c3b322d4 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -167,37 +167,26 @@ int ApplyNormal( CouplingScheme* cs ) const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; - switch ( pen_enfrc_options.kinematic_calculation ) { - case KINEMATIC_CONSTANT: { - // pre-multiply each spring stiffness by each mesh's penalty scale - auto stiffness1 = pen_scale1 * mesh1.getElementData().m_penalty_stiffness; - auto stiffness2 = pen_scale2 * mesh2.getElementData().m_penalty_stiffness; - // compute the equivalent contact penalty spring stiffness per area - penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - break; - } - case KINEMATIC_ELEMENT: { - // add tiny_length to element thickness to avoid division by zero - auto t1 = mesh1.getElementData().m_thickness[index1] + pen_enfrc_options.tiny_length; - auto t2 = mesh2.getElementData().m_thickness[index2] + pen_enfrc_options.tiny_length; - - if ( t1 < 0. || t2 < 0. ) { - neg_thickness[0] = true; - err[0] = 1; - } - - // compute each element spring stiffness. Pre-multiply the material modulus - // (i.e. material stiffness) by each mesh's penalty scale - auto stiffness1 = pen_scale1 * mesh1.getElementData().m_mat_mod[index1] / t1; - auto stiffness2 = pen_scale2 * mesh2.getElementData().m_mat_mod[index2] / t2; - // compute the equivalent contact penalty spring stiffness per area - penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - break; + RealT stiffness1 = 0.; + RealT stiffness2 = 0.; + RealT t1 = mesh1.getElementData().m_thickness[index1]; + RealT t2 = mesh2.getElementData().m_thickness[index2]; + + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, pen_scale1, + mesh1.getElementData().m_mat_mod[index1], mesh1.getElementData().m_penalty_stiffness, + stiffness1 ); + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, pen_scale2, + mesh2.getElementData().m_mat_mod[index2], mesh2.getElementData().m_penalty_stiffness, + stiffness2 ); + + if ( pen_enfrc_options.kinematic_calculation == KINEMATIC_ELEMENT ) { + if ( t1 <= 0. || t2 <= 0. ) { + neg_thickness[0] = true; + err[0] = 1; } - default: - // no-op, quiet compiler - break; - } // end switch on kinematic penalty calculation option + } + + penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); //////////////////////////////////////////////////// // Compute contact pressure(s) on current overlap // diff --git a/src/tribol/physics/CommonPlane.hpp b/src/tribol/physics/CommonPlane.hpp index a7b0ef8c..11bd3d9b 100644 --- a/src/tribol/physics/CommonPlane.hpp +++ b/src/tribol/physics/CommonPlane.hpp @@ -7,8 +7,44 @@ #define SRC_TRIBOL_PHYSICS_COMMONPLANE_HPP_ #include "Physics.hpp" +#include "tribol/common/Parameters.hpp" namespace tribol { +/*! + * + * \brief computes the individual (uncoupled) kinematic penalty spring stiffness for a single face + * + * \details This routine calculates the stiffness of a single contact spring on one side of + * the contact pair (either mortar or non-mortar), completely independent of (uncoupled from) + * the other side. The resulting uncoupled stiffnesses from both sides are later combined + * in series (harmonic mean) via ComputePenaltyStiffnessPerArea() to obtain the coupled + * equivalent interface stiffness per unit area. + * + * \param [in] kinematic_calc calculation option (constant or element-based) + * \param [in] thickness element thickness + * \param [in] tiny_length numeric safeguard offset to avoid division by zero + * \param [in] pen_scale mesh penalty scale + * \param [in] mat_mod material modulus (bulk modulus) + * \param [in] const_penalty constant penalty value + * \param [out] stiffness computed uncoupled spring stiffness for this face/side + * + */ +TRIBOL_HOST_DEVICE inline void ComputeUncoupledStiffness( const KinematicPenaltyCalculation kinematic_calc, + const RealT thickness, const RealT tiny_length, + const RealT pen_scale, const RealT mat_mod, + const RealT const_penalty, RealT& stiffness ) +{ + if ( kinematic_calc == KINEMATIC_CONSTANT ) { + stiffness = pen_scale * const_penalty; + } else if ( kinematic_calc == KINEMATIC_ELEMENT ) { + RealT denom = thickness; + if ( denom < tiny_length ) { + denom = tiny_length; + } + stiffness = pen_scale * mat_mod / denom; + } +} + /*! * * \brief computes penalty stiffness for Common Plane + Penalty diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 6da885f4..4003ba7c 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1109,8 +1109,11 @@ int TestMesh::tribolSetupAndUpdate( ContactMethod method, EnforcementMethod enfo model, enforcement, BINNING_GRID, ExecutionMode::Sequential ); enableTimestepVote( csIndex, params.enable_timestep_vote ); + enableTimestepStabilityLimits( csIndex, params.enable_timestep_stability_limits ); setTimestepPenFrac( csIndex, params.timestep_pen_frac ); setTimestepScale( csIndex, params.timestep_scale ); + setTimestepChatterFactor( csIndex, params.timestep_chatter_factor ); + setTimestepEnergyFactor( csIndex, params.timestep_energy_factor ); // if enforcement is penalty, register penalty parameters if ( enforcement == PENALTY ) { diff --git a/src/tribol/utils/TestUtils.hpp b/src/tribol/utils/TestUtils.hpp index 52272275..967924c5 100644 --- a/src/tribol/utils/TestUtils.hpp +++ b/src/tribol/utils/TestUtils.hpp @@ -48,8 +48,11 @@ struct TestControlParameters { rate_penalty_ratio( 0.0 ), const_penalty( 1.0 ), enable_timestep_vote( false ), + enable_timestep_stability_limits( false ), timestep_pen_frac( 0.30 ), - timestep_scale( 1.0 ) + timestep_scale( 1.0 ), + timestep_chatter_factor( 0.1 ), + timestep_energy_factor( 0.05 ) { } @@ -69,8 +72,11 @@ struct TestControlParameters { RealT rate_penalty_ratio; RealT const_penalty; bool enable_timestep_vote; + bool enable_timestep_stability_limits; RealT timestep_pen_frac; RealT timestep_scale; + RealT timestep_chatter_factor; + RealT timestep_energy_factor; }; /*! From c5e745a971173f6f660f770aac354f789fa664e7 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 29 May 2026 22:24:35 -0700 Subject: [PATCH 2/5] just keep the cfl limit --- src/tests/tribol_timestep_vote.cpp | 282 ----------------------------- src/tribol/common/Parameters.hpp | 2 - src/tribol/interface/tribol.cpp | 28 --- src/tribol/interface/tribol.hpp | 2 - src/tribol/mesh/CouplingScheme.cpp | 37 +--- src/tribol/mesh/CouplingScheme.hpp | 14 -- src/tribol/utils/TestUtils.cpp | 4 +- src/tribol/utils/TestUtils.hpp | 6 +- 8 files changed, 4 insertions(+), 371 deletions(-) diff --git a/src/tests/tribol_timestep_vote.cpp b/src/tests/tribol_timestep_vote.cpp index e69fd109..0721ee61 100644 --- a/src/tests/tribol_timestep_vote.cpp +++ b/src/tests/tribol_timestep_vote.cpp @@ -1008,288 +1008,6 @@ TEST_F( CommonPlaneTest, critical_timestep_cfl_limit ) tribol::finalize(); } -TEST_F( CommonPlaneTest, critical_timestep_chatter_limit ) -{ - // This test validates the Contact Chatter limit calculation. - // Under the default coefficients, the Energy limit (prefactor 0.1118) is always - // 2.81 times more restrictive than the Chatter limit (prefactor 0.3141). - // We verify that the timestep behaves exactly as predicted by the companion limits. - 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. ); - - 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 = 1.0; - parameters.timestep_chatter_factor = 0.01; // gamma = 0.01 - - 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 Chatter limit: - // 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 ); - - // Chatter limit: dt_chatter = pi * gamma * dt * sqrt( stiffness / penalty_stiff_per_area ) - constexpr RealT pi = 3.14159265358979323846; - const RealT dt_chatter1 = - pi * parameters.timestep_chatter_factor * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); - const RealT dt_chatter2 = - pi * parameters.timestep_chatter_factor * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); - const RealT expected_dt = std::min( dt_chatter1, dt_chatter2 ); - RealT dt_tol = 1.e-6; - EXPECT_NEAR( parameters.dt, expected_dt, dt_tol ); - - tribol::finalize(); -} - -TEST_F( CommonPlaneTest, critical_timestep_energy_limit ) -{ - // This test isolates and validates the Spurious Energy Generation limit. - // Under stationary contact (zero velocity), the impact limit is bypassed, - // and the energy limit naturally 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 = 1.0; - - 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 Energy limit 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 ); - - // Energy limit: dt_energy = (sqrt(beta) / 2) * dt * sqrt( stiffness / penalty_stiff_per_area ) - const RealT dt_energy1 = - ( std::sqrt( parameters.timestep_energy_factor ) / 2.0 ) * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); - const RealT dt_energy2 = - ( std::sqrt( parameters.timestep_energy_factor ) / 2.0 ) * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); - const RealT expected_dt = std::min( dt_energy1, dt_energy2 ); - RealT dt_tol = 1.e-8; - EXPECT_NEAR( parameters.dt, expected_dt, dt_tol ); - - tribol::finalize(); -} - -TEST_F( CommonPlaneTest, critical_timestep_impact_limit ) -{ - // This test isolates and validates the Impact Velocity Jump limit. - // By specifying a dynamic interpenetrating velocity, we ensure that the impact - // limit is calculated, and because it is smaller than the energy limit, it governs. - 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; - - // Setup overlapping bounding boxes (interpenetration) - 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. ); - - // Specify a non-zero velocity in the z-direction to trigger the 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.05; // active interpenetrating velocity component - RealT velX2 = 0.; - RealT velY2 = 0.; - RealT velZ2 = 0.05; - - 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; - - 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 values: - // 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 ); - - // The actual gap computed at the common plane centroid by Tribol is: - constexpr RealT g_val = 0.065803604340438485; - // The relative velocity between the two meshes is: - constexpr RealT v_rel = 0.1; - - // Impact limit: dt_impact = 0.025 * dt * (dt * v_rel / g_val) * (stiffness / penalty_stiff_per_area) - constexpr RealT expected_dt = 0.025 * dt * ( ( dt * v_rel ) / g_val ) * ( stiffness2 / penalty_stiff_per_area ); - 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; diff --git a/src/tribol/common/Parameters.hpp b/src/tribol/common/Parameters.hpp index 5e3707ee..aa661f31 100644 --- a/src/tribol/common/Parameters.hpp +++ b/src/tribol/common/Parameters.hpp @@ -479,8 +479,6 @@ struct Parameters { 3.0e-1; ///! Max allowable interpenetration as percent of element thickness prior to triggering timestep vote RealT timestep_scale = 1.0; ///! Scale factor (>0) applied to the timestep vote giving users some control over the vote - RealT timestep_chatter_factor = 0.1; ///! Safety factor for contact chatter limit (gamma) - RealT timestep_energy_factor = 0.05; ///! Allowable energy error fraction for energy preservation limit (beta) int vis_cycle_incr = 100; ///! Frequency for visualizations dumps VisType vis_type = VIS_NONE; ///! Type of interface physics visualization output diff --git a/src/tribol/interface/tribol.cpp b/src/tribol/interface/tribol.cpp index fb0f9e13..46f39336 100644 --- a/src/tribol/interface/tribol.cpp +++ b/src/tribol/interface/tribol.cpp @@ -187,34 +187,6 @@ void setTimestepScale( IndexT cs_id, RealT scale ) cs->getParameters().timestep_scale = scale; } //------------------------------------------------------------------------------ -void setTimestepChatterFactor( IndexT cs_id, RealT factor ) -{ - if ( factor <= 0. ) { - return; - } - - auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); - - SLIC_ERROR_ROOT_IF( !cs, "tribol::setTimestepChatterFactor(): call tribol::registerCouplingScheme() " - << "prior to calling this routine." ); - - cs->getParameters().timestep_chatter_factor = factor; -} -//------------------------------------------------------------------------------ -void setTimestepEnergyFactor( IndexT cs_id, RealT factor ) -{ - if ( factor <= 0. ) { - return; - } - - auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); - - SLIC_ERROR_ROOT_IF( !cs, "tribol::setTimestepEnergyFactor(): call tribol::registerCouplingScheme() " - << "prior to calling this routine." ); - - cs->getParameters().timestep_energy_factor = factor; -} -//------------------------------------------------------------------------------ void setContactAreaFrac( IndexT cs_id, RealT frac ) { auto cs = CouplingSchemeManager::getInstance().findData( cs_id ); diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index e8ab7166..98d75af4 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -142,8 +142,6 @@ void setTimestepPenFrac( IndexT cs_id, RealT frac ); * */ void setTimestepScale( IndexT cs_id, RealT scale ); -void setTimestepChatterFactor( IndexT cs_id, RealT factor ); -void setTimestepEnergyFactor( IndexT cs_id, RealT factor ); /*! * \brief Sets the area fraction for inclusion of a contact overlap diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 8e7ce8a4..64008d7a 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -1699,43 +1699,10 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) if ( stiffness1 > 0. && stiffness2 > 0. ) { RealT penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); if ( penalty_stiff_per_area > 0. ) { - // A. Courant (CFL) Stability Limit + // Rescale the host timestep by the added interface stiffness. RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) ); RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) ); - RealT dt_CFL = alpha * dt * axom::utilities::min( f_scale1, f_scale2 ); - - // B. Contact Chatter Limit (gamma) - RealT gamma_val = cs_view.getTimestepChatterFactor(); - RealT dt_chatter1 = - 3.14159265358979323846 * gamma_val * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); - RealT dt_chatter2 = - 3.14159265358979323846 * gamma_val * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); - RealT dt_chatter = axom::utilities::min( dt_chatter1, dt_chatter2 ); - - // C. Spurious Energy Generation Limit (beta) - RealT beta_val = cs_view.getTimestepEnergyFactor(); - RealT dt_energy1 = - ( std::sqrt( beta_val ) / 2.0 ) * dt * std::sqrt( stiffness1 / penalty_stiff_per_area ); - RealT dt_energy2 = - ( std::sqrt( beta_val ) / 2.0 ) * dt * std::sqrt( stiffness2 / penalty_stiff_per_area ); - RealT dt_energy = axom::utilities::min( dt_energy1, dt_energy2 ); - - // D. Impact Velocity Jump Limit (eta = 0.1) - RealT dt_impact1 = 1.e6; - RealT dt_impact2 = 1.e6; - RealT v_rel = std::abs( v1_dot_n1 ) + std::abs( v2_dot_n2 ); - RealT g_val = std::abs( plane.m_gap ) + 1.e-12; - if ( dt1_vel_check ) { - dt_impact1 = 0.025 * dt * ( ( dt * v_rel ) / g_val ) * ( stiffness1 / penalty_stiff_per_area ); - } - if ( dt2_vel_check ) { - dt_impact2 = 0.025 * dt * ( ( dt * v_rel ) / g_val ) * ( stiffness2 / penalty_stiff_per_area ); - } - RealT dt_impact = axom::utilities::min( dt_impact1, dt_impact2 ); - - // Combine all contact dynamics stability limits - RealT dt_crit = axom::utilities::min( - dt_CFL, axom::utilities::min( dt_chatter, axom::utilities::min( dt_energy, dt_impact ) ) ); + RealT dt_crit = alpha * dt * axom::utilities::min( f_scale1, f_scale2 ); if ( dt_crit > 0. ) { #ifdef TRIBOL_USE_RAJA diff --git a/src/tribol/mesh/CouplingScheme.hpp b/src/tribol/mesh/CouplingScheme.hpp index 5f2c7b5c..bdbc8b11 100644 --- a/src/tribol/mesh/CouplingScheme.hpp +++ b/src/tribol/mesh/CouplingScheme.hpp @@ -197,20 +197,6 @@ class CouplingScheme { return m_parameters.enable_timestep_stability_limits; } - /** - * @brief Get the contact chatter timestep limit factor (gamma) - * - * @return contact chatter factor - */ - TRIBOL_HOST_DEVICE RealT getTimestepChatterFactor() const { return m_parameters.timestep_chatter_factor; } - - /** - * @brief Get the energy preservation timestep limit factor (beta) - * - * @return energy preservation factor - */ - TRIBOL_HOST_DEVICE RealT getTimestepEnergyFactor() const { return m_parameters.timestep_energy_factor; } - /** * @brief Get the gap tolerance that determines in contact face-pairs * diff --git a/src/tribol/utils/TestUtils.cpp b/src/tribol/utils/TestUtils.cpp index 4003ba7c..9bda6fee 100644 --- a/src/tribol/utils/TestUtils.cpp +++ b/src/tribol/utils/TestUtils.cpp @@ -1112,8 +1112,6 @@ int TestMesh::tribolSetupAndUpdate( ContactMethod method, EnforcementMethod enfo enableTimestepStabilityLimits( csIndex, params.enable_timestep_stability_limits ); setTimestepPenFrac( csIndex, params.timestep_pen_frac ); setTimestepScale( csIndex, params.timestep_scale ); - setTimestepChatterFactor( csIndex, params.timestep_chatter_factor ); - setTimestepEnergyFactor( csIndex, params.timestep_energy_factor ); // if enforcement is penalty, register penalty parameters if ( enforcement == PENALTY ) { @@ -1920,4 +1918,4 @@ mfem::Vector ExplicitMechanics::ComputeInvMass( mfem::ParFiniteElementSpace& fes #endif -} // namespace mfem_ext \ No newline at end of file +} // namespace mfem_ext diff --git a/src/tribol/utils/TestUtils.hpp b/src/tribol/utils/TestUtils.hpp index 967924c5..d0fc109d 100644 --- a/src/tribol/utils/TestUtils.hpp +++ b/src/tribol/utils/TestUtils.hpp @@ -50,9 +50,7 @@ struct TestControlParameters { enable_timestep_vote( false ), enable_timestep_stability_limits( false ), timestep_pen_frac( 0.30 ), - timestep_scale( 1.0 ), - timestep_chatter_factor( 0.1 ), - timestep_energy_factor( 0.05 ) + timestep_scale( 1.0 ) { } @@ -75,8 +73,6 @@ struct TestControlParameters { bool enable_timestep_stability_limits; RealT timestep_pen_frac; RealT timestep_scale; - RealT timestep_chatter_factor; - RealT timestep_energy_factor; }; /*! From 3b5d900dca69cde84b228f95d5c0b87ccbe79bc2 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Fri, 29 May 2026 22:55:36 -0700 Subject: [PATCH 3/5] clear up magic numbers; fix indexing bug --- src/tests/tribol_common_plane_penalty.cpp | 90 +++ src/tribol/mesh/CouplingScheme.cpp | 664 +++++++++++----------- src/tribol/physics/CommonPlane.cpp | 23 +- 3 files changed, 439 insertions(+), 338 deletions(-) diff --git a/src/tests/tribol_common_plane_penalty.cpp b/src/tests/tribol_common_plane_penalty.cpp index 5f264725..d847d0df 100644 --- a/src/tests/tribol_common_plane_penalty.cpp +++ b/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( cs )->getMesh1().getView(); + const auto mesh2 = const_cast( 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; diff --git a/src/tribol/mesh/CouplingScheme.cpp b/src/tribol/mesh/CouplingScheme.cpp index 64008d7a..235729b0 100644 --- a/src/tribol/mesh/CouplingScheme.cpp +++ b/src/tribol/mesh/CouplingScheme.cpp @@ -1357,364 +1357,370 @@ void CouplingScheme::computeCommonPlaneTimeStep( RealT& dt ) static_cast( false ) }, getAllocatorId() ); ArrayViewT msg = msg_data; - forAllExec( getExecutionMode(), getNumActivePairs(), - [cs_view, dim, proj_ratio, msg, dt_temp, dt] TRIBOL_HOST_DEVICE( IndexT i ) { - auto& cg_view = cs_view.getCompGeomView(); - auto& plane = cg_view.getCommonPlane( i ); - - auto& mesh1 = cs_view.getMesh1View(); - auto& mesh2 = cs_view.getMesh2View(); - - // get pair indices - IndexT index1 = plane.getCpElementId1(); - IndexT index2 = plane.getCpElementId2(); - - constexpr int max_dim = 3; - constexpr int max_nodes_per_elem = 4; - StackArrayT x1; - StackArrayT v1; - mesh1.getFaceCoords( index1, x1 ); - mesh1.getFaceVelocities( index1, v1 ); - - StackArrayT x2; - StackArrayT v2; - mesh2.getFaceCoords( index2, x2 ); - mesh2.getFaceVelocities( index2, v2 ); - - ///////////////////////////////////////////////////////////// - // calculate face velocities at projected overlap centroid // - ///////////////////////////////////////////////////////////// - StackArrayT vel_f1; - StackArrayT vel_f2; - initRealArray( vel_f1, dim, 0.0 ); - initRealArray( vel_f2, dim, 0.0 ); - - // interpolate nodal velocity at overlap centroid as projected - // onto face 1 - RealT cXf1 = plane.m_cXf1; - RealT cYf1 = plane.m_cYf1; - RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; - GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); - // interpolate nodal velocity at overlap centroid as projected - // onto face 2 - RealT cXf2 = plane.m_cXf2; - RealT cYf2 = plane.m_cYf2; - RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; - GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); - - //////////////////////////////////////////////// - // // - // Compute Timestep Vote Based on a Few Cases // - // // - //////////////////////////////////////////////// - - /////////////////////////////////////////////// - // compute data common to all timestep votes // - /////////////////////////////////////////////// - - // compute velocity projections: - // compute the dot product between the face velocities - // at the overlap-centroid-to-face projected centroid and each - // face's outward unit normal AND the overlap normal. - RealT v1_dot_n, v2_dot_n, v1_dot_n1, v2_dot_n2; - RealT overlapNormal[max_dim]; - overlapNormal[0] = plane.m_nX; - overlapNormal[1] = plane.m_nY; - if ( dim == 3 ) { - overlapNormal[2] = plane.m_nZ; - } + forAllExec( + getExecutionMode(), getNumActivePairs(), + [cs_view, dim, proj_ratio, msg, dt_temp, dt] TRIBOL_HOST_DEVICE( IndexT i ) { + auto& cg_view = cs_view.getCompGeomView(); + auto& plane = cg_view.getCommonPlane( i ); + + auto& mesh1 = cs_view.getMesh1View(); + auto& mesh2 = cs_view.getMesh2View(); + + // get pair indices + IndexT index1 = plane.getCpElementId1(); + IndexT index2 = plane.getCpElementId2(); + + constexpr int max_dim = 3; + constexpr int max_nodes_per_elem = 4; + StackArrayT x1; + StackArrayT v1; + mesh1.getFaceCoords( index1, x1 ); + mesh1.getFaceVelocities( index1, v1 ); + + StackArrayT x2; + StackArrayT v2; + mesh2.getFaceCoords( index2, x2 ); + mesh2.getFaceVelocities( index2, v2 ); + + ///////////////////////////////////////////////////////////// + // calculate face velocities at projected overlap centroid // + ///////////////////////////////////////////////////////////// + StackArrayT vel_f1; + StackArrayT vel_f2; + initRealArray( vel_f1, dim, 0.0 ); + initRealArray( vel_f2, dim, 0.0 ); + + // interpolate nodal velocity at overlap centroid as projected + // onto face 1 + RealT cXf1 = plane.m_cXf1; + RealT cYf1 = plane.m_cYf1; + RealT cZf1 = ( dim == 3 ) ? plane.m_cZf1 : 0.; + GalerkinEval( x1, cXf1, cYf1, cZf1, LINEAR, PHYSICAL, dim, dim, v1, vel_f1 ); + // interpolate nodal velocity at overlap centroid as projected + // onto face 2 + RealT cXf2 = plane.m_cXf2; + RealT cYf2 = plane.m_cYf2; + RealT cZf2 = ( dim == 3 ) ? plane.m_cZf2 : 0.; + GalerkinEval( x2, cXf2, cYf2, cZf2, LINEAR, PHYSICAL, dim, dim, v2, vel_f2 ); + + //////////////////////////////////////////////// + // // + // Compute Timestep Vote Based on a Few Cases // + // // + //////////////////////////////////////////////// + + /////////////////////////////////////////////// + // compute data common to all timestep votes // + /////////////////////////////////////////////// + + // compute velocity projections: + // compute the dot product between the face velocities + // at the overlap-centroid-to-face projected centroid and each + // face's outward unit normal AND the overlap normal. + RealT v1_dot_n, v2_dot_n, v1_dot_n1, v2_dot_n2; + RealT overlapNormal[max_dim]; + overlapNormal[0] = plane.m_nX; + overlapNormal[1] = plane.m_nY; + if ( dim == 3 ) { + overlapNormal[2] = plane.m_nZ; + } - // get face normals - RealT fn1[max_dim], fn2[max_dim]; - mesh1.getFaceNormal( index1, fn1 ); - mesh2.getFaceNormal( index2, fn2 ); - - // compute projections - v1_dot_n = dotProd( vel_f1, overlapNormal, dim ); - v2_dot_n = dotProd( vel_f2, overlapNormal, dim ); - v1_dot_n1 = dotProd( vel_f1, fn1, dim ); - v2_dot_n2 = dotProd( vel_f2, fn2, dim ); - - // Velocity tolerance for divide-by-zero threshold guarding - constexpr RealT tiny = 1.e-12; - - // get volume element thicknesses associated with each face in this pair - RealT t1 = mesh1.getElementData().m_thickness[index1]; - RealT t2 = mesh2.getElementData().m_thickness[index2]; - - // compute the existing gap vector (recall gap is x1-x2 by convention) - RealT gapVec[max_dim]; - gapVec[0] = plane.m_cXf1 - plane.m_cXf2; - gapVec[1] = plane.m_cYf1 - plane.m_cYf2; - if ( dim == 3 ) { - gapVec[2] = plane.m_cZf1 - plane.m_cZf2; - } + // get face normals + RealT fn1[max_dim], fn2[max_dim]; + mesh1.getFaceNormal( index1, fn1 ); + mesh2.getFaceNormal( index2, fn2 ); + + // compute projections + v1_dot_n = dotProd( vel_f1, overlapNormal, dim ); + v2_dot_n = dotProd( vel_f2, overlapNormal, dim ); + v1_dot_n1 = dotProd( vel_f1, fn1, dim ); + v2_dot_n2 = dotProd( vel_f2, fn2, dim ); + + // Velocity tolerance for divide-by-zero threshold guarding + constexpr RealT tiny = 1.e-12; + // Large fallback value used when the time-to-threshold is effectively infinite. + constexpr RealT fallback_dt = 1.e6; + + // get volume element thicknesses associated with each face in this pair + RealT t1 = mesh1.getElementData().m_thickness[index1]; + RealT t2 = mesh2.getElementData().m_thickness[index2]; + + // compute the existing gap vector (recall gap is x1-x2 by convention) + RealT gapVec[max_dim]; + gapVec[0] = plane.m_cXf1 - plane.m_cXf2; + gapVec[1] = plane.m_cYf1 - plane.m_cYf2; + if ( dim == 3 ) { + gapVec[2] = plane.m_cZf1 - plane.m_cZf2; + } - // compute the dot product between gap vector and the outward unit face normals. - RealT gap_f1_n1 = dotProd( gapVec, fn1, dim ); - RealT gap_f2_n2 = dotProd( gapVec, fn2, dim ); - - RealT dt1 = 1.e6; // initialize as large number - RealT dt2 = 1.e6; // initialize as large number - RealT alpha = cs_view.getTimestepScale(); // multiplier on timestep estimate - bool dt1_check1 = false; - bool dt2_check1 = false; - bool dt1_vel_check = false; - bool dt2_vel_check = false; - - // maximum allowable interpenetration in the normal direction of each element - RealT max_delta1 = proj_ratio * t1; - RealT max_delta2 = proj_ratio * t2; - - // Separation or interpenetration trigger for check 1 and 2: - // check if there is further interpen or separation based on the - // velocity projection in the direction of the common-plane normal, - // which is in the direction of face-2 normal. - // The two cases are: - // if v1*n < 0 there is interpen - // if v2*n > 0 there is interpen - dt1_vel_check = ( v1_dot_n < 0. ) ? true : false; - dt2_vel_check = ( v2_dot_n > 0. ) ? true : false; - - ////////////////////////////////////////////////////////////////////////// - // Check 1. Current interpenetration gap exceeds max allowable interpen // - ////////////////////////////////////////////////////////////////////////// - - // check if face-pair is in contact (i.e. gap < gap_tol), which is determined - // in Common Plane ApplyNormal<>() routine - if ( plane.m_inContact ) { - // compute the difference between the 'face-gaps' and the max allowable - // interpen as a function of element thickness. Note, we have to use the - // gap projected onto the outward unit face-normal to check against the - // max allowable gap as a factor of the thickness in the element normal - // direction - RealT delta1 = max_delta1 - gap_f1_n1; // >0 not exceeding max allowable - RealT delta2 = max_delta2 + gap_f2_n2; // >0 not exceeding max allowable - - auto exceed_max_gap1 = ( delta1 < 0. ) ? true : false; - auto exceed_max_gap2 = ( delta2 < 0. ) ? true : false; - - // if velocity projection indicates further interpenetration, and the gaps - // EXCEED max allowable, then compute time step estimates to reduce overlap - dt1_check1 = ( dt1_vel_check ) ? exceed_max_gap1 : false; - dt2_check1 = ( dt2_vel_check ) ? exceed_max_gap2 : false; + // compute the dot product between gap vector and the outward unit face normals. + RealT gap_f1_n1 = dotProd( gapVec, fn1, dim ); + RealT gap_f2_n2 = dotProd( gapVec, fn2, dim ); + + RealT dt1 = fallback_dt; + RealT dt2 = fallback_dt; + RealT alpha = cs_view.getTimestepScale(); // multiplier on timestep estimate + bool dt1_check1 = false; + bool dt2_check1 = false; + bool dt1_vel_check = false; + bool dt2_vel_check = false; + + // maximum allowable interpenetration in the normal direction of each element + RealT max_delta1 = proj_ratio * t1; + RealT max_delta2 = proj_ratio * t2; + + // Separation or interpenetration trigger for check 1 and 2: + // check if there is further interpen or separation based on the + // velocity projection in the direction of the common-plane normal, + // which is in the direction of face-2 normal. + // The two cases are: + // if v1*n < 0 there is interpen + // if v2*n > 0 there is interpen + dt1_vel_check = ( v1_dot_n < 0. ) ? true : false; + dt2_vel_check = ( v2_dot_n > 0. ) ? true : false; + + ////////////////////////////////////////////////////////////////////////// + // Check 1. Current interpenetration gap exceeds max allowable interpen // + ////////////////////////////////////////////////////////////////////////// + + // check if face-pair is in contact (i.e. gap < gap_tol), which is determined + // in Common Plane ApplyNormal<>() routine + if ( plane.m_inContact ) { + // compute the difference between the 'face-gaps' and the max allowable + // interpen as a function of element thickness. Note, we have to use the + // gap projected onto the outward unit face-normal to check against the + // max allowable gap as a factor of the thickness in the element normal + // direction + RealT delta1 = max_delta1 - gap_f1_n1; // >0 not exceeding max allowable + RealT delta2 = max_delta2 + gap_f2_n2; // >0 not exceeding max allowable + + auto exceed_max_gap1 = ( delta1 < 0. ) ? true : false; + auto exceed_max_gap2 = ( delta2 < 0. ) ? true : false; + + // if velocity projection indicates further interpenetration, and the gaps + // EXCEED max allowable, then compute time step estimates to reduce overlap + dt1_check1 = ( dt1_vel_check ) ? exceed_max_gap1 : false; + dt2_check1 = ( dt2_vel_check ) ? exceed_max_gap2 : false; #ifdef TRIBOL_USE_RAJA - RAJA::atomicMax( &msg[0], static_cast( exceed_max_gap1 ) ); - RAJA::atomicMax( &msg[1], static_cast( exceed_max_gap2 ) ); + RAJA::atomicMax( &msg[0], static_cast( exceed_max_gap1 ) ); + RAJA::atomicMax( &msg[1], static_cast( exceed_max_gap2 ) ); #else - msg[0] = exceed_max_gap1; - msg[1] = exceed_max_gap2; + msg[0] = exceed_max_gap1; + msg[1] = exceed_max_gap2; #endif - // compute dt for face 1 and 2 based on the velocity and gap projections onto - // the face-normals for faces where currect gap exceeds max allowable gap. - // - // NOTE: - // - // This calculation RESETS the current gap to be g = 0, and computes a timestep - // such that the velocity projection of the overlap-to-face projected overlap - // centroid does not exceed the max allowable gap. - // - // This avoid a timestep crash in the case that the current gap barely exceeds - // the max allowable and also allows a soft contact response with interpen - // in excess of the max allowable gap without causing timestep crashes. - // - // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen - if ( dt1_check1 ) { - // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. - dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : 1.e6; - } - if ( dt2_check1 ) { - // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. - dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : 1.e6; - } - - // update dt_temp1 only for positive dt1 and/or dt2 - if ( dt1 > 0. ) { + // compute dt for face 1 and 2 based on the velocity and gap projections onto + // the face-normals for faces where currect gap exceeds max allowable gap. + // + // NOTE: + // + // This calculation RESETS the current gap to be g = 0, and computes a timestep + // such that the velocity projection of the overlap-to-face projected overlap + // centroid does not exceed the max allowable gap. + // + // This avoid a timestep crash in the case that the current gap barely exceeds + // the max allowable and also allows a soft contact response with interpen + // in excess of the max allowable gap without causing timestep crashes. + // + // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen + if ( dt1_check1 ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : fallback_dt; + } + if ( dt2_check1 ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : fallback_dt; + } + + // update dt_temp1 only for positive dt1 and/or dt2 + if ( dt1 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[0], axom::utilities::min( dt1, 1.e6 ) ); + RAJA::atomicMin( &dt_temp[0], axom::utilities::min( dt1, fallback_dt ) ); #else - dt_temp[0] = axom::utilities::min(dt_temp[0], axom::utilities::min(dt1, 1.e6)); + dt_temp[0] = axom::utilities::min( dt_temp[0], axom::utilities::min( dt1, fallback_dt ) ); #endif - } - if ( dt2 > 0. ) { + } + if ( dt2 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[0], axom::utilities::min( 1.e6, dt2 ) ); + RAJA::atomicMin( &dt_temp[0], axom::utilities::min( fallback_dt, dt2 ) ); #else - dt_temp[0] = axom::utilities::min(dt_temp[0], axom::utilities::min(1.e6, dt2)); + dt_temp[0] = axom::utilities::min( dt_temp[0], axom::utilities::min( fallback_dt, dt2 ) ); #endif - } + } - if ( dt1 < 0. || dt2 < 0. ) { + if ( dt1 < 0. || dt2 < 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMax( &msg[2], static_cast( true ) ); + RAJA::atomicMax( &msg[2], static_cast( true ) ); #else - msg[2] = true; + msg[2] = true; #endif - } - } // end case 1 - - //////////////////////////////////////////////////////////////////////// - // 2. Velocity projection exceeds max interpenetration // - // // - // Note: This is performed for all contact candidates even if they // - // are not 'in contact' per the common-plane method. Every // - // contact candidate has a contact plane // - //////////////////////////////////////////////////////////////////////// - - { - // compute the delta between the velocity projection of each face-projected - // common plane centroid location - // - // First project each face-projected common plane centroid using linear velocity - // projection as approximation of configuration next cycle - RealT proj_delta_x1 = plane.m_cXf1 + dt * vel_f1[0]; - RealT proj_delta_y1 = plane.m_cYf1 + dt * vel_f1[1]; - RealT proj_delta_z1 = 0.; - - RealT proj_delta_x2 = plane.m_cXf2 + dt * vel_f2[0]; - RealT proj_delta_y2 = plane.m_cYf2 + dt * vel_f2[1]; - RealT proj_delta_z2 = 0.; - - // Second compute the amount of interpenetration as the difference between the two - // velocity projected points - RealT proj_delta_x1_fixed = proj_delta_x1; - RealT proj_delta_y1_fixed = proj_delta_y1; - - proj_delta_x1 -= proj_delta_x2; - proj_delta_y1 -= proj_delta_y2; - - proj_delta_x2 -= proj_delta_x1_fixed; - proj_delta_y2 -= proj_delta_y1_fixed; - - // compute the dot product between each face's delta and the OTHER - // face's outward unit normal. This is the magnitude of interpenetration - // of one face's projected overlap-centroid in the 'thickness-direction' - // of the other face (with whom in may be in contact currently, or in - // a velocity projected sense). - RealT proj_delta_n_1 = proj_delta_x1 * fn2[0] + proj_delta_y1 * fn2[1]; - RealT proj_delta_n_2 = proj_delta_x2 * fn1[0] + proj_delta_y2 * fn1[1]; - - if ( dim == 3 ) { - // project the z-component - proj_delta_z1 = plane.m_cZf1 + dt * vel_f1[2]; - proj_delta_z2 = plane.m_cZf2 + dt * vel_f2[2]; - - RealT proj_delta_z1_fixed = proj_delta_z1; - - // compute difference between each projected point - proj_delta_z1 -= proj_delta_z2; - proj_delta_z2 -= proj_delta_z1_fixed; - - // add the z-component of the projection onto the face normal - proj_delta_n_1 += proj_delta_z1 * fn2[2]; - proj_delta_n_2 += proj_delta_z2 * fn1[2]; - } - - // Reset the dt velocity check only for faces with continued interpen that exceeds the - // max allowable gap AND where the current gap did NOT exceed that face's max allowable - // gap per check 1 (would result in same dt calc). - // - // Note: - // If proj_delta_n_i < 0, (i=1,2) there is interpen from the velocity projection. - // Check this interpen against the maximum allowable to determine if a velocity projection - // timestep estimate is still required. - if ( dt1_vel_check && !dt1_check1 ) // continued interpen - { - dt1_vel_check = ( proj_delta_n_1 < 0. ) - ? ( ( std::abs( proj_delta_n_1 ) > max_delta1 ) ? true : false ) - : false; - } - - if ( dt2_vel_check && !dt2_check1 ) // continued interpen - { - dt2_vel_check = ( proj_delta_n_2 < 0. ) - ? ( ( std::abs( proj_delta_n_2 ) > max_delta2 ) ? true : false ) - : false; - } - - // compute velocity projection based dt (check 2) using a RESET gap (g=0) such that - // the velocity projected gap does not exceed the max allowable gap. This avoid timestep - // crashes for velocity projected gaps slightly in excess of the max allowable and still - // allows for a soft contact response without a timestep crash. - // - // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen - if ( dt1_vel_check ) { - // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. - dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : 1.e6; - } - if ( dt2_vel_check ) { - // If velocity is effectively zero, the time to exceed gap is infinite; default to 1.e6. - dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : 1.e6; - } - - // update dt_temp2 only for positive dt1 and/or dt2 - if ( dt1 > 0. ) { + } + } // end case 1 + + //////////////////////////////////////////////////////////////////////// + // 2. Velocity projection exceeds max interpenetration // + // // + // Note: This is performed for all contact candidates even if they // + // are not 'in contact' per the common-plane method. Every // + // contact candidate has a contact plane // + //////////////////////////////////////////////////////////////////////// + + { + // compute the delta between the velocity projection of each face-projected + // common plane centroid location + // + // First project each face-projected common plane centroid using linear velocity + // projection as approximation of configuration next cycle + RealT proj_delta_x1 = plane.m_cXf1 + dt * vel_f1[0]; + RealT proj_delta_y1 = plane.m_cYf1 + dt * vel_f1[1]; + RealT proj_delta_z1 = 0.; + + RealT proj_delta_x2 = plane.m_cXf2 + dt * vel_f2[0]; + RealT proj_delta_y2 = plane.m_cYf2 + dt * vel_f2[1]; + RealT proj_delta_z2 = 0.; + + // Second compute the amount of interpenetration as the difference between the two + // velocity projected points + RealT proj_delta_x1_fixed = proj_delta_x1; + RealT proj_delta_y1_fixed = proj_delta_y1; + + proj_delta_x1 -= proj_delta_x2; + proj_delta_y1 -= proj_delta_y2; + + proj_delta_x2 -= proj_delta_x1_fixed; + proj_delta_y2 -= proj_delta_y1_fixed; + + // compute the dot product between each face's delta and the OTHER + // face's outward unit normal. This is the magnitude of interpenetration + // of one face's projected overlap-centroid in the 'thickness-direction' + // of the other face (with whom in may be in contact currently, or in + // a velocity projected sense). + RealT proj_delta_n_1 = proj_delta_x1 * fn2[0] + proj_delta_y1 * fn2[1]; + RealT proj_delta_n_2 = proj_delta_x2 * fn1[0] + proj_delta_y2 * fn1[1]; + + if ( dim == 3 ) { + // project the z-component + proj_delta_z1 = plane.m_cZf1 + dt * vel_f1[2]; + proj_delta_z2 = plane.m_cZf2 + dt * vel_f2[2]; + + RealT proj_delta_z1_fixed = proj_delta_z1; + + // compute difference between each projected point + proj_delta_z1 -= proj_delta_z2; + proj_delta_z2 -= proj_delta_z1_fixed; + + // add the z-component of the projection onto the face normal + proj_delta_n_1 += proj_delta_z1 * fn2[2]; + proj_delta_n_2 += proj_delta_z2 * fn1[2]; + } + + // Reset the dt velocity check only for faces with continued interpen that exceeds the + // max allowable gap AND where the current gap did NOT exceed that face's max allowable + // gap per check 1 (would result in same dt calc). + // + // Note: + // If proj_delta_n_i < 0, (i=1,2) there is interpen from the velocity projection. + // Check this interpen against the maximum allowable to determine if a velocity projection + // timestep estimate is still required. + if ( dt1_vel_check && !dt1_check1 ) // continued interpen + { + dt1_vel_check = + ( proj_delta_n_1 < 0. ) ? ( ( std::abs( proj_delta_n_1 ) > max_delta1 ) ? true : false ) : false; + } + + if ( dt2_vel_check && !dt2_check1 ) // continued interpen + { + dt2_vel_check = + ( proj_delta_n_2 < 0. ) ? ( ( std::abs( proj_delta_n_2 ) > max_delta2 ) ? true : false ) : false; + } + + // compute velocity projection based dt (check 2) using a RESET gap (g=0) such that + // the velocity projected gap does not exceed the max allowable gap. This avoid timestep + // crashes for velocity projected gaps slightly in excess of the max allowable and still + // allows for a soft contact response without a timestep crash. + // + // v1_dot_n1 > 0 and v2_dot_n2 > 0 for further interpen + if ( dt1_vel_check ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt1 = ( std::abs( v1_dot_n1 ) > tiny ) ? ( alpha * max_delta1 / v1_dot_n1 ) : fallback_dt; + } + if ( dt2_vel_check ) { + // If velocity is effectively zero, the time to exceed gap is effectively infinite. + dt2 = ( std::abs( v2_dot_n2 ) > tiny ) ? ( alpha * max_delta2 / v2_dot_n2 ) : fallback_dt; + } + + // update dt_temp2 only for positive dt1 and/or dt2 + if ( dt1 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[1], axom::utilities::min( dt1, 1.e6 ) ); + RAJA::atomicMin( &dt_temp[1], axom::utilities::min( dt1, fallback_dt ) ); #else - dt_temp[1] = axom::utilities::min(dt_temp[1], axom::utilities::min(dt1, 1.e6)); + dt_temp[1] = axom::utilities::min( dt_temp[1], axom::utilities::min( dt1, fallback_dt ) ); #endif - } - if ( dt2 > 0. ) { + } + if ( dt2 > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[1], axom::utilities::min( 1.e6, dt2 ) ); + RAJA::atomicMin( &dt_temp[1], axom::utilities::min( fallback_dt, dt2 ) ); #else - dt_temp[1] = axom::utilities::min(dt_temp[1], axom::utilities::min(1.e6, dt2)); + dt_temp[1] = axom::utilities::min( dt_temp[1], axom::utilities::min( fallback_dt, dt2 ) ); #endif - } - if ( dt1 < 0. || dt2 < 0. ) { + } + if ( dt1 < 0. || dt2 < 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMax( &msg[3], static_cast( true ) ); + RAJA::atomicMax( &msg[3], static_cast( true ) ); #else - msg[3] = true; + msg[3] = true; #endif - } - - } // end check 2 - - //////////////////////////////////////////////////////////////////////// - // 3. Contact Dynamics and Stability Limits // - // // - // These limits constrain the explicit time integration based on // - // the added stiffness of active contact springs. // - //////////////////////////////////////////////////////////////////////// - - if ( plane.m_inContact && cs_view.enableTimestepStabilityLimits() ) { - RealT stiffness1 = 0.; - RealT stiffness2 = 0.; - auto& enforcement_options = cs_view.getEnforcementOptions(); - const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; - RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; - RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; - - ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, - pen_scale1, mesh1.getElementData().m_mat_mod[index1], - mesh1.getElementData().m_penalty_stiffness, stiffness1 ); - ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, - pen_scale2, mesh2.getElementData().m_mat_mod[index2], - mesh2.getElementData().m_penalty_stiffness, stiffness2 ); - - if ( stiffness1 > 0. && stiffness2 > 0. ) { - RealT penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); - if ( penalty_stiff_per_area > 0. ) { - // Rescale the host timestep by the added interface stiffness. - RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) ); - RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) ); - RealT dt_crit = alpha * dt * axom::utilities::min( f_scale1, f_scale2 ); - - if ( dt_crit > 0. ) { + } + + } // end check 2 + + //////////////////////////////////////////////////////////////////////// + // 3. Contact Dynamics and Stability Limits // + // // + // These limits constrain the explicit time integration based on // + // the added stiffness of active contact springs. // + //////////////////////////////////////////////////////////////////////// + + if ( plane.m_inContact && cs_view.enableTimestepStabilityLimits() ) { + RealT stiffness1 = 0.; + RealT stiffness2 = 0.; + auto& enforcement_options = cs_view.getEnforcementOptions(); + const PenaltyEnforcementOptions& pen_enfrc_options = enforcement_options.penalty_options; + RealT pen_scale1 = mesh1.getElementData().m_penalty_scale; + RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; + RealT mat_mod1 = 0.; + RealT mat_mod2 = 0.; + + if ( pen_enfrc_options.kinematic_calculation == KINEMATIC_ELEMENT ) { + mat_mod1 = mesh1.getElementData().m_mat_mod[index1]; + mat_mod2 = mesh2.getElementData().m_mat_mod[index2]; + } + + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, + pen_scale1, mat_mod1, mesh1.getElementData().m_penalty_stiffness, stiffness1 ); + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, + pen_scale2, mat_mod2, mesh2.getElementData().m_penalty_stiffness, stiffness2 ); + + if ( stiffness1 > 0. && stiffness2 > 0. ) { + RealT penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); + if ( penalty_stiff_per_area > 0. ) { + // Rescale the host timestep by the added interface stiffness. + RealT f_scale1 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness1 ) ); + RealT f_scale2 = 1.0 / std::sqrt( 1.0 + ( penalty_stiff_per_area / stiffness2 ) ); + RealT dt_crit = alpha * dt * axom::utilities::min( f_scale1, f_scale2 ); + + if ( dt_crit > 0. ) { #ifdef TRIBOL_USE_RAJA - RAJA::atomicMin( &dt_temp[2], axom::utilities::min( dt_crit, 1.e6 ) ); + RAJA::atomicMin( &dt_temp[2], axom::utilities::min( dt_crit, fallback_dt ) ); #else - dt_temp[2] = axom::utilities::min( dt_temp[2], axom::utilities::min( dt_crit, 1.e6 ) ); + dt_temp[2] = axom::utilities::min( dt_temp[2], axom::utilities::min( dt_crit, fallback_dt ) ); #endif - } - } - } - } // end contact dynamics and stability limits - } ); + } + } + } + } // end contact dynamics and stability limits + } ); // print general messages once // Can we output this message on root? SRW diff --git a/src/tribol/physics/CommonPlane.cpp b/src/tribol/physics/CommonPlane.cpp index c3b322d4..d976f04a 100644 --- a/src/tribol/physics/CommonPlane.cpp +++ b/src/tribol/physics/CommonPlane.cpp @@ -169,23 +169,28 @@ int ApplyNormal( CouplingScheme* cs ) RealT pen_scale2 = mesh2.getElementData().m_penalty_scale; RealT stiffness1 = 0.; RealT stiffness2 = 0.; - RealT t1 = mesh1.getElementData().m_thickness[index1]; - RealT t2 = mesh2.getElementData().m_thickness[index2]; - - ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, pen_scale1, - mesh1.getElementData().m_mat_mod[index1], mesh1.getElementData().m_penalty_stiffness, - stiffness1 ); - ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, pen_scale2, - mesh2.getElementData().m_mat_mod[index2], mesh2.getElementData().m_penalty_stiffness, - stiffness2 ); + RealT t1 = 0.; + RealT t2 = 0.; + RealT mat_mod1 = 0.; + RealT mat_mod2 = 0.; if ( pen_enfrc_options.kinematic_calculation == KINEMATIC_ELEMENT ) { + t1 = mesh1.getElementData().m_thickness[index1]; + t2 = mesh2.getElementData().m_thickness[index2]; + mat_mod1 = mesh1.getElementData().m_mat_mod[index1]; + mat_mod2 = mesh2.getElementData().m_mat_mod[index2]; if ( t1 <= 0. || t2 <= 0. ) { neg_thickness[0] = true; err[0] = 1; + return; } } + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t1, pen_enfrc_options.tiny_length, pen_scale1, + mat_mod1, mesh1.getElementData().m_penalty_stiffness, stiffness1 ); + ComputeUncoupledStiffness( pen_enfrc_options.kinematic_calculation, t2, pen_enfrc_options.tiny_length, pen_scale2, + mat_mod2, mesh2.getElementData().m_penalty_stiffness, stiffness2 ); + penalty_stiff_per_area = ComputePenaltyStiffnessPerArea( stiffness1, stiffness2 ); //////////////////////////////////////////////////// From a58ea173bc7b11b99020bfafb55e15534fe6e80c Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Mon, 1 Jun 2026 14:11:23 -0700 Subject: [PATCH 4/5] add doxygen comment --- src/tribol/interface/tribol.hpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/tribol/interface/tribol.hpp b/src/tribol/interface/tribol.hpp index 98d75af4..9a7e55ff 100644 --- a/src/tribol/interface/tribol.hpp +++ b/src/tribol/interface/tribol.hpp @@ -225,6 +225,16 @@ 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 ); /** From 748bae8dcda20c3f69dc9e29908e79caa3390d73 Mon Sep 17 00:00:00 2001 From: "Eric B. Chin" Date: Mon, 1 Jun 2026 14:16:49 -0700 Subject: [PATCH 5/5] add a plan to change into a direct CFL calc --- ...common_plane_contact_mass_timestep_plan.md | 75 +++++++++++++++++++ 1 file changed, 75 insertions(+) create mode 100644 docs/common_plane_contact_mass_timestep_plan.md diff --git a/docs/common_plane_contact_mass_timestep_plan.md b/docs/common_plane_contact_mass_timestep_plan.md new file mode 100644 index 00000000..aea1fc0d --- /dev/null +++ b/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.