Tupek/contact adjoint

src/serac/physics/contact/contact_data.cpp

@@ -226,7 +226,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 +239,7 @@ 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
@@ -253,14 +253,8 @@ void ContactData::residualFunction(const mfem::Vector& u, mfem::Vector& r)
   g_blk.Set(1.0, mergedGaps(true));
 }
 
-std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(const mfem::Vector&   u,
-                                                                   mfem::HypreParMatrix* orig_J) const
+std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(mfem::HypreParMatrix* orig_J) const
 {
-  // u_const should not change in this method; const cast is to create vector views which are used to compute the
-  // (non-contact) Jacobian
-  auto&              u_const = const_cast<mfem::Vector&>(u);
-  const mfem::Vector u_blk(u_const, 0, reference_nodes_->ParFESpace()->GetTrueVSize());
-
   auto J_contact = mergedJacobian();
   if (J_contact->IsZeroBlock(0, 0)) {
     J_contact->SetBlock(0, 0, orig_J);
@@ -295,10 +289,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
@@ -379,16 +377,13 @@ 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) {}
-
-std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(const mfem::Vector&   u,
-                                                                   mfem::HypreParMatrix* orig_J) const
+void ContactData::residualFunction([[maybe_unused]] const mfem::Vector& u_shape, [[maybe_unused]] const mfem::Vector& u,
+                                   [[maybe_unused]] mfem::Vector& r)
 {
-  // u_const should not change in this method; const cast is to create vector views which are used to compute the
-  // (non-contact) Jacobian
-  auto&              u_const = const_cast<mfem::Vector&>(u);
-  const mfem::Vector u_blk(u_const, 0, reference_nodes_->ParFESpace()->GetTrueVSize());
+}
 
+std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(mfem::HypreParMatrix* orig_J) const
+{
   auto J_contact = mergedJacobian();
   if (J_contact->IsZeroBlock(0, 0)) {
     J_contact->SetBlock(0, 0, orig_J);
@@ -402,7 +397,10 @@ std::unique_ptr<mfem::BlockOperator> ContactData::jacobianFunction(const mfem::V
 
 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

@@ -121,20 +121,20 @@ 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
    */
-  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
    *
-   * @param u Solution vector ([displacement; pressure] block vector)
    * @param orig_J The non-contact terms of the Jacobian, not including essential boundary conditions
    * @return Jacobian with contact terms, not including essential boundary conditions
    */
-  std::unique_ptr<mfem::BlockOperator> jacobianFunction(const mfem::Vector& u, mfem::HypreParMatrix* orig_J) const;
+  std::unique_ptr<mfem::BlockOperator> jacobianFunction(mfem::HypreParMatrix* orig_J) const;
 
   /**
    * @brief Set the pressure field
@@ -152,9 +152,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

@@ -15,10 +15,9 @@ namespace detail {
 void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat, mfem::HypreParMatrix* k_mat,
                        mfem::HypreParVector& disp_adjoint_load_vector, mfem::HypreParVector& velo_adjoint_load_vector,
                        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::Solver& lin_solver)
+                       mfem::HypreParVector&     implicit_sensitivity_displacement_start_of_step_,
+                       mfem::HypreParVector&     implicit_sensitivity_velocity_start_of_step_,
+                       BoundaryConditionManager& bcs_, mfem::Solver& lin_solver)
 {
   // there are hard-coded here for now
   static constexpr double beta  = 0.25;
@@ -39,6 +38,10 @@ void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat,
   // recall that temperature_adjoint_load_vector and d_temperature_dt_adjoint_load_vector were already multiplied by
   // -1 above
 
+  // By default, use a homogeneous essential boundary condition
+  mfem::HypreParVector adjoint_essential(adjoint_displacement_);
+  adjoint_essential = 0.0;
+
   mfem::HypreParVector adjoint_rhs(accel_adjoint_load_vector);
   adjoint_rhs.Add(fac4 * dt_n, velo_adjoint_load_vector);
   adjoint_rhs.Add(fac3 * dt_n * dt_n, disp_adjoint_load_vector);

src/serac/physics/solid_mechanics.hpp

@@ -36,10 +36,9 @@ namespace detail {
 void adjoint_integrate(double dt_n, double dt_np1, mfem::HypreParMatrix* m_mat, mfem::HypreParMatrix* k_mat,
                        mfem::HypreParVector& disp_adjoint_load_vector, mfem::HypreParVector& velo_adjoint_load_vector,
                        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::Solver& lin_solver);
+                       mfem::HypreParVector&     implicit_sensitivity_displacement_start_of_step_,
+                       mfem::HypreParVector&     implicit_sensitivity_velocity_start_of_step_,
+                       BoundaryConditionManager& bcs_, mfem::Solver& lin_solver);
 }  // namespace detail
 
 /**
@@ -1182,6 +1181,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);
           },
 
@@ -1340,12 +1340,6 @@ class SolidMechanics<order, dim, Parameters<parameter_space...>, std::integer_se
   /// @overload
   void reverseAdjointTimestep() override
   {
-    auto& lin_solver = nonlin_solver_->linearSolver();
-
-    // By default, use a homogeneous essential boundary condition
-    mfem::HypreParVector adjoint_essential(displacement_adjoint_load_);
-    adjoint_essential = 0.0;
-
     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");
@@ -1359,20 +1353,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...);
-      auto jacobian  = assemble(drdu);
-      auto J_T       = std::unique_ptr<mfem::HypreParMatrix>(jacobian->Transpose());
-
-      for (const auto& bc : bcs_.essentials()) {
-        bc.apply(*J_T, displacement_adjoint_load_, adjoint_essential);
-      }
-
-      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.");
@@ -1393,10 +1374,11 @@ 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_,
-          implicit_sensitivity_velocity_start_of_step_, adjoint_essential, bcs_, lin_solver);
+          implicit_sensitivity_velocity_start_of_step_, bcs_, lin_solver);
     }
 
     time_end_step_ = time_;
@@ -1629,6 +1611,30 @@ 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*/)
+  {
+    // By default, use a homogeneous essential boundary condition
+    mfem::HypreParVector adjoint_essential(displacement_adjoint_load_);
+    adjoint_essential = 0.0;
+
+    auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
+                                  *parameters_[parameter_indices].state...);
+    auto jacobian  = assemble(drdu);
+    auto J_T       = std::unique_ptr<mfem::HypreParMatrix>(jacobian->Transpose());
+
+    for (const auto& bc : bcs_.essentials()) {
+      bc.apply(*J_T, displacement_adjoint_load_, adjoint_essential);
+    }
+
+    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

src/serac/physics/solid_mechanics_contact.hpp

@@ -119,7 +119,8 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
       // See https://github.com/mfem/mfem/issues/3531
       mfem::Vector r_blk(r, 0, displacement_.Size());
       r_blk = res;
-      contact_.residualFunction(u, r);
+
+      contact_.residualFunction(shape_displacement_, u, r);
       r_blk.SetSubVector(bcs_.allEssentialTrueDofs(), 0.0);
     };
     // This if-block below breaks up building the Jacobian operator depending if there is Lagrange multiplier
@@ -138,7 +139,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
             J_             = assemble(drdu);
 
             // create block operator holding jacobian contributions
-            J_constraint_ = contact_.jacobianFunction(u, J_.release());
+            J_constraint_ = contact_.jacobianFunction(J_.release());
 
             // take ownership of blocks
             J_constraint_->owns_blocks = false;
@@ -174,7 +175,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
             J_             = assemble(drdu);
 
             // get 11-block holding jacobian contributions
-            auto block_J         = contact_.jacobianFunction(u, J_.release());
+            auto block_J         = contact_.jacobianFunction(J_.release());
             block_J->owns_blocks = false;
             J_ = std::unique_ptr<mfem::HypreParMatrix>(static_cast<mfem::HypreParMatrix*>(&block_J->GetBlock(0, 0)));
 
@@ -253,6 +254,52 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
     forces_.SetVector(contact_.forces(), 0);
   }
 
+  /// @brief Solve the Quasi-static Newton system
+  void quasiStaticAdjointSolve(double /*dt*/) override
+  {
+    SLIC_ERROR_ROOT_IF(contact_.haveLagrangeMultipliers(),
+                       "Lagrange multiplier contact does not currently support sensitivities/adjoints.");
+
+    // By default, use a homogeneous essential boundary condition
+    mfem::HypreParVector adjoint_essential(displacement_adjoint_load_);
+    adjoint_essential = 0.0;
+
+    auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_,
+                                  *parameters_[parameter_indices].state...);
+    auto jacobian  = assemble(drdu);
+
+    auto block_J         = contact_.jacobianFunction(jacobian.release());
+    block_J->owns_blocks = false;
+    jacobian = std::unique_ptr<mfem::HypreParMatrix>(static_cast<mfem::HypreParMatrix*>(&block_J->GetBlock(0, 0)));
+    auto J_T = std::unique_ptr<mfem::HypreParMatrix>(jacobian->Transpose());
+
+    for (const auto& bc : bcs_.essentials()) {
+      bc.apply(*J_T, displacement_adjoint_load_, adjoint_essential);
+    }
+
+    auto& lin_solver = nonlin_solver_->linearSolver();
+    lin_solver.SetOperator(*J_T);
+    lin_solver.Mult(displacement_adjoint_load_, adjoint_displacement_);
+  }
+
+  /// @overload
+  FiniteElementDual& computeTimestepShapeSensitivity() override
+  {
+    auto drdshape =
+        serac::get<DERIVATIVE>((*residual_)(time_end_step_, differentiate_wrt(shape_displacement_), displacement_,
+                                            acceleration_, *parameters_[parameter_indices].state...));
+
+    auto drdshape_mat = assemble(drdshape);
+
+    auto block_J         = contact_.jacobianFunction(drdshape_mat.release());
+    block_J->owns_blocks = false;
+    drdshape_mat = std::unique_ptr<mfem::HypreParMatrix>(static_cast<mfem::HypreParMatrix*>(&block_J->GetBlock(0, 0)));
+
+    drdshape_mat->MultTranspose(adjoint_displacement_, *shape_displacement_sensitivity_);
+
+    return *shape_displacement_sensitivity_;
+  }
+
   using BasePhysics::bcs_;
   using BasePhysics::cycle_;
   using BasePhysics::duals_;
@@ -261,19 +308,23 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
   using BasePhysics::name_;
   using BasePhysics::parameters_;
   using BasePhysics::shape_displacement_;
+  using BasePhysics::shape_displacement_sensitivity_;
   using BasePhysics::states_;
   using BasePhysics::time_;
   using SolidMechanicsBase::acceleration_;
+  using SolidMechanicsBase::adjoint_displacement_;
   using SolidMechanicsBase::d_residual_d_;
   using SolidMechanicsBase::DERIVATIVE;
   using SolidMechanicsBase::displacement_;
+  using SolidMechanicsBase::displacement_adjoint_load_;
   using SolidMechanicsBase::du_;
   using SolidMechanicsBase::J_;
   using SolidMechanicsBase::J_e_;
   using SolidMechanicsBase::nonlin_solver_;
   using SolidMechanicsBase::ode_time_point_;
   using SolidMechanicsBase::residual_;
   using SolidMechanicsBase::residual_with_bcs_;
+  using SolidMechanicsBase::time_end_step_;
   using SolidMechanicsBase::warmStartDisplacement;
 
   /// Pointer to the Jacobian operator (J_ if no Lagrange multiplier contact, J_constraint_ otherwise)

src/serac/physics/tests/CMakeLists.txt

@@ -40,7 +40,7 @@ set(physics_parallel_test_sources
     solid_reaction_adjoint.cpp
     thermal_nonlinear_solve.cpp
     )
-blt_list_append(TO physics_parallel_test_sources ELEMENTS contact_patch.cpp contact_beam.cpp IF TRIBOL_FOUND)
+blt_list_append(TO physics_parallel_test_sources ELEMENTS contact_patch.cpp contact_beam.cpp contact_solid_adjoint.cpp IF TRIBOL_FOUND)
 
 serac_add_tests(SOURCES       ${physics_parallel_test_sources}
                 DEPENDS_ON    ${physics_test_depends}