Tupek/enable reaction adjoints with contact

src/smith/physics/solid_mechanics_contact.hpp

@@ -393,10 +393,6 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
     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_, BasePhysics::shapeDisplacement(), differentiate_wrt(displacement_),
                                   acceleration_, *parameters_[parameter_indices].state...);
 
@@ -406,7 +402,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
     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);
+      bc.apply(*J_T, displacement_adjoint_load_, reactions_adjoint_bcs_);
     }
 
     auto& lin_solver = nonlin_solver_->linearSolver();
@@ -450,6 +446,7 @@ class SolidMechanicsContact<order, dim, Parameters<parameter_space...>,
   using SolidMechanicsBase::J_;
   using SolidMechanicsBase::J_e_;
   using SolidMechanicsBase::nonlin_solver_;
+  using SolidMechanicsBase::reactions_adjoint_bcs_;
   using SolidMechanicsBase::residual_;
   using SolidMechanicsBase::residual_with_bcs_;
   using SolidMechanicsBase::shape_displacement_dual_;

src/smith/physics/tests/contact_solid_adjoint.cpp

@@ -60,6 +60,68 @@ void computeStepAdjointLoad(const FiniteElementState& displacement, FiniteElemen
   d_qoi_d_displacement.Add(1.0, displacement);
 }
 
+FiniteElementState createReactionDirection(const BasePhysics& solid_solver, int direction, std::shared_ptr<Mesh> mesh)
+{
+  const FiniteElementDual& reactions = solid_solver.dual("reactions");
+
+  FiniteElementState reactionDirections(reactions.space(), "reaction_directions");
+  reactionDirections = 0.0;
+
+  mfem::VectorFunctionCoefficient func(dim, [direction](const mfem::Vector& /*x*/, mfem::Vector& u) {
+    u = 0.0;
+    u[direction] = 1.0;
+  });
+
+  reactionDirections.project(func, mesh->domain("two"));
+
+  return reactionDirections;
+}
+
+double computeContactReactionQoi(BasePhysics& solid_solver, std::shared_ptr<Mesh> mesh)
+{
+  solid_solver.resetStates();
+  solid_solver.advanceTimestep(1.0);
+
+  const FiniteElementDual& reactions = solid_solver.dual("reactions");
+  auto reactionDirections = createReactionDirection(solid_solver, 1, mesh);
+
+  return innerProduct(reactions, reactionDirections);
+}
+
+auto computeContactReactionQoiSensitivities(BasePhysics& solid_solver, std::shared_ptr<Mesh> mesh)
+{
+  EXPECT_EQ(0, solid_solver.cycle());
+
+  double qoi = computeContactReactionQoi(solid_solver, mesh);
+
+  // Verify nonzero reaction forces
+  EXPECT_NE(qoi, 0.0);
+
+  FiniteElementDual shape_sensitivity(solid_solver.shapeDisplacement().space(), "shape sensitivity");
+
+  auto reaction_adjoint_load = createReactionDirection(solid_solver, 1, mesh);
+
+  EXPECT_EQ(1, solid_solver.cycle());
+  solid_solver.setDualAdjointBcs({{"reactions", reaction_adjoint_load}});
+  solid_solver.reverseAdjointTimestep();
+  shape_sensitivity = solid_solver.computeTimestepShapeSensitivity();
+  EXPECT_EQ(0, solid_solver.cycle());
+
+  return std::make_tuple(qoi, shape_sensitivity);
+}
+
+double computeContactReactionQoiAdjustingShape(SolidMechanics<p, dim>& solid_solver,
+                                               const FiniteElementState& shape_derivative_direction, double pertubation,
+                                               std::shared_ptr<Mesh> mesh)
+{
+  FiniteElementState shape_disp(shape_derivative_direction.space(), "input_shape_displacement");
+
+  shape_disp.Add(pertubation, shape_derivative_direction);
+  solid_solver.setShapeDisplacement(shape_disp);
+
+  return computeContactReactionQoi(solid_solver, mesh);
+}
+
 using SolidMechT = smith::SolidMechanicsContact<p, dim>;
 // using SolidMechT = smith::SolidMechanics<p, dim>;
 
@@ -220,6 +282,22 @@ TEST_F(ContactSensitivityFixture, QuasiStaticShapeSensitivities)
   EXPECT_NEAR(directional_deriv, directional_deriv_fd, eps);
 }
 
+TEST_F(ContactSensitivityFixture, ReactionShapeSensitivities)
+{
+  auto solid_solver = createContactSolver(mesh, nonlinear_opts, dyn_opts, mat);
+  auto [qoi_base, shape_sensitivity] = computeContactReactionQoiSensitivities(*solid_solver, mesh);
+
+  solid_solver->resetStates();
+  FiniteElementState derivative_direction(shape_sensitivity.space(), "derivative_direction");
+  fillDirection(*solid_solver, derivative_direction);
+
+  double qoi_plus = computeContactReactionQoiAdjustingShape(*solid_solver, derivative_direction, eps, mesh);
+
+  double directional_deriv = innerProduct(derivative_direction, shape_sensitivity);
+  double directional_deriv_fd = (qoi_plus - qoi_base) / eps;
+  EXPECT_NEAR(directional_deriv, directional_deriv_fd, eps);
+}
+
 }  // namespace smith
 
 int main(int argc, char* argv[])

src/smith/physics/tests/test_solid_weak_form.cpp

@@ -30,7 +30,6 @@
 #include "smith/physics/state/finite_element_state.hpp"
 
 auto element_shape = mfem::Element::QUADRILATERAL;
-
 namespace smith {
 
 struct NeoHookeanWithFieldWithRateForTesting {