diff --git a/src/serac/physics/base_physics.hpp b/src/serac/physics/base_physics.hpp index 20690107a..3d15c57ec 100644 --- a/src/serac/physics/base_physics.hpp +++ b/src/serac/physics/base_physics.hpp @@ -214,7 +214,7 @@ class BasePhysics { * @param dual_name The name of the Finite Element State dual (reaction adjoint load) solution to retrieve * @return The named adjoint Finite Element State */ - virtual const FiniteElementDual& dualAdjoint(const std::string& dual_name) const + virtual const FiniteElementState& dualAdjoint(const std::string& dual_name) const { SLIC_ERROR_ROOT(axom::fmt::format( "dualAdjoint '{}' requested from physics module '{}' which does not support duals", dual_name, name_)); @@ -338,58 +338,6 @@ class BasePhysics { return *shape_displacement_sensitivity_; } - /** - * @brief computes and sets the adjoint load due to reaction load sensitivities - * - * @param dual_name Name for the physics and residual specific dual (reactions) - * @param reaction_direction A FiniteElementState which specifies how the reactions dofs are weighted for the reaction - * qoi - */ - virtual void computeDualAdjointLoad(const std::string& dual_name, const serac::FiniteElementState& reaction_direction) - { - (void)reaction_direction; - SLIC_ERROR_ROOT(axom::fmt::format("computeDualAdjointLoad not enabled in physics module {}, dual name {} requested", - name_, dual_name)); - } - - /** - * @brief computes the partial sensitivity of the dual (reaction) loads in specified direction with respect to - * parameter - * - * @param reaction_direction A FiniteElementState which specifies how the reactions dofs are weighted for the reaction - * qoi - * @param parameter_index the index of the parameter - * @return reaction sensitivity field - * - * @pre `computeDualAdjointLoad' for the desired dual (reaction) and `reverseAdjointTimestep` must be called before - * this - */ - virtual const serac::FiniteElementDual& computeDualSensitivity(const serac::FiniteElementState& reaction_direction, - size_t parameter_index) - { - (void)reaction_direction; - SLIC_ERROR_ROOT(axom::fmt::format("computeDualSensitivity not enabled in physics module {}", name_)); - return *parameters_[parameter_index].sensitivity; - }; - - /** - * @brief computes the partial sensitivity of the reaction loads (in specified direction) with respect to shape - * - * @param reaction_direction A FiniteElementState which specifies how the reactions dofs are weighted for the reaction - * qoi - * @return reaction sensitivity field - * - * @pre `computeDualAdjointLoad' for the desired dual (reaction) and `reverseAdjointTimestep` must be called before - * this - */ - virtual const serac::FiniteElementDual& computeDualShapeSensitivity( - const serac::FiniteElementState& reaction_direction) - { - (void)reaction_direction; - SLIC_ERROR_ROOT(axom::fmt::format("computeDualShapeSensitivity not enabled in physics module {}", name_)); - return *shape_displacement_sensitivity_; - }; - /** * @brief Compute the implicit sensitivity of the quantity of interest with respect to the initial condition fields * @@ -421,6 +369,15 @@ class BasePhysics { SLIC_ERROR_ROOT(axom::fmt::format("Adjoint analysis not defined for physics module {}", name_)); } + /** + * @brief Set the dual loads (dirichlet values) for the adjoint reverse timestep solve + * This must be called after setAdjointLoad + */ + virtual void setDualAdjointBcs(std::unordered_map) + { + SLIC_ERROR_ROOT(axom::fmt::format("Adjoint analysis not defined for physics module {}", name_)); + } + /** * @brief Solve the adjoint reverse timestep problem * @pre It is expected that the forward analysis is complete and the current states are valid @@ -443,12 +400,27 @@ class BasePhysics { * @brief Accessor for getting a single named finite element state primal solution from the physics modules at a given * checkpointed cycle index * - * @param cycle The cycle to retrieve state from * @param state_name The name of the state to retrieve (e.g. "temperature", "displacement") + * @param cycle The cycle to retrieve state from * @return The named primal Finite Element State */ FiniteElementState loadCheckpointedState(const std::string& state_name, int cycle) const; + /** + * @brief Accessor for getting a single named finite element dual solution from the physics modules at a given + * checkpointed cycle index + * + * @param state_name The name of the state to retrieve (e.g. "reaction") + * @param cycle The cycle to retrieve state from + * @return The named Finite Element Dual + */ + virtual FiniteElementDual loadCheckpointedDual([[maybe_unused]] const std::string& state_name, + [[maybe_unused]] int cycle) const + { + SLIC_ERROR_ROOT(axom::fmt::format("loadCheckpointedDual not enabled in physics module {}", name_)); + return *duals_[0]; + } + /** * @brief Get a timestep increment which has been previously checkpointed at the give cycle * @param cycle The previous 'timestep' number where the timestep increment is requested @@ -550,7 +522,7 @@ class BasePhysics { /** * @brief List of adjoint finite element duals associated with this physics module */ - std::vector dual_adjoints_; + std::vector dual_adjoints_; /// @brief The information needed for the physics parameters stored as Finite Element State fields struct ParameterInfo { diff --git a/src/serac/physics/solid_mechanics.hpp b/src/serac/physics/solid_mechanics.hpp index 3ba8818ab..acb9c281f 100644 --- a/src/serac/physics/solid_mechanics.hpp +++ b/src/serac/physics/solid_mechanics.hpp @@ -182,7 +182,7 @@ class SolidMechanics, std::integer_se implicit_sensitivity_displacement_start_of_step_(displacement_.space(), "total_deriv_wrt_displacement."), implicit_sensitivity_velocity_start_of_step_(displacement_.space(), "total_deriv_wrt_velocity."), reactions_(StateManager::newDual(H1{}, detail::addPrefix(physics_name, "reactions"), mesh_tag_)), - reactions_adjoint_load_(reactions_.space(), "reactions_shape_sensitivity"), + reactions_adjoint_bcs_(reactions_.space(), "reactions_shape_sensitivity"), nonlin_solver_(std::move(solver)), ode2_(displacement_.space().TrueVSize(), {.time = time_, .c0 = c0_, .c1 = c1_, .u = u_, .du_dt = v_, .d2u_dt2 = acceleration_}, *nonlin_solver_, @@ -216,7 +216,7 @@ class SolidMechanics, std::integer_se adjoints_.push_back(&adjoint_displacement_); duals_.push_back(&reactions_); - dual_adjoints_.push_back(&reactions_adjoint_load_); + dual_adjoints_.push_back(&reactions_adjoint_bcs_); // Create a pack of the primal field and parameter finite element spaces mfem::ParFiniteElementSpace* test_space = &displacement_.space(); @@ -385,7 +385,8 @@ class SolidMechanics, std::integer_se implicit_sensitivity_displacement_start_of_step_ = 0.0; implicit_sensitivity_velocity_start_of_step_ = 0.0; - reactions_ = 0.0; + reactions_ = 0.0; + reactions_adjoint_bcs_ = 0.0; u_ = 0.0; v_ = 0.0; @@ -753,15 +754,40 @@ class SolidMechanics, std::integer_se } /// @overload - const FiniteElementDual& dualAdjoint(const std::string& dual_name) const override + const FiniteElementState& dualAdjoint(const std::string& dual_name) const override { if (dual_name == "reactions") { - return reactions_adjoint_load_; + return reactions_adjoint_bcs_; } SLIC_ERROR_ROOT(axom::fmt::format( "dualAdjoint '{}' requested from solid mechanics module '{}', but it doesn't exist", dual_name, name_)); - return reactions_adjoint_load_; + return reactions_adjoint_bcs_; + } + + /// @overload + FiniteElementDual loadCheckpointedDual(const std::string& dual_name, int cycle) const override + { + if (dual_name == "reactions") { + auto checkpointed_sol = getCheckpointedStates(cycle); + + FiniteElementDual reactions(reactions_.space(), "reactions_tmp"); + + if (is_quasistatic_) { + reactions = (*residual_)(time_, shape_displacement_, checkpointed_sol.at("displacement"), acceleration_, + *parameters_[parameter_indices].state...); + } else { + reactions = (*residual_)(time_, shape_displacement_, checkpointed_sol.at("displacement"), + checkpointed_sol.at("acceleration"), *parameters_[parameter_indices].state...); + } + + return reactions; + } + + SLIC_ERROR_ROOT( + axom::fmt::format("loadCheckpointedDual '{}' requested from solid mechanics module '{}', but it doesn't exist", + dual_name, name_)); + return reactions_; } /** @@ -1134,8 +1160,10 @@ class SolidMechanics, std::integer_se SERAC_MARK_FUNCTION; auto [r, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(u), acceleration_, *parameters_[parameter_indices].state...); - J_ = assemble(drdu); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); + J_.reset(); + J_ = assemble(drdu); + J_e_.reset(); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); return *J_; }); } @@ -1201,7 +1229,9 @@ class SolidMechanics, std::integer_se std::unique_ptr m_mat(assemble(M)); // J = M + c0 * K + J_.reset(); J_.reset(mfem::Add(1.0, *m_mat, c0_, *k_mat)); + J_e_.reset(); J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); return *J_; @@ -1316,9 +1346,11 @@ class SolidMechanics, std::integer_se SLIC_ERROR_ROOT_IF(disp_adjoint_load == loads.end(), "Adjoint load for \"displacement\" not found."); - displacement_adjoint_load_ = disp_adjoint_load->second; - // Add the sign correction to move the term to the RHS - displacement_adjoint_load_ *= -1.0; + if (disp_adjoint_load != loads.end()) { + displacement_adjoint_load_ = disp_adjoint_load->second; + // Add the sign correction to move the term to the RHS + displacement_adjoint_load_ *= -1.0; + } auto velo_adjoint_load = loads.find("velocity"); @@ -1337,15 +1369,24 @@ class SolidMechanics, std::integer_se } } + virtual void setDualAdjointBcs(std::unordered_map bcs) override + { + SLIC_ERROR_ROOT_IF(bcs.size() == 0, "Adjoint load container size must be greater than 0 in the solid mechanics."); + + auto reaction_adjoint_load = bcs.find("reactions"); + + SLIC_ERROR_ROOT_IF(reaction_adjoint_load == bcs.end(), "Adjoint load for \"reaction\" not found."); + + if (reaction_adjoint_load != bcs.end()) { + reactions_adjoint_bcs_ = reaction_adjoint_load->second; + } + } + /// @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"); @@ -1361,11 +1402,20 @@ class SolidMechanics, std::integer_se 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(jacobian->Transpose()); + 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(); - for (const auto& bc : bcs_.essentials()) { - bc.apply(*J_T, displacement_adjoint_load_, adjoint_essential); + 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); @@ -1396,7 +1446,7 @@ class SolidMechanics, std::integer_se 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_, reactions_adjoint_bcs_, bcs_, lin_solver); } time_end_step_ = time_; @@ -1462,47 +1512,6 @@ class SolidMechanics, std::integer_se /// @brief getter for nodal forces (before zeroing-out essential dofs) const serac::FiniteElementDual& reactions() const { return reactions_; }; - /// @overload - void computeDualAdjointLoad(const std::string& dual_name, - const serac::FiniteElementState& reaction_direction) override - { - SLIC_ERROR_ROOT_IF(dual_name != "reactions", "Solid mechanics has reactions as its only dual"); - - auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, - *parameters_[parameter_indices].state...); - std::unique_ptr jacobian = assemble(drdu); - reactions_adjoint_load_ = 0.0; - jacobian->MultTranspose(reaction_direction, reactions_adjoint_load_); - setAdjointLoad({{"displacement", reactions_adjoint_load_}}); - } - - /// @overload - const serac::FiniteElementDual& computeDualSensitivity(const serac::FiniteElementState& reaction_direction, - size_t parameter_field) override - { - SLIC_ASSERT_MSG(parameter_field < sizeof...(parameter_indices), - axom::fmt::format("Invalid parameter index '{}' requested for reaction sensitivity.")); - - auto drdparam = serac::get(d_residual_d_[parameter_field](time_end_step_)); - auto drdparam_mat = assemble(drdparam); - - drdparam_mat->MultTranspose(reaction_direction, *parameters_[parameter_field].sensitivity); - - return *parameters_[parameter_field].sensitivity; - }; - - /// @overload - const serac::FiniteElementDual& computeDualShapeSensitivity( - const serac::FiniteElementState& reaction_direction) override - { - auto drdshape = - serac::get((*residual_)(time_end_step_, differentiate_wrt(shape_displacement_), displacement_, - acceleration_, *parameters_[parameter_indices].state...)); - auto drdshape_mat = assemble(drdshape); - drdshape_mat->MultTranspose(reaction_direction, *shape_displacement_sensitivity_); - return *shape_displacement_sensitivity_; - }; - protected: /// The compile-time finite element trial space for displacement and velocity (H1 of order p) using trial = H1; @@ -1546,7 +1555,7 @@ class SolidMechanics, std::integer_se FiniteElementDual reactions_; /// sensitivity of qoi with respect to reaction forces - FiniteElementDual reactions_adjoint_load_; + FiniteElementState reactions_adjoint_bcs_; /// serac::Functional that is used to calculate the residual and its derivatives std::unique_ptr> residual_; @@ -1748,8 +1757,10 @@ class SolidMechanics, std::integer_se // use the most recently evaluated Jacobian auto [_, drdu] = (*residual_)(time_, shape_displacement_, differentiate_wrt(displacement_), acceleration_, *parameters_[parameter_indices].previous_state...); - J_ = assemble(drdu); - J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); + J_.reset(); + J_ = assemble(drdu); + J_e_.reset(); + J_e_ = bcs_.eliminateAllEssentialDofsFromMatrix(*J_); r *= -1.0; diff --git a/src/serac/physics/tests/dynamic_solid_adjoint.cpp b/src/serac/physics/tests/dynamic_solid_adjoint.cpp index fdb2ecda1..910ca3380 100644 --- a/src/serac/physics/tests/dynamic_solid_adjoint.cpp +++ b/src/serac/physics/tests/dynamic_solid_adjoint.cpp @@ -17,6 +17,7 @@ #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 { @@ -43,20 +44,24 @@ constexpr double boundary_disp = 0.013; constexpr double initial_interior_disp = 0.03; constexpr double initial_interior_velo = 0.04; -// MRT: add explicit velocity dependence -double computeStepQoi(const FiniteElementState& displacement, double dt) +double computeStepQoi(const FiniteElementState& displacement, const FiniteElementDual& reactions, double dt) { FiniteElementState displacement_error(displacement); displacement_error = disp_target; displacement_error.Add(1.0, displacement); - return 0.5 * dt * innerProduct(displacement_error, displacement_error); + return 0.5 * dt * innerProduct(displacement_error, displacement_error) + + 0.05 * dt * innerProduct(reactions, reactions); } -void computeStepAdjointLoad(const FiniteElementState& displacement, FiniteElementDual& d_qoi_d_displacement, double dt) +void computeStepAdjointLoad(const FiniteElementState& displacement, const FiniteElementDual& reactions, + FiniteElementDual& d_qoi_d_displacement, FiniteElementState& d_qoi_d_reactions, double dt) { d_qoi_d_displacement = disp_target; d_qoi_d_displacement.Add(1.0, displacement); d_qoi_d_displacement *= dt; + + d_qoi_d_reactions = reactions; + d_qoi_d_reactions *= 0.1 * dt; } void applyInitialAndBoundaryConditions(SolidMechanics& solid_solver) @@ -110,14 +115,18 @@ double computeSolidMechanicsQoi(BasePhysics& solid_solver, const TimeSteppingInf auto dts = ts_info.dts; solid_solver.advanceTimestep(dts(0)); // advance by 0.0 seconds to get initial acceleration solid_solver.outputStateToDisk(); - FiniteElementState dispForObjective = solid_solver.state("displacement"); - double qoi = computeStepQoi(dispForObjective, 0.5 * (dts(0) + dts(1))); + FiniteElementState dispForObjective = solid_solver.state("displacement"); + FiniteElementDual reactionsForObjective = solid_solver.dual("reactions"); + double qoi = computeStepQoi(dispForObjective, reactionsForObjective, 0.5 * (dts(0) + dts(1))); + for (int i = 1; i <= ts_info.numTimesteps(); ++i) { solid_solver.advanceTimestep(dts(i)); solid_solver.outputStateToDisk(); - dispForObjective = solid_solver.state("displacement"); - qoi += computeStepQoi(dispForObjective, 0.5 * (dts(i) + dts(i + 1))); + + dispForObjective = solid_solver.state("displacement"); + reactionsForObjective = solid_solver.dual("reactions"); + qoi += computeStepQoi(dispForObjective, reactionsForObjective, 0.5 * (dts(i) + dts(i + 1))); } return qoi; } @@ -137,18 +146,22 @@ std::tuple comp FiniteElementDual shape_sensitivity(solid_solver.shapeDisplacement().space(), "shape sensitivity"); shape_sensitivity = 0.0; - FiniteElementDual adjoint_load(solid_solver.state("displacement").space(), "adjoint_displacement_load"); + FiniteElementDual adjoint_load(solid_solver.state("displacement").space(), "adjoint_displacement_load"); + FiniteElementState adjoint_bcs(solid_solver.dual("reactions").space(), "adjoint_reaction_bcs"); // for solids, we go back to time = 0, because there is an extra hidden implicit solve at the start // consider unifying the interface between solids and thermal for (int i = solid_solver.cycle(); i > 0; --i) { auto previous_displacement = solid_solver.loadCheckpointedState("displacement", solid_solver.cycle()); + auto previous_reactions = solid_solver.loadCheckpointedDual("reactions", solid_solver.cycle()); computeStepAdjointLoad( - previous_displacement, adjoint_load, + previous_displacement, previous_reactions, adjoint_load, adjoint_bcs, 0.5 * (solid_solver.getCheckpointedTimestep(i - 1) + solid_solver.getCheckpointedTimestep(i))); EXPECT_EQ(i, solid_solver.cycle()); solid_solver.setAdjointLoad({{"displacement", adjoint_load}}); + solid_solver.setDualAdjointBcs({{"reactions", adjoint_bcs}}); solid_solver.reverseAdjointTimestep(); + shape_sensitivity += solid_solver.computeTimestepShapeSensitivity(); EXPECT_EQ(i - 1, solid_solver.cycle()); } @@ -336,12 +349,8 @@ TEST_F(SolidMechanicsSensitivityFixture, WhenShapeSensitivitiesCalledTwice_GetSa int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); - MPI_Init(&argc, &argv); - - axom::slic::SimpleLogger logger; - std::cout << std::setprecision(16); + serac::initialize(argc, argv); int result = RUN_ALL_TESTS(); - MPI_Finalize(); - + serac::exitGracefully(result); return result; } diff --git a/src/serac/physics/tests/quasistatic_solid_adjoint.cpp b/src/serac/physics/tests/quasistatic_solid_adjoint.cpp index bf13b2868..90f88b00d 100644 --- a/src/serac/physics/tests/quasistatic_solid_adjoint.cpp +++ b/src/serac/physics/tests/quasistatic_solid_adjoint.cpp @@ -17,6 +17,7 @@ #include "serac/mesh/mesh_utils.hpp" #include "serac/physics/state/state_manager.hpp" #include "serac/serac_config.hpp" +#include "serac/infrastructure/terminator.hpp" struct ParameterizedLinearIsotropicSolid { using State = ::serac::Empty; ///< this material has no internal variables @@ -243,12 +244,9 @@ TEST(quasistatic, finiteDifference) int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); - MPI_Init(&argc, &argv); - - axom::slic::SimpleLogger logger; - std::cout << std::setprecision(16); + serac::initialize(argc, argv); int result = RUN_ALL_TESTS(); - MPI_Finalize(); + serac::exitGracefully(result); return result; } diff --git a/src/serac/physics/tests/solid_reaction_adjoint.cpp b/src/serac/physics/tests/solid_reaction_adjoint.cpp index 2c03485d9..b154be211 100644 --- a/src/serac/physics/tests/solid_reaction_adjoint.cpp +++ b/src/serac/physics/tests/solid_reaction_adjoint.cpp @@ -18,6 +18,7 @@ #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 { @@ -30,7 +31,8 @@ const std::string mesh_tag = "mesh"; const std::string physics_prefix = "solid"; using SolidMaterial = solid_mechanics::ParameterizedNeoHookeanSolid; -auto geoNonlinear = GeometricNonlinearities::On; +// using SolidMaterial = solid_mechanics::ParameterizedLinearIsotropicSolid; +auto geoNonlinear = GeometricNonlinearities::Off; constexpr double boundary_disp = 0.013; constexpr double shear_modulus_value = 1.0; @@ -101,12 +103,15 @@ FiniteElementState createReactionDirection(const BasePhysics& solid_solver, int double computeSolidMechanicsQoi(BasePhysics& solid_solver) { - solid_solver.advanceTimestep(0.0); + solid_solver.resetStates(); + solid_solver.advanceTimestep(0.1); - const FiniteElementDual& reactions = solid_solver.dual("reactions"); + const FiniteElementDual& reactions = solid_solver.dual("reactions"); + auto reactionDirections = createReactionDirection(solid_solver, 0); + // reactionDirections = solid_solver.dual("reactions"); - auto reactionDirections = createReactionDirection(solid_solver, 0); - return innerProduct(reactions, reactionDirections); + const FiniteElementState& displacements = solid_solver.state("displacement"); + return innerProduct(reactions, reactionDirections) + 0.05 * innerProduct(displacements, displacements); } auto computeSolidMechanicsQoiSensitivities(BasePhysics& solid_solver) @@ -119,16 +124,18 @@ auto computeSolidMechanicsQoiSensitivities(BasePhysics& solid_solver) FiniteElementDual shear_modulus_sensitivity(StateManager::mesh(mesh_tag), H1

{}, "shear modulus sensitivity"); shear_modulus_sensitivity = 0.0; - auto reactionDirections = createReactionDirection(solid_solver, 0); + auto reaction_adjoint_load = createReactionDirection(solid_solver, 0); + + FiniteElementDual displacement_adjoint_load(solid_solver.state("displacement").space(), "displacement_adjoint_load"); + displacement_adjoint_load = solid_solver.state("displacement"); + displacement_adjoint_load *= 0.1; - solid_solver.computeDualAdjointLoad(solid_solver.dualNames()[0], reactionDirections); + solid_solver.setAdjointLoad({{"displacement", displacement_adjoint_load}}); + solid_solver.setDualAdjointBcs({{"reactions", reaction_adjoint_load}}); solid_solver.reverseAdjointTimestep(); shear_modulus_sensitivity += solid_solver.computeTimestepSensitivity(0); - shear_modulus_sensitivity += solid_solver.computeDualSensitivity(reactionDirections, 0); - shape_sensitivity += solid_solver.computeTimestepShapeSensitivity(); - shape_sensitivity += solid_solver.computeDualShapeSensitivity(reactionDirections); return std::make_tuple(qoi, shear_modulus_sensitivity, shape_sensitivity); } @@ -177,7 +184,6 @@ struct SolidMechanicsSensitivityFixture : public ::testing::Test { void fillDirection(FiniteElementState& direction) const { auto sz = direction.Size(); - std::cout << "sizes = " << sz << std::endl; for (int i = 0; i < sz; ++i) { direction(i) = -1.2 + 2.02 * (double(i) / sz); } @@ -207,7 +213,6 @@ TEST_F(SolidMechanicsSensitivityFixture, ReactionShapeSensitivities) double qoi_plus = computeSolidMechanicsQoiAdjustingShape(*solid_solver, derivative_direction, eps); double directional_deriv = innerProduct(derivative_direction, shape_sensitivity); - EXPECT_NEAR(qoi_base, -0.35, 1e-14); EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); } @@ -223,7 +228,6 @@ TEST_F(SolidMechanicsSensitivityFixture, ReactionParameterSensitivities) double qoi_plus = computeSolidMechanicsQoiAdjustingShearModulus(*solid_solver, derivative_direction, eps); double directional_deriv = innerProduct(derivative_direction, shear_modulus_sensitivity); - EXPECT_NEAR(qoi_base, -0.35, 1e-14); EXPECT_NEAR(directional_deriv, (qoi_plus - qoi_base) / eps, eps); } @@ -232,12 +236,8 @@ TEST_F(SolidMechanicsSensitivityFixture, ReactionParameterSensitivities) int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); - MPI_Init(&argc, &argv); - - axom::slic::SimpleLogger logger; - std::cout << std::setprecision(16); + serac::initialize(argc, argv); int result = RUN_ALL_TESTS(); - MPI_Finalize(); - + serac::exitGracefully(result); return result; }