(WIP) -- Contact-Homotopy Updates (Cache results of Tribol calls, Non-tied contact, etc).

examples/inertia_relief/inertia_relief_example.cpp

@@ -427,9 +427,9 @@ mfem::Vector InertialReliefProblem::constraint(const mfem::Vector& u, bool fresh
 // Jacobian of the constraint
 mfem::HypreParMatrix* InertialReliefProblem::constraintJacobian(const mfem::Vector& u, bool fresh_evaluation)
 {
-  int dim_constraints = GetMultiplierDim();
-  int glbdim_displacement = GetGlobalDisplacementDim();
   if (fresh_evaluation) {
+    int dim_constraints = GetMultiplierDim();
+    int glbdim_displacement = GetGlobalDisplacementDim();
     obj_states_[FIELD::DISP]->Set(1.0, u);
     // dense rows
     int nentries = glbdim_displacement;

src/smith/physics/constraint.hpp

@@ -54,11 +54,14 @@ class Constraint {
    * @param fields vector of smith::FiniteElementState*
    * @param direction index for which field to take the gradient with respect to
    * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation
+   * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached
+   * evaluation
    * @return std::unique_ptr<mfem::HypreParMatrix>
    */
   virtual std::unique_ptr<mfem::HypreParMatrix> jacobian(double time, double dt,
                                                          const std::vector<ConstFieldPtr>& fields, int direction,
-                                                         bool fresh_evaluation = true) const = 0;
+                                                         bool fresh_evaluation = true,
+                                                         bool fresh_derivative = true) const = 0;
 
   /** @brief Virtual interface for computing constraint Jacobian_tilde from a vector of smith::FiniteElementState*
    *         Jacobian_tilde is an optional approximation of the true Jacobian
@@ -68,13 +71,16 @@ class Constraint {
    * @param fields vector of smith::FiniteElementState*
    * @param direction index for which field to take the gradient with respect to
    * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation
+   * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached
+   * evaluation
    * @return std::unique_ptr<mfem::HypreParMatrix>
    */
   virtual std::unique_ptr<mfem::HypreParMatrix> jacobian_tilde(double time, double dt,
                                                                const std::vector<ConstFieldPtr>& fields, int direction,
-                                                               bool fresh_evaluation = true) const
+                                                               bool fresh_evaluation = true,
+                                                               bool fresh_derivative = true) const
   {
-    return jacobian(time, dt, fields, direction, fresh_evaluation);
+    return jacobian(time, dt, fields, direction, fresh_evaluation, fresh_derivative);
   };
 
   /** @brief Virtual interface for computing residual contribution Jacobian_tilde^(Transpose) * (Lagrange multiplier)
@@ -86,13 +92,16 @@ class Constraint {
    * @param multipliers mfem::Vector of Lagrange multipliers
    * @param direction index for which field to take the gradient with respect to
    * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation
+   * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached
+   * evaluation
    * @return std::Vector
    */
   virtual mfem::Vector residual_contribution(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                                              const mfem::Vector& multipliers, int direction,
-                                             bool fresh_evaluation = true) const
+                                             bool fresh_evaluation = true, bool fresh_derivative = true) const
   {
-    std::unique_ptr<mfem::HypreParMatrix> jac = jacobian_tilde(time, dt, fields, direction, fresh_evaluation);
+    std::unique_ptr<mfem::HypreParMatrix> jac =
+        jacobian_tilde(time, dt, fields, direction, fresh_evaluation, fresh_derivative);
     mfem::Vector y(jac->Width());
     y = 0.0;
     SLIC_ERROR_ROOT_IF(jac->Height() != multipliers.Size(), "Incompatible matrix and vector sizes.");
@@ -109,12 +118,15 @@ class Constraint {
    * @param multipliers mfem::Vector of Lagrange multipliers
    * @param direction index for which field to take the gradient with respect to
    * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation
+   * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached
+   * evaluation
    * @return std::unique_ptr<mfem::HypreParMatrix>
    */
   virtual std::unique_ptr<mfem::HypreParMatrix> residual_contribution_jacobian(
       [[maybe_unused]] double time, [[maybe_unused]] double dt,
       [[maybe_unused]] const std::vector<ConstFieldPtr>& fields, [[maybe_unused]] const mfem::Vector& multipliers,
-      [[maybe_unused]] int direction, [[maybe_unused]] bool fresh_evaluation = true) const
+      [[maybe_unused]] int direction, [[maybe_unused]] bool fresh_evaluation = true,
+      [[maybe_unused]] bool fresh_derivative = true) const
   {
     SLIC_ERROR_ROOT(axom::fmt::format("Base class must override residual_contribution_jacobian before usage"));
     std::unique_ptr<mfem::HypreParMatrix> res_contr_jacobian = nullptr;

src/smith/physics/contact_constraint.hpp

@@ -44,41 +44,23 @@ enum ContactFields
 
 /** @brief Interface for extracting the iblock, jblock block from a std::unique_ptr<mfem::BlockOperator>
  *         said block is returned as a std::unique_ptr<mfem::HypreParMatrix> if possible
- *         All other blocks are deleted
  *
  * @param block_operator the block operator
  * @param iblock row block index
  * @param jblock column block index
  * @return std::unique_ptr<mfem::HypreParMatrix> The requested block of block_operator
  */
-static std::unique_ptr<mfem::HypreParMatrix> safelyObtainBlock(mfem::BlockOperator* block_operator, int iblock,
-                                                               int jblock, bool own_blocks = false)
+static std::unique_ptr<mfem::HypreParMatrix> obtainBlock(mfem::BlockOperator* block_operator, int iblock, int jblock)
 {
   SLIC_ERROR_IF(iblock < 0 || jblock < 0, "block indicies must be non-negative");
   SLIC_ERROR_IF(iblock > block_operator->NumRowBlocks() || jblock > block_operator->NumColBlocks(),
                 "one or more block indicies are too large and the requested block does not exist");
   SLIC_ERROR_IF(block_operator->IsZeroBlock(iblock, jblock), "attempting to extract a null block");
-  block_operator->owns_blocks = false;
-  for (int i = 0; i < block_operator->NumRowBlocks(); i++) {
-    for (int j = 0; j < block_operator->NumColBlocks(); j++) {
-      if (i == iblock && j == jblock) {
-        continue;
-      }
-      if (!block_operator->IsZeroBlock(i, j) && !own_blocks) {
-        delete &block_operator->GetBlock(i, j);
-      }
-    }
-  }
   auto Ablk = dynamic_cast<mfem::HypreParMatrix*>(&block_operator->GetBlock(iblock, jblock));
   SLIC_ERROR_IF(!Ablk, "failed cast block to mfem::HypreParMatrix");
-  if (own_blocks) {
-    // deep copy --> unique_ptr
-    auto Ablk_unique = std::make_unique<mfem::HypreParMatrix>(*Ablk);
-    return Ablk_unique;
-  } else {
-    std::unique_ptr<mfem::HypreParMatrix> Ablk_unique(Ablk);
-    return Ablk_unique;
-  }
+  // deep copy --> unique_ptr
+  auto Ablk_unique = std::make_unique<mfem::HypreParMatrix>(*Ablk);
+  return Ablk_unique;
 };
 
 class FiniteElementState;
@@ -111,6 +93,7 @@ class ContactConstraint : public Constraint {
     contact_opts_.enforcement = ContactEnforcement::LagrangeMultiplier;
     contact_.addContactInteraction(interaction_id, bdry_attr_surf1, bdry_attr_surf2, contact_opts_);
     interaction_id_ = interaction_id;
+    tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP);
   }
 
   /// @brief destructor
@@ -128,20 +111,18 @@ class ContactConstraint : public Constraint {
   mfem::Vector evaluate(double time, double dt, const std::vector<ConstFieldPtr>& fields,
                         bool fresh_evaluation = true) const override
   {
-    contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP);
-
     if (fresh_evaluation) {
+      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
       // note: Tribol does not use cycle.
       int cycle = 0;
       contact_.update(cycle, time, dt);
-      pressures_set_ = false;
+      auto gaps_hpv = contact_.mergedGaps(false);
+      // Note: this copy is needed to prevent the HypreParVector pointer from going out of scope.  see
+      // https://github.com/mfem/mfem/issues/5029
+      cached_gap_.SetSize(gaps_hpv.Size());
+      cached_gap_.Set(1.0, gaps_hpv);
     }
-    auto gaps_hpv = contact_.mergedGaps(false);
-    // Note: this copy is needed to prevent the HypreParVector pointer from going out of scope.  see
-    // https://github.com/mfem/mfem/issues/5029
-    mfem::Vector gaps = gaps_hpv;
-    return gaps;
+    return cached_gap_;
   };
 
   /** @brief Interface for computing contact gap constraint Jacobian from a vector of smith::FiniteElementState*
@@ -151,24 +132,27 @@ class ContactConstraint : public Constraint {
    * @param fields vector of smith::FiniteElementState*
    * @param direction index for which field to take the gradient with respect to
    * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation
+   * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached
+   * evaluation
    * @return std::unique_ptr<mfem::HypreParMatrix> The true Jacobian
    */
   std::unique_ptr<mfem::HypreParMatrix> jacobian(double time, double dt, const std::vector<ConstFieldPtr>& fields,
-                                                 int direction, bool fresh_evaluation = true) const override
+                                                 int direction, bool fresh_evaluation = true,
+                                                 [[maybe_unused]] bool fresh_derivative = true) const override
   {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
-    contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN);
-
+    /* If this is not a new evaluation point e.g.,
+       if the constraint was evaluated at this point
+       then the derivatives have already been computed
+       and we do not need to recompute them.*/
     if (fresh_evaluation) {
+      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
       int cycle = 0;
       contact_.update(cycle, time, dt);
       J_contact_ = contact_.mergedJacobian();
-      pressures_set_ = false;
     }
-    int iblock = 1;
-    int jblock = 0;
-    auto dgdu = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_);
+    // obtain (1, 0) block entry from the 2 x 2 block contact linear system
+    auto dgdu = obtainBlock(J_contact_.get(), 1, 0);
     return dgdu;
   };
 
@@ -181,35 +165,32 @@ class ContactConstraint : public Constraint {
    * @param multipliers mfem::Vector of Lagrange multipliers
    * @param direction index for which field to take the gradient with respect to
    * @param fresh_evaluation boolean indicating if we re-evaluate or use previously cached evaluation
+   * @param fresh_derivative boolean indicating with fresh_evaluation if we re-evaluate or use previously cached
+   * evaluation
    * @return std::Vector
    */
   mfem::Vector residual_contribution(double time, double dt, const std::vector<ConstFieldPtr>& fields,
-                                     const mfem::Vector& multipliers, int direction,
-                                     bool fresh_evaluation = true) const override
+                                     const mfem::Vector& multipliers, int direction, bool fresh_evaluation = true,
+                                     bool fresh_derivative = true) const override
   {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
-    contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_GAP);
-
+    // TODO: include more fine scale caching control
     int cycle = 0;
-    if (fresh_evaluation) {
-      // we need to call update first to update gaps
+    if (fresh_evaluation || fresh_derivative) {
+      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
+      // first update gaps
       for (auto& interaction : contact_.getContactInteractions()) {
         interaction.evalJacobian(false);
       }
       contact_.update(cycle, time, dt);
-      pressures_set_ = false;
     }
-    if (!pressures_set_) {
-      // with updated gaps, we can update pressure for contact interactions with penalty enforcement
-      contact_.setPressures(multipliers);
-      // call update again with the right pressures
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(true);
-      }
-      contact_.update(cycle, time, dt);
-      pressures_set_ = true;
+    // with updated gaps, then update pressure for contact interactions with penalty enforcement
+    contact_.setPressures(multipliers);
+    // call update again with the right pressures
+    for (auto& interaction : contact_.getContactInteractions()) {
+      interaction.evalJacobian(true);
     }
+    contact_.update(cycle, time, dt);
 
     return contact_.forces();
   };
@@ -228,35 +209,32 @@ class ContactConstraint : public Constraint {
   std::unique_ptr<mfem::HypreParMatrix> residual_contribution_jacobian(double time, double dt,
                                                                        const std::vector<ConstFieldPtr>& fields,
                                                                        const mfem::Vector& multipliers, int direction,
-                                                                       bool fresh_evaluation = true) const override
+                                                                       bool fresh_evaluation = true,
+                                                                       bool fresh_derivative = true) const override
   {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
-    contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN);
 
     int cycle = 0;
-    if (fresh_evaluation) {
-      // we need to call update first to update gaps
+    // TODO: include more fine scale caching control
+    if (fresh_evaluation || fresh_derivative) {
+      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
+      // first update gaps
       for (auto& interaction : contact_.getContactInteractions()) {
         interaction.evalJacobian(false);
       }
       contact_.update(cycle, time, dt);
-      pressures_set_ = false;
     }
-    if (!pressures_set_) {
-      // with updated gaps, we can update pressure for contact interactions with penalty enforcement
-      contact_.setPressures(multipliers);
-      // call update again with the right pressures
-      for (auto& interaction : contact_.getContactInteractions()) {
-        interaction.evalJacobian(true);
-      }
-      contact_.update(cycle, time, dt);
-      J_contact_ = contact_.mergedJacobian();
-      pressures_set_ = true;
+    // with updated gaps, we can update pressure for contact interactions with penalty enforcement
+    contact_.setPressures(multipliers);
+    // call update again with the right pressures
+    for (auto& interaction : contact_.getContactInteractions()) {
+      interaction.evalJacobian(true);
     }
-    int iblock = 0;
-    int jblock = 0;
-    auto Hessian = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_);
+    contact_.update(cycle, time, dt);
+    J_contact_ = contact_.mergedJacobian();
+
+    // obtain (0, 0) block entry from the 2 x 2 block contact linear system
+    auto Hessian = obtainBlock(J_contact_.get(), 0, 0);
     return Hessian;
   };
 
@@ -270,22 +248,22 @@ class ContactConstraint : public Constraint {
    * @return std::unique_ptr<mfem::HypreParMatrix>
    */
   std::unique_ptr<mfem::HypreParMatrix> jacobian_tilde(double time, double dt, const std::vector<ConstFieldPtr>& fields,
-                                                       int direction, bool fresh_evaluation = true) const override
+                                                       int direction, bool fresh_evaluation = true,
+                                                       [[maybe_unused]] bool fresh_derivative = true) const override
   {
     SLIC_ERROR_IF(direction != ContactFields::DISP, "requesting a non displacement-field derivative");
-    contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
-    tribol::setLagrangeMultiplierOptions(interaction_id_, tribol::ImplicitEvalMode::MORTAR_JACOBIAN);
 
     int cycle = 0;
     if (fresh_evaluation) {
+      contact_.setDisplacements(*fields[ContactFields::SHAPE], *fields[ContactFields::DISP]);
       contact_.update(cycle, time, dt);
       J_contact_.reset();
       J_contact_ = contact_.mergedJacobian();
-      pressures_set_ = false;
     }
-    int iblock = 0;
-    int jblock = 1;
-    auto dgduT = safelyObtainBlock(J_contact_.get(), iblock, jblock, own_blocks_);
+    // obtain (0, 1) block entry from the 2 x 2 block contact linear system
+    auto dgduT = obtainBlock(J_contact_.get(), 0, 1);
+    // TODO: do we really need to do transpose here? Is another transpose done in the
+    // other HomotopySolver or problem class codes? Would a utility function be helpful?
     std::unique_ptr<mfem::HypreParMatrix> dgdu(dgduT->Transpose());
     return dgdu;
   };
@@ -314,14 +292,16 @@ class ContactConstraint : public Constraint {
   mutable std::unique_ptr<mfem::BlockOperator> J_contact_;
 
   /**
-   * @brief own_blocks_ temporary boolean
+   * @brief own_blocks_ boolean, specifying whether the ContactConstraint class
+   *                    will own the derivative blocks of J_contact_ and be responsible
+   *                    for managing the associated memory.
    */
   const bool own_blocks_ = true;
 
   /**
-   * @brief pressures_set_ are the Lagrange multipliers set
+   * @brief cached_gap_ gap function cached in order to not redo any unnecessary computations
    */
-  mutable bool pressures_set_ = false;
+  mutable mfem::Vector cached_gap_;
 };
 
 }  // namespace smith