diff --git a/data/meshes/contact_two_blocks.g b/data/meshes/contact_two_blocks.g new file mode 100644 index 000000000..1355b6306 Binary files /dev/null and b/data/meshes/contact_two_blocks.g differ diff --git a/src/serac/physics/contact/contact_data.cpp b/src/serac/physics/contact/contact_data.cpp index c7f595bd0..05fc38d9f 100644 --- a/src/serac/physics/contact/contact_data.cpp +++ b/src/serac/physics/contact/contact_data.cpp @@ -41,6 +41,15 @@ void ContactData::addContactInteraction(int interaction_id, const std::set& } } +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 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 ContactData::mergedJacobian() const return std::make_unique(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 ContactData::jacobianFunction(mfem::HypreParMatrix* orig_J) const { @@ -390,7 +407,10 @@ std::unique_ptr 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 diff --git a/src/serac/physics/contact/contact_data.hpp b/src/serac/physics/contact/contact_data.hpp index 5f9a83714..e2ad3dcb8 100644 --- a/src/serac/physics/contact/contact_data.hpp +++ b/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? diff --git a/src/serac/physics/solid_mechanics.cpp b/src/serac/physics/solid_mechanics.cpp index ac1b34df7..23038d99d 100644 --- a/src/serac/physics/solid_mechanics.cpp +++ b/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); diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index b6b9b7c7a..f4a4f36e1 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/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, 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, 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, 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(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, std::integer_se *parameters_[parameter_indices].state...)); std::unique_ptr 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, 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(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 diff --git a/src/serac/physics/solid_mechanics_contact.hpp b/src/serac/physics/solid_mechanics_contact.hpp index a85547389..7b47aeeb7 100644 --- a/src/serac/physics/solid_mechanics_contact.hpp +++ b/src/serac/physics/solid_mechanics_contact.hpp @@ -110,6 +110,17 @@ class SolidMechanicsContact, duals_.push_back(&forces_); } + /// @overload + void resetStates(int cycle = 0, double time = 0.0) override + { + SolidMechanicsBase::resetStates(cycle, time); + forces_ = 0.0; + contact_.setDisplacements(shape_displacement_, displacement_); + contact_.reset(); + double dt = 0.0; + contact_.update(cycle, time, dt); + } + /// @brief Build the quasi-static operator corresponding to the total Lagrangian formulation std::unique_ptr buildQuasistaticOperator() override { @@ -123,7 +134,8 @@ class SolidMechanicsContact, // 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 @@ -177,7 +189,8 @@ class SolidMechanicsContact, *parameters_[parameter_indices].state...); J_ = assemble(drdu); - // get 11-block holding jacobian contributions + // contact_.setDisplacements(u_shape, u_blk); + // get 11-block holding jacobian contributions auto block_J = contact_.jacobianFunction(J_.release()); block_J->owns_blocks = false; J_ = std::unique_ptr(static_cast(&block_J->GetBlock(0, 0))); @@ -295,13 +308,29 @@ class SolidMechanicsContact, } if (use_warm_start_) { + // Update the system residual + mfem::Vector augmented_residual(displacement_.Size() + contact_.numPressureDofs()); + augmented_residual = 0.0; + const mfem::Vector res = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, + *parameters_[parameter_indices].state...); + + contact_.setPressures(mfem::Vector(augmented_residual, displacement_.Size(), contact_.numPressureDofs())); + contact_.update(cycle_, time_, dt); + // TODO this copy is required as the sundials solvers do not allow move assignments because of their memory + // tracking strategy + // See https://github.com/mfem/mfem/issues/3531 + mfem::Vector r_blk(augmented_residual, 0, displacement_.space().TrueVSize()); + r_blk = res; + + contact_.residualFunction(shape_displacement_, displacement_, augmented_residual); + r_blk.SetSubVector(bcs_.allEssentialTrueDofs(), 0.0); + // use the most recently evaluated Jacobian auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, *parameters_[parameter_indices].previous_state...); J_.reset(); J_ = assemble(drdu); - contact_.update(cycle_, time_, dt); if (contact_.haveLagrangeMultipliers()) { J_offsets_ = mfem::Array({0, displacement_.Size(), displacement_.Size() + contact_.numPressureDofs()}); J_constraint_ = contact_.jacobianFunction(J_.release()); @@ -333,29 +362,16 @@ class SolidMechanicsContact, J_operator_ = J_.get(); } - // Update the linearized Jacobian matrix - mfem::Vector augmented_residual(displacement_.space().TrueVSize() + contact_.numPressureDofs()); - augmented_residual = 0.0; - const mfem::Vector res = (*residual_)(time_ + dt, shape_displacement_, displacement_, acceleration_, - *parameters_[parameter_indices].state...); - - // TODO this copy is required as the sundials solvers do not allow move assignments because of their memory - // tracking strategy - // See https://github.com/mfem/mfem/issues/3531 - mfem::Vector r(augmented_residual, 0, displacement_.space().TrueVSize()); - r = res; - r += contact_.forces(); - augmented_residual *= -1.0; mfem::Vector augmented_solution(displacement_.space().TrueVSize() + contact_.numPressureDofs()); augmented_solution = 0.0; mfem::Vector du(augmented_solution, 0, displacement_.space().TrueVSize()); du = du_; - mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du, r); + mfem::EliminateBC(*J_, *J_e_, constrained_dofs, du, r_blk); for (int i = 0; i < constrained_dofs.Size(); i++) { - int j = constrained_dofs[i]; - r[j] = du[j]; + int j = constrained_dofs[i]; + r_blk[j] = du[j]; } auto& lin_solver = nonlin_solver_->linearSolver(); @@ -370,6 +386,52 @@ class SolidMechanicsContact, displacement_ += du_; } + /// @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(static_cast(&block_J->GetBlock(0, 0))); + auto J_T = std::unique_ptr(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((*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(static_cast(&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_; @@ -378,19 +440,24 @@ class SolidMechanicsContact, 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::residual_; using SolidMechanicsBase::residual_with_bcs_; + using SolidMechanicsBase::time_end_step_; using SolidMechanicsBase::use_warm_start_; + using SolidMechanicsBase::warmStartDisplacement; /// Pointer to the Jacobian operator (J_ if no Lagrange multiplier contact, J_constraint_ otherwise) mfem::Operator* J_operator_; diff --git a/src/serac/physics/tests/CMakeLists.txt b/src/serac/physics/tests/CMakeLists.txt index f7505e492..01b6f0437 100644 --- a/src/serac/physics/tests/CMakeLists.txt +++ b/src/serac/physics/tests/CMakeLists.txt @@ -21,6 +21,11 @@ set(physics_serial_test_sources solid_multi_material.cpp ) +blt_list_append(TO physics_serial_test_sources + ELEMENTS + contact_solid_adjoint.cpp + IF TRIBOL_FOUND) + serac_add_tests(SOURCES ${physics_serial_test_sources} DEPENDS_ON ${physics_test_depends} NUM_MPI_TASKS 1) diff --git a/src/serac/physics/tests/contact_beam.cpp b/src/serac/physics/tests/contact_beam.cpp index daf071884..3ccc18cf5 100644 --- a/src/serac/physics/tests/contact_beam.cpp +++ b/src/serac/physics/tests/contact_beam.cpp @@ -42,7 +42,7 @@ TEST_P(ContactTest, beam) auto mesh = mesh::refineAndDistribute(buildMeshFromFile(filename), 1, 0); auto& pmesh = serac::StateManager::setMesh(std::move(mesh), "beam_mesh"); - LinearSolverOptions linear_options{.linear_solver = LinearSolver::Strumpack, .print_level = 1}; + LinearSolverOptions linear_options{.linear_solver = LinearSolver::Strumpack, .print_level = 0}; #ifndef MFEM_USE_STRUMPACK SLIC_INFO_ROOT("Contact requires MFEM built with strumpack."); return; diff --git a/src/serac/physics/tests/contact_patch_tied.cpp b/src/serac/physics/tests/contact_patch_tied.cpp index 0eb767cf8..fbee17b8d 100644 --- a/src/serac/physics/tests/contact_patch_tied.cpp +++ b/src/serac/physics/tests/contact_patch_tied.cpp @@ -55,7 +55,7 @@ TEST_P(ContactPatchTied, patch) // }; // #elif defined(MFEM_USE_STRUMPACK) #ifdef MFEM_USE_STRUMPACK - LinearSolverOptions linear_options{.linear_solver = LinearSolver::Strumpack, .print_level = 1}; + LinearSolverOptions linear_options{.linear_solver = LinearSolver::Strumpack, .print_level = 0}; #else LinearSolverOptions linear_options{}; SLIC_INFO_ROOT("Contact requires MFEM built with strumpack."); @@ -64,8 +64,8 @@ TEST_P(ContactPatchTied, patch) NonlinearSolverOptions nonlinear_options{.nonlin_solver = NonlinearSolver::Newton, .relative_tol = 1.0e-10, - .absolute_tol = 1.0e-10, - .max_iterations = 20, + .absolute_tol = 5.0e-10, + .max_iterations = 25, .print_level = 1}; ContactOptions contact_options{.method = ContactMethod::SingleMortar, diff --git a/src/serac/physics/tests/contact_solid_adjoint.cpp b/src/serac/physics/tests/contact_solid_adjoint.cpp new file mode 100644 index 000000000..6f9559abf --- /dev/null +++ b/src/serac/physics/tests/contact_solid_adjoint.cpp @@ -0,0 +1,213 @@ +// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) +#include +#include +#include + +#include "serac/physics/solid_mechanics_contact.hpp" +#include "serac/physics/materials/solid_material.hpp" + +#include "axom/slic/core/SimpleLogger.hpp" +#include +#include "mfem.hpp" + +#include "serac/mesh/mesh_utils.hpp" +#include "serac/physics/state/state_manager.hpp" +#include "serac/serac_config.hpp" +#include "serac/infrastructure/terminator.hpp" + +namespace serac { + +constexpr int dim = 3; +constexpr int p = 1; + +const std::string mesh_tag = "mesh"; +const std::string physics_prefix = "solid"; + +using SolidMaterial = solid_mechanics::NeoHookean; + +struct TimeSteppingInfo { + TimeSteppingInfo() : dts({0.0, 0.2, 0.4, 0.24, 0.12, 0.0}) {} + + int numTimesteps() const { return dts.Size() - 2; } + + mfem::Vector dts; +}; + +constexpr double disp_target = -0.34; + +// MRT: add explicit velocity dependence +double computeStepQoi(const FiniteElementState& displacement) +{ + FiniteElementState displacement_error(displacement); + displacement_error = -disp_target; + displacement_error.Add(1.0, displacement); + return 0.5 * innerProduct(displacement_error, displacement_error); +} + +void computeStepAdjointLoad(const FiniteElementState& displacement, FiniteElementDual& d_qoi_d_displacement) +{ + d_qoi_d_displacement = -disp_target; + d_qoi_d_displacement.Add(1.0, displacement); +} + +using SolidMechT = serac::SolidMechanicsContact; + +std::unique_ptr createContactSolver(const NonlinearSolverOptions& nonlinear_opts, + const TimesteppingOptions& dyn_opts, const SolidMaterial& mat) +{ + static int iter = 0; + + auto solid = std::make_unique(nonlinear_opts, solid_mechanics::direct_linear_options, dyn_opts, + physics_prefix + std::to_string(iter++), mesh_tag, + std::vector{}, 0, 0.0, false, true); + + Domain whole_mesh = EntireDomain(StateManager::mesh(mesh_tag)); + solid->setMaterial(mat, whole_mesh); + + solid->setDisplacementBCs({2}, [](const mfem::Vector& /*X*/, double /*t*/, mfem::Vector& disp) { disp = 0.0; }); + solid->setDisplacementBCs({4}, [](const mfem::Vector& /*X*/, double /*t*/, mfem::Vector& disp) { + disp = 0.0; + disp[1] = -0.1; + }); + + auto contact_type = serac::ContactEnforcement::Penalty; + double element_length = 1.0; + double penalty = 105.1 * mat.K / element_length; + + serac::ContactOptions contact_options{.method = serac::ContactMethod::SingleMortar, + .enforcement = contact_type, + .type = serac::ContactType::Frictionless, + .penalty = penalty}; + auto contact_interaction_id = 0; + solid->addContactInteraction(contact_interaction_id, {3}, {5}, contact_options); + + solid->completeSetup(); + + return solid; +} + +double computeSolidMechanicsQoi(BasePhysics& solid_solver, const TimeSteppingInfo& ts_info) +{ + auto dts = ts_info.dts; + solid_solver.resetStates(); + solid_solver.outputStateToDisk("paraview_contact"); + solid_solver.advanceTimestep(1.0); + solid_solver.outputStateToDisk("paraview_contact"); + return computeStepQoi(solid_solver.state("displacement")); +} + +auto computeContactQoiSensitivities(BasePhysics& solid_solver, const TimeSteppingInfo& ts_info) +{ + EXPECT_EQ(0, solid_solver.cycle()); + + double qoi = computeSolidMechanicsQoi(solid_solver, ts_info); + + FiniteElementDual shape_sensitivity(solid_solver.shapeDisplacement().space(), "shape sensitivity"); + FiniteElementDual adjoint_load(solid_solver.state("displacement").space(), "adjoint_displacement_load"); + + auto previous_displacement = solid_solver.loadCheckpointedState("displacement", solid_solver.cycle()); + computeStepAdjointLoad(previous_displacement, adjoint_load); + EXPECT_EQ(1, solid_solver.cycle()); + solid_solver.setAdjointLoad({{"displacement", adjoint_load}}); + solid_solver.reverseAdjointTimestep(); + shape_sensitivity = solid_solver.computeTimestepShapeSensitivity(); + EXPECT_EQ(0, solid_solver.cycle()); + + return std::make_tuple(qoi, shape_sensitivity); +} + +double computeSolidMechanicsQoiAdjustingShape(SolidMechanics& solid_solver, const TimeSteppingInfo& ts_info, + const FiniteElementState& shape_derivative_direction, double pertubation) +{ + FiniteElementState shape_disp(shape_derivative_direction.space(), "input_shape_displacement"); + + shape_disp.Add(pertubation, shape_derivative_direction); + solid_solver.setShapeDisplacement(shape_disp); + + return computeSolidMechanicsQoi(solid_solver, ts_info); +} + +struct ContactSensitivityFixture : public ::testing::Test { + void SetUp() override + { + MPI_Barrier(MPI_COMM_WORLD); + StateManager::initialize(dataStore, "contact_solve"); + std::string filename = std::string(SERAC_REPO_DIR) + "/data/meshes/contact_two_blocks.g"; + mesh = &StateManager::setMesh(mesh::refineAndDistribute(buildMeshFromFile(filename), 0), mesh_tag); + + mat.density = 1.0; + mat.K = 1.0; + mat.G = 0.1; + } + + void fillDirection(FiniteElementState& direction) const + { + auto sz = direction.Size(); + for (int i = 0; i < sz; ++i) { + direction(i) = -1.2 + 2.02 * (double(i) / sz); + } + } + + axom::sidre::DataStore dataStore; + mfem::ParMesh* mesh; + + NonlinearSolverOptions nonlinear_opts{.relative_tol = 1.0e-10, .absolute_tol = 1.0e-12}; + + bool dispBc = true; + TimesteppingOptions dyn_opts{.timestepper = TimestepMethod::QuasiStatic}; + + SolidMaterial mat; + TimeSteppingInfo tsInfo; + + static constexpr double eps = 2e-7; +}; + +TEST_F(ContactSensitivityFixture, WhenShapeSensitivitiesCalledTwice_GetSameObjectiveAndGradient) +{ + auto solid_solver = createContactSolver(nonlinear_opts, dyn_opts, mat); + auto [qoi1, shape_sensitivity1] = computeContactQoiSensitivities(*solid_solver, tsInfo); + auto [qoi2, shape_sensitivity2] = computeContactQoiSensitivities(*solid_solver, tsInfo); + + EXPECT_EQ(qoi1, qoi2); + + solid_solver->resetStates(); + FiniteElementState derivative_direction(shape_sensitivity1.space(), "derivative_direction"); + fillDirection(derivative_direction); + + double directional_deriv1 = innerProduct(derivative_direction, shape_sensitivity1); + double directional_deriv2 = innerProduct(derivative_direction, shape_sensitivity2); + EXPECT_EQ(directional_deriv1, directional_deriv2); +} + +TEST_F(ContactSensitivityFixture, QuasiStaticShapeSensitivities) +{ + auto solid_solver = createContactSolver(nonlinear_opts, dyn_opts, mat); + auto [qoi_base, shape_sensitivity] = computeContactQoiSensitivities(*solid_solver, tsInfo); + + solid_solver->resetStates(); + FiniteElementState derivative_direction(shape_sensitivity.space(), "derivative_direction"); + fillDirection(derivative_direction); + + double qoi_plus = computeSolidMechanicsQoiAdjustingShape(*solid_solver, tsInfo, derivative_direction, eps); + + double directional_deriv = innerProduct(derivative_direction, shape_sensitivity); + double directional_deriv_fd = (qoi_plus - qoi_base) / eps; + // std::cout << "qoi, derivs = " << qoi_base << " " << directional_deriv << " " << directional_deriv_fd << std::endl; + EXPECT_NEAR(directional_deriv, directional_deriv_fd, 0.2); // These are very large tolerances +} + +} // namespace serac + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + serac::initialize(argc, argv); + int result = RUN_ALL_TESTS(); + serac::exitGracefully(result); + + return result; +}