Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Support for non-symmetric contact constraints #1315

Draft
wants to merge 17 commits into
base: develop
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
14 changes: 14 additions & 0 deletions cmake/thirdparty/SetupSeracThirdParty.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,20 @@ if (NOT SERAC_THIRD_PARTY_LIBRARIES_FOUND)
set(serac_device_depends "" CACHE STRING "" FORCE)
endif()

#------------------------------------------------------------------------------
# Enzyme
#------------------------------------------------------------------------------
if (ENZYME_DIR)
message(STATUS "Setting up external Enzyme TPL...")

set(Enzyme_DIR ${ENZYME_DIR})
find_dependency(Enzyme REQUIRED PATHS "${ENZYME_DIR}")

set(TRIBOL_USE_ENZYME TRUE)
else()
message(STATUS "Enzyme support is OFF")
endif()

#------------------------------------------------------------------------------
# Camp
#------------------------------------------------------------------------------
Expand Down
3 changes: 2 additions & 1 deletion examples/contact/beam_bending.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,8 @@ int main(int argc, char* argv[])
serac::ContactOptions contact_options{.method = serac::ContactMethod::SingleMortar,
.enforcement = serac::ContactEnforcement::Penalty,
.type = serac::ContactType::Frictionless,
.penalty = 1.0e3};
.penalty = 1.0e3,
.jacobian = serac::ContactJacobian::Exact};

serac::SolidMechanicsContact<p, dim, serac::Parameters<serac::L2<0>, serac::L2<0>>> solid_solver(
nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, "beam_mesh",
Expand Down
7 changes: 4 additions & 3 deletions examples/contact/ironing.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -47,15 +47,16 @@ int main(int argc, char* argv[])
#endif

serac::NonlinearSolverOptions nonlinear_options{.nonlin_solver = serac::NonlinearSolver::Newton,
.relative_tol = 1.0e-7,
.absolute_tol = 1.0e-7,
.relative_tol = 1.0e-10,
.absolute_tol = 1.0e-10,
.max_iterations = 200,
.print_level = 1};

serac::ContactOptions contact_options{.method = serac::ContactMethod::SingleMortar,
.enforcement = serac::ContactEnforcement::Penalty,
.type = serac::ContactType::TiedNormal,
.penalty = 1.0e3};
.penalty = 1.0e3,
.jacobian = serac::ContactJacobian::Approximate};

serac::SolidMechanicsContact<p, dim, serac::Parameters<serac::L2<0>, serac::L2<0>>> solid_solver(
nonlinear_options, linear_options, serac::solid_mechanics::default_quasistatic_options, name, "ironing_mesh",
Expand Down
12 changes: 12 additions & 0 deletions src/serac/physics/contact/contact_config.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,15 @@ enum class ContactType
Frictionless /**< Enforce gap >= 0, pressure <= 0, gap * pressure = 0 in the normal direction */
};

/**
* @brief Method for computing Jacobian of contact terms
*/
enum class ContactJacobian
{
Approximate, /**< Ignore higher order contributions to the Jacobian */
Exact /**< Compute exact Jacobian (needed for quadratic convergence with Newton) */
};

/**
* @brief Stores the options for a contact pair
*/
Expand All @@ -55,6 +64,9 @@ struct ContactOptions {

/// Penalty parameter (only used when enforcement == ContactEnforcement::Penalty)
double penalty = 1.0e3;

/// The method to use for Jacobian calculations
ContactJacobian jacobian = ContactJacobian::Approximate;
};

} // namespace serac
54 changes: 35 additions & 19 deletions src/serac/physics/contact/contact_data.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -114,7 +114,12 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
block_J->owns_blocks = true;
// rather than returning different blocks for each contact interaction with Lagrange multipliers, merge them all into
// a single block
mfem::Array2D<const mfem::HypreParMatrix*> constraint_matrices(static_cast<int>(interactions_.size()), 1);
mfem::Array2D<const mfem::HypreParMatrix*> dgdu_blocks(static_cast<int>(interactions_.size()), 1);
mfem::Array2D<const mfem::HypreParMatrix*> dfdp_blocks(1, static_cast<int>(interactions_.size()));
for (size_t i{0}; i < interactions_.size(); ++i) {
dgdu_blocks(static_cast<int>(i), 0) = nullptr;
dfdp_blocks(0, static_cast<int>(i)) = nullptr;
}

for (size_t i{0}; i < interactions_.size(); ++i) {
// this is the BlockOperator for one of the contact interactions
Expand All @@ -139,16 +144,20 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
}
// add the contact interaction's (other) contribution to df_(contact)/dx (for penalty) or to df_(contact)/dp and
// dg/dx (for Lagrange multipliers)
if (!interaction_J->IsZeroBlock(1, 0)) {
auto B = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(1, 0));
SLIC_ERROR_ROOT_IF(!B, "Only HypreParMatrix constraint matrix blocks are currently supported.");
// zero out rows not in the active set
B->EliminateRows(interactions_[i].inactiveDofs());
if (!interaction_J->IsZeroBlock(1, 0) && !interaction_J->IsZeroBlock(0, 1)) {
auto dgdu = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(1, 0));
auto dfdp = dynamic_cast<mfem::HypreParMatrix*>(&interaction_J->GetBlock(0, 1));
SLIC_ERROR_ROOT_IF(!dgdu, "Only HypreParMatrix constraint matrix blocks are currently supported.");
SLIC_ERROR_ROOT_IF(!dfdp, "Only HypreParMatrix constraint matrix blocks are currently supported.");
// zero out rows and cols not in the active set
auto inactive_dofs = interactions_[i].inactiveDofs();
dgdu->EliminateRows(inactive_dofs);
dfdp->EliminateCols(inactive_dofs);
if (interactions_[i].getContactOptions().enforcement == ContactEnforcement::Penalty) {
// compute contribution to df_(contact)/dx (the 0, 0 block) for penalty
std::unique_ptr<mfem::HypreParMatrix> BTB(
mfem::ParMult(std::unique_ptr<mfem::HypreParMatrix>(B->Transpose()).get(), B, true));
std::unique_ptr<mfem::HypreParMatrix> BTB(mfem::ParMult(dfdp, dgdu, true));
delete &interaction_J->GetBlock(1, 0);
delete &interaction_J->GetBlock(0, 1);
if (block_J->IsZeroBlock(0, 0)) {
mfem::Vector penalty(reference_nodes_->ParFESpace()->GetTrueVSize());
penalty = interactions_[i].getContactOptions().penalty;
Expand All @@ -159,16 +168,12 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
mfem::Add(1.0, static_cast<mfem::HypreParMatrix&>(block_J->GetBlock(0, 0)),
interactions_[i].getContactOptions().penalty, *BTB));
}
constraint_matrices(static_cast<int>(i), 0) = nullptr;
} else // enforcement == ContactEnforcement::LagrangeMultiplier
{
// compute contribution to off-diagonal blocks for Lagrange multiplier
constraint_matrices(static_cast<int>(i), 0) = static_cast<mfem::HypreParMatrix*>(B);
}
if (interaction_J->IsZeroBlock(0, 1)) {
SLIC_ERROR_ROOT("Only symmetric constraint matrices are currently supported.");
dgdu_blocks(static_cast<int>(i), 0) = dgdu;
dfdp_blocks(0, static_cast<int>(i)) = dfdp;
}
delete &interaction_J->GetBlock(0, 1);
if (!interaction_J->IsZeroBlock(1, 1)) {
// we track our own active set, so get rid of the tribol inactive dof block
delete &interaction_J->GetBlock(1, 1);
Expand All @@ -177,9 +182,14 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
}
if (haveLagrangeMultipliers()) {
// merge all of the contributions from all of the contact interactions
block_J->SetBlock(1, 0, mfem::HypreParMatrixFromBlocks(constraint_matrices));
block_J->SetBlock(1, 0, mfem::HypreParMatrixFromBlocks(dgdu_blocks));
// store the transpose explicitly (rather than as a TransposeOperator) for solvers that need HypreParMatrixs
block_J->SetBlock(0, 1, static_cast<mfem::HypreParMatrix&>(block_J->GetBlock(1, 0)).Transpose());
block_J->SetBlock(0, 1, mfem::HypreParMatrixFromBlocks(dfdp_blocks));
// explicitly delete the blocks
for (size_t i{0}; i < interactions_.size(); ++i) {
delete dgdu_blocks(static_cast<int>(i), 0);
delete dfdp_blocks(0, static_cast<int>(i));
}
// build I_(inactive): a diagonal matrix with ones on inactive dofs and zeros elsewhere
mfem::Array<const mfem::Array<int>*> inactive_tdofs_vector(static_cast<int>(interactions_.size()));
int inactive_tdofs_ct = 0;
Expand Down Expand Up @@ -216,9 +226,9 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
// if the size of ones is zero, SparseMatrix creates its own memory which it
// owns. explicitly prevent this...
inactive_diag.SetDataOwner(false);
auto& block_1_0 = static_cast<mfem::HypreParMatrix&>(block_J->GetBlock(1, 0));
auto block_1_1 = new mfem::HypreParMatrix(block_1_0.GetComm(), block_1_0.GetGlobalNumRows(),
block_1_0.GetRowStarts(), &inactive_diag);
auto block_1_1 =
new mfem::HypreParMatrix(mesh_.GetComm(), global_pressure_dof_offsets_[global_pressure_dof_offsets_.Size() - 1],
global_pressure_dof_offsets_, &inactive_diag);
block_1_1->SetOwnerFlags(3, 3, 1);
block_J->SetBlock(1, 1, block_1_1);
// end building I_(inactive)
Expand All @@ -241,10 +251,16 @@ void ContactData::residualFunction(const mfem::Vector& u, mfem::Vector& r)

setDisplacements(u_blk);
// we need to call update first to update gaps
for (auto& interaction : interactions_) {
interaction.evalJacobian(false);
}
update(cycle_, time_, dt_);
// with updated gaps, we can update pressure for contact interactions with penalty enforcement
setPressures(p_blk);
// call update again with the right pressures
for (auto& interaction : interactions_) {
interaction.evalJacobian(true);
}
update(cycle_, time_, dt_);

r_blk += forces();
Expand Down
36 changes: 31 additions & 5 deletions src/serac/physics/contact/contact_interaction.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -62,6 +62,18 @@ ContactInteraction::ContactInteraction(int interaction_id, const mfem::ParMesh&
pressure_space.GetRestrictionMatrix()->BooleanMult(dof_markers, tdof_markers);
mfem::FiniteElementSpace::MarkerToList(tdof_markers, inactive_tdofs_);
}

// make sure we can compute exact Jacobian if requested
if (getContactOptions().jacobian == ContactJacobian::Exact) {
#ifdef TRIBOL_USE_ENZYME
tribol::enableEnzyme(interaction_id, true);
tribol::registerMfemReferenceCoords(interaction_id, static_cast<const mfem::ParGridFunction&>(*mesh.GetNodes()));
#else
SLIC_INFO_ROOT(
"Exact contact Jacobian contributions require Enzyme support in Tribol. Using approximate Jacobian.");
contact_opts_.jacobian = ContactJacobian::Approximate;
#endif
}
}

FiniteElementDual ContactInteraction::forces() const
Expand Down Expand Up @@ -112,15 +124,20 @@ void ContactInteraction::setPressure(const FiniteElementState& pressure) const
tribol::getMfemPressure(getInteractionId()) = pressure.gridFunction();
}

const mfem::Array<int>& ContactInteraction::inactiveDofs() const
const mfem::Array<int>& ContactInteraction::inactiveDofs() const { return inactiveDofs(pressure()); }

const mfem::Array<int>& ContactInteraction::inactiveDofs(const FiniteElementState& pressure) const
{
if (getContactOptions().type == ContactType::Frictionless) {
auto p = pressure();
auto g = gaps();
std::vector<int> inactive_tdofs_vector;
inactive_tdofs_vector.reserve(static_cast<size_t>(p.Size()));
for (int d{0}; d < p.Size(); ++d) {
if (p[d] >= 0.0 && g[d] >= -1.0e-14) {
inactive_tdofs_vector.reserve(static_cast<size_t>(pressure.Size()));
for (int d{0}; d < pressure.Size(); ++d) {
// don't check pressure for penalty; it's just the gap * a constant
if (getContactOptions().enforcement == ContactEnforcement::Penalty && g[d] >= -1.0e-14) {
inactive_tdofs_vector.push_back(d);
} else if (getContactOptions().enforcement == ContactEnforcement::LagrangeMultiplier && pressure[d] >= 0.0 &&
g[d] >= -1.0e-14) {
inactive_tdofs_vector.push_back(d);
}
}
Expand All @@ -130,6 +147,15 @@ const mfem::Array<int>& ContactInteraction::inactiveDofs() const
return inactive_tdofs_;
}

void ContactInteraction::evalJacobian(bool eval) const
{
if (eval) {
tribol::setLagrangeMultiplierOptions(getInteractionId(), tribol::ImplicitEvalMode::MORTAR_RESIDUAL_JACOBIAN);
} else {
tribol::setLagrangeMultiplierOptions(getInteractionId(), tribol::ImplicitEvalMode::MORTAR_GAP);
}
}

tribol::ContactMethod ContactInteraction::getMethod() const
{
switch (contact_opts_.method) {
Expand Down
15 changes: 15 additions & 0 deletions src/serac/physics/contact/contact_interaction.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -135,6 +135,21 @@ class ContactInteraction {
*/
const mfem::Array<int>& inactiveDofs() const;

/**
* @brief List of pressure/gap DOFs that are not active
*
* @param pressure A custom pressure field
* @return Array of inactive DOFs
*/
const mfem::Array<int>& inactiveDofs(const FiniteElementState& pressure) const;

/**
* @brief Turn on/off Jacobian evaluation
*
* @param eval True if Jacobian should be evaluated, false otherwise
*/
void evalJacobian(bool eval) const;

private:
/**
* @brief Get the Tribol enforcement method given a serac enforcement method
Expand Down
34 changes: 34 additions & 0 deletions src/serac/physics/solid_mechanics_contact.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -165,6 +165,40 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
J_constraint_e_->SetBlock(0, 0, J_e_.get());
J_constraint_e_->SetBlock(1, 0, J_e_21_.get());

// temp debug prints begin
mfem::SparseMatrix J11_sparse;
J_->MergeDiagAndOffd(J11_sparse);
mfem::DenseMatrix J11_dense;
J11_sparse.ToDenseMatrix(J11_dense);
std::ofstream dense11_mat("dense_J11.mat");
J11_dense.PrintMatlab(dense11_mat);
dense11_mat.close();

mfem::SparseMatrix J12_sparse;
J_12_->MergeDiagAndOffd(J12_sparse);
mfem::DenseMatrix J12_dense;
J12_sparse.ToDenseMatrix(J12_dense);
std::ofstream dense12_mat("dense_J12.mat");
J12_dense.PrintMatlab(dense12_mat);
dense12_mat.close();

mfem::SparseMatrix J21_sparse;
J_21_->MergeDiagAndOffd(J21_sparse);
mfem::DenseMatrix J21_dense;
J21_sparse.ToDenseMatrix(J21_dense);
std::ofstream dense21_mat("dense_J21.mat");
J21_dense.PrintMatlab(dense21_mat);
dense21_mat.close();

mfem::SparseMatrix J22_sparse;
J_22_->MergeDiagAndOffd(J22_sparse);
mfem::DenseMatrix J22_dense;
J22_sparse.ToDenseMatrix(J22_dense);
std::ofstream dense22_mat("dense_J22.mat");
J22_dense.PrintMatlab(dense22_mat);
dense22_mat.close();
// temp debug prints end

J_operator_ = J_constraint_.get();
return *J_constraint_;
});
Expand Down
38 changes: 25 additions & 13 deletions src/serac/physics/tests/contact_beam.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,8 @@

namespace serac {

class ContactTest : public testing::TestWithParam<std::tuple<ContactEnforcement, ContactType, std::string>> {};
class ContactTest
: public testing::TestWithParam<std::tuple<ContactEnforcement, ContactType, ContactJacobian, std::string>> {};

TEST_P(ContactTest, beam)
{
Expand All @@ -33,7 +34,7 @@ TEST_P(ContactTest, beam)
MPI_Barrier(MPI_COMM_WORLD);

// Create DataStore
std::string name = "contact_beam_" + std::get<2>(GetParam());
std::string name = "contact_beam_" + std::get<3>(GetParam());
axom::sidre::DataStore datastore;
StateManager::initialize(datastore, name + "_data");

Expand All @@ -54,17 +55,18 @@ TEST_P(ContactTest, beam)
.absolute_tol = 1.0e-12,
.max_iterations = 200,
.print_level = 1};
#ifdef SERAC_USE_SUNDIALS
// KINFullStep is preferred, but has issues when active set is enabled
if (std::get<1>(GetParam()) == ContactType::TiedNormal) {
nonlinear_options.nonlin_solver = NonlinearSolver::KINFullStep;
}
#endif
// #ifdef SERAC_USE_SUNDIALS
// // KINFullStep is preferred, but has issues when active set is enabled
// if (std::get<1>(GetParam()) == ContactType::TiedNormal) {
// nonlinear_options.nonlin_solver = NonlinearSolver::KINFullStep;
// }
// #endif

ContactOptions contact_options{.method = ContactMethod::SingleMortar,
.enforcement = std::get<0>(GetParam()),
.type = std::get<1>(GetParam()),
.penalty = 1.0e2};
.penalty = 1.0e2,
.jacobian = std::get<2>(GetParam())};

SolidMechanicsContact<p, dim> solid_solver(nonlinear_options, linear_options,
solid_mechanics::default_quasistatic_options, name, "beam_mesh");
Expand Down Expand Up @@ -114,12 +116,22 @@ TEST_P(ContactTest, beam)
// NOTE: if Penalty is first and Lagrange Multiplier is second, SuperLU gives a zero diagonal error
INSTANTIATE_TEST_SUITE_P(
tribol, ContactTest,
testing::Values(std::make_tuple(ContactEnforcement::Penalty, ContactType::TiedNormal, "penalty_tiednormal"),
std::make_tuple(ContactEnforcement::Penalty, ContactType::Frictionless, "penalty_frictionless"),
testing::Values(std::make_tuple(ContactEnforcement::Penalty, ContactType::TiedNormal, ContactJacobian::Approximate,
"penalty_tiednormal_Japprox"),
std::make_tuple(ContactEnforcement::Penalty, ContactType::Frictionless,
ContactJacobian::Approximate, "penalty_frictionless_Japprox"),
std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactType::TiedNormal,
ContactJacobian::Approximate, "lagrange_multiplier_tiednormal_Japprox"),
std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactType::Frictionless,
ContactJacobian::Approximate, "lagrange_multiplier_frictionless_Japprox"),
std::make_tuple(ContactEnforcement::Penalty, ContactType::TiedNormal, ContactJacobian::Exact,
"penalty_tiednormal_Jexact"),
std::make_tuple(ContactEnforcement::Penalty, ContactType::Frictionless, ContactJacobian::Exact,
"penalty_frictionless_Jexact"),
std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactType::TiedNormal,
"lagrange_multiplier_tiednormal"),
ContactJacobian::Exact, "lagrange_multiplier_tiednormal_Jexact"),
std::make_tuple(ContactEnforcement::LagrangeMultiplier, ContactType::Frictionless,
"lagrange_multiplier_frictionless")));
ContactJacobian::Exact, "lagrange_multiplier_frictionless_Jexact")));

} // namespace serac

Expand Down
Loading
Loading