Tupek/contact adjoint

src/serac/physics/contact/contact_data.cpp

@@ -41,6 +41,15 @@ void ContactData::addContactInteraction(int interaction_id, const std::set<int>&
   }
 }
 
+void ContactData::reset()
+{
+  for (auto& interaction : interactions_) {
+    FiniteElementState zero = interaction.pressure();
+    zero                    = 0.0;
+    interaction.setPressure(zero);
+  }
+}
+
 void ContactData::update(int cycle, double time, double& dt)
 {
   cycle_ = cycle;
@@ -226,7 +235,7 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
   return block_J;
 }
 
-void ContactData::residualFunction(const mfem::Vector& u, mfem::Vector& r)
+void ContactData::residualFunction(const mfem::Vector& u_shape, const mfem::Vector& u, mfem::Vector& r)
 {
   const int disp_size = reference_nodes_->ParFESpace()->GetTrueVSize();
 
@@ -239,7 +248,8 @@ void ContactData::residualFunction(const mfem::Vector& u, mfem::Vector& r)
   mfem::Vector r_blk(r, 0, disp_size);
   mfem::Vector g_blk(r, disp_size, numPressureDofs());
 
-  setDisplacements(u_blk);
+  setDisplacements(u_shape, u_blk);
+
   // we need to call update first to update gaps
   update(cycle_, time_, dt_);
   // with updated gaps, we can update pressure for contact interactions with penalty enforcement
@@ -289,10 +299,14 @@ void ContactData::setPressures(const mfem::Vector& merged_pressures) const
   }
 }
 
-void ContactData::setDisplacements(const mfem::Vector& u)
+void ContactData::setDisplacements(const mfem::Vector& shape_u, const mfem::Vector& u)
 {
+  mfem::ParGridFunction prolonged_shape_disp{current_coords_};
   reference_nodes_->ParFESpace()->GetProlongationMatrix()->Mult(u, current_coords_);
+  reference_nodes_->ParFESpace()->GetProlongationMatrix()->Mult(shape_u, prolonged_shape_disp);
+
   current_coords_ += *reference_nodes_;
+  current_coords_ += prolonged_shape_disp;
 }
 
 void ContactData::updateDofOffsets() const
@@ -373,7 +387,10 @@ std::unique_ptr<mfem::BlockOperator> ContactData::mergedJacobian() const
   return std::make_unique<mfem::BlockOperator>(jacobian_offsets_);
 }
 
-void ContactData::residualFunction([[maybe_unused]] const mfem::Vector& u, [[maybe_unused]] mfem::Vector& r) {}
+void ContactData::residualFunction([[maybe_unused]] const mfem::Vector& u_shape, [[maybe_unused]] const mfem::Vector& u,
+                                   [[maybe_unused]] mfem::Vector& r)
+{
+}
 
 std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(mfem::HypreParMatrix* orig_J) const
 {
@@ -390,7 +407,10 @@ std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(mfem::HyprePa
 
 void ContactData::setPressures([[maybe_unused]] const mfem::Vector& true_pressures) const {}
 
-void ContactData::setDisplacements([[maybe_unused]] const mfem::Vector& true_displacement) {}
+void ContactData::setDisplacements([[maybe_unused]] const mfem::Vector& u_shape,
+                                   [[maybe_unused]] const mfem::Vector& true_displacement)
+{
+}
 
 #endif
 

src/serac/physics/contact/contact_data.hpp

@@ -74,6 +74,12 @@ class ContactData {
    */
   void update(int cycle, double time, double& dt);
 
+  /**
+   * @brief Updates the positions, forces, and Jacobian contributions associated with contact
+   *
+   */
+  void reset();
+
   /**
    * @brief Get the contact constraint residual (i.e. nodal forces) from all contact interactions
    *
@@ -125,6 +131,7 @@ class ContactData {
   /**
    * @brief Computes the residual including contact terms
    *
+   * @param [in] u_shape Shape displacement vector (size of [displacement] block)
    * @param [in] u Solution vector ([displacement; pressure] block vector)
    * @param [in,out] r Residual vector ([force; gap] block vector); takes in initialized residual force vector and adds
    * contact contributions
@@ -133,7 +140,7 @@ class ContactData {
    *
    * @note This method calls update() to compute residual and Jacobian contributions based on the current configuration
    */
-  void residualFunction(const mfem::Vector& u, mfem::Vector& r);
+  void residualFunction(const mfem::Vector& u_shape, const mfem::Vector& u, mfem::Vector& r);
 
   /**
    * @brief Computes the Jacobian including contact terms, given the non-contact Jacobian terms
@@ -161,9 +168,10 @@ class ContactData {
   /**
    * @brief Update the current coordinates based on the new displacement field
    *
+   * @param u_shape Shape displacement vector
    * @param u Current displacement dof values
    */
-  void setDisplacements(const mfem::Vector& u);
+  void setDisplacements(const mfem::Vector& u_shape, const mfem::Vector& u);
 
   /**
    * @brief Have there been contact interactions added?

src/serac/physics/solid_mechanics.cpp

@@ -17,7 +17,7 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
                        mfem::HypreParVector& accel_adjoint_load_vector, mfem::HypreParVector& adjoint_displacement_,
                        mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_,
                        mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_,
-                       mfem::HypreParVector& adjoint_essential, BoundaryConditionManager& bcs_,
+                       mfem::HypreParVector& adjoint_essential_, BoundaryConditionManager& bcs_,
                        mfem::Solver& lin_solver)
 {
   // there are hard-coded here for now
@@ -48,7 +48,7 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
                   implicit_sensitivity_displacement_start_of_step_);
 
   for (const auto& bc : bcs_.essentials()) {
-    bc.apply(*J_T, adjoint_rhs, adjoint_essential);
+    bc.apply(*J_T, adjoint_rhs, adjoint_essential_);
   }
 
   lin_solver.SetOperator(*J_T);

src/serac/physics/solid_mechanics.hpp

@@ -38,7 +38,7 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
                        mfem::HypreParVector& accel_adjoint_load_vector, mfem::HypreParVector& adjoint_displacement_,
                        mfem::HypreParVector& implicit_sensitivity_displacement_start_of_step_,
                        mfem::HypreParVector& implicit_sensitivity_velocity_start_of_step_,
-                       mfem::HypreParVector& adjoint_essential, BoundaryConditionManager& bcs_,
+                       mfem::HypreParVector& adjoint_essential_, BoundaryConditionManager& bcs_,
                        mfem::Solver& lin_solver);
 }  // namespace detail
 
@@ -1209,6 +1209,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
             // tracking strategy
             // See https://github.com/mfem/mfem/issues/3531
             r = res;
+
             r.SetSubVector(bcs_.allEssentialTrueDofs(), 0.0);
           },
 
@@ -1385,8 +1386,6 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   void reverseAdjointTimestep() override
   {
     SERAC_MARK_FUNCTION;
-    auto& lin_solver = nonlin_solver_->linearSolver();
-
     SLIC_ERROR_ROOT_IF(cycle_ <= min_cycle_,
                        "Maximum number of adjoint timesteps exceeded! The number of adjoint timesteps must equal the "
                        "number of forward timesteps");
@@ -1401,29 +1400,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
     displacement_ = end_step_solution.at("displacement");
 
     if (is_quasistatic_) {
-      auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
-                                    *parameters_[parameter_indices].state...);
-      J_.reset();
-      J_ = assemble(drdu);
-
-      auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
-
-      J_e_.reset();
-      J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_T);
-
-      auto& constrained_dofs = bcs_.allEssentialTrueDofs();
-
-      mfem::EliminateBC(*J_T, *J_e_, constrained_dofs, reactions_adjoint_bcs_, displacement_adjoint_load_);
-      for (int i = 0; i < constrained_dofs.Size(); i++) {
-        int j                         = constrained_dofs[i];
-        displacement_adjoint_load_[j] = reactions_adjoint_bcs_[j];
-      }
-
-      lin_solver.SetOperator(*J_T);
-      lin_solver.Mult(displacement_adjoint_load_, adjoint_displacement_);
-
-      // Reset the equation solver to use the full nonlinear residual operator.  MRT, is this needed?
-      nonlin_solver_->setOperator(*residual_with_bcs_);
+      quasiStaticAdjointSolve(dt_n_to_np1);
     } else {
       SLIC_ERROR_ROOT_IF(ode2_.GetTimestepper() != TimestepMethod::Newmark,
                          "Only Newmark implemented for transient adjoint solid mechanics.");
@@ -1444,6 +1421,7 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
                                                    *parameters_[parameter_indices].state...));
       std::unique_ptr<mfem::HypreParMatrix> m_mat(assemble(M));
 
+      auto& lin_solver = nonlin_solver_->linearSolver();
       solid_mechanics::detail::adjoint_integrate(
           dt_n_to_np1, dt_np1_to_np2, m_mat.get(), k_mat.get(), displacement_adjoint_load_, velocity_adjoint_load_,
           acceleration_adjoint_load_, adjoint_displacement_, implicit_sensitivity_displacement_start_of_step_,
@@ -1636,6 +1614,35 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
     nonlin_solver_->solve(displacement_);
   }
 
+  /// @brief Solve the Quasi-static adjoint linear
+  virtual void quasiStaticAdjointSolve(double /*dt*/)
+  {
+    auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
+                                  *parameters_[parameter_indices].state...);
+    J_.reset();
+    J_ = assemble(drdu);
+
+    auto J_T = std::unique_ptr<mfem::HypreParMatrix>(J_->Transpose());
+
+    J_e_.reset();
+    J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_T);
+
+    auto& constrained_dofs = bcs_.allEssentialTrueDofs();
+
+    mfem::EliminateBC(*J_T, *J_e_, constrained_dofs, reactions_adjoint_bcs_, displacement_adjoint_load_);
+    for (int i = 0; i < constrained_dofs.Size(); i++) {
+      int j                         = constrained_dofs[i];
+      displacement_adjoint_load_[j] = reactions_adjoint_bcs_[j];
+    }
+
+    auto& lin_solver = nonlin_solver_->linearSolver();
+    lin_solver.SetOperator(*J_T);
+    lin_solver.Mult(displacement_adjoint_load_, adjoint_displacement_);
+
+    // Reset the equation solver to use the full nonlinear residual operator.  MRT, is this needed?
+    nonlin_solver_->setOperator(*residual_with_bcs_);
+  }
+
   /**
    * @brief Calculate a list of constrained dofs in the true displacement vector from a function that
    * returns true if a physical coordinate is in the constrained set